Graphviz  2.41.20170921.2350
constrained_majorization_ipsep.c
Go to the documentation of this file.
1 /* $Id$ $Revision$ */
2 /* vim:set shiftwidth=4 ts=8: */
3 
20 /**********************************************************
21  * Based on constrained_majorization.c
22  *
23  * Perform stress majorization subject
24  * to separation constraints, for background see the paper:
25  * "IPSep-CoLa: An Incremental Procedure for Separation Constraint Layout of Graphs"
26  * by Tim Dwyer, Yehuda Koren and Kim Marriott
27  *
28  * Available separation constraints so far are:
29  * o Directed edge constraints
30  * o Node non-overlap constraints
31  * o Cluster containment constraints
32  * o Cluster/node non-overlap constraints
33  *
34  * Tim Dwyer, 2006
35  **********************************************************/
36 
37 #include "digcola.h"
38 #ifdef IPSEPCOLA
39 #include <math.h>
40 #include <stdlib.h>
41 #include <time.h>
42 #include <stdio.h>
43 #include <float.h>
44 #include "stress.h"
45 #include "dijkstra.h"
46 #include "bfs.h"
47 #include "matrix_ops.h"
48 #include "kkutils.h"
49 #include "conjgrad.h"
50 #include <csolve_VPSC.h>
51 #include "quad_prog_vpsc.h"
52 #include "quad_prog_solver.h"
53 #include "matrix_ops.h"
54 
55 #define localConstrMajorIterations 1000
56 
57 int stress_majorization_cola(vtx_data * graph, /* Input graph in sparse representation */
58  int n, /* Number of nodes */
59  int nedges_graph, /* Number of edges */
60  double **d_coords, /* Coordinates of nodes (output layout) */
61  node_t ** nodes, /* Original nodes */
62  int dim, /* Dimemsionality of layout */
63  int model, /* difference model */
64  int maxi, /* max iterations */
65  ipsep_options * opt)
66 {
67  int iterations = 0; /* Output: number of iteration of the process */
68 
69  /*************************************************
70  ** Computation of full, dense, unrestricted k-D **
71  ** stress minimization by majorization **
72  ** This function imposes HIERARCHY CONSTRAINTS **
73  *************************************************/
74 
75  int i, j, k;
76  float *lap1 = NULL;
77  float *dist_accumulator = NULL;
78  float *tmp_coords = NULL;
79  float **b = NULL;
80  double *degrees = NULL;
81  float *lap2 = NULL;
82  int lap_length;
83  float *f_storage = NULL;
84  float **coords = NULL;
85  int orig_n = n;
86 
87  /*double conj_tol=tolerance_cg; *//* tolerance of Conjugate Gradient */
88  CMajEnvVPSC *cMajEnvHor = NULL;
89  CMajEnvVPSC *cMajEnvVrt = NULL;
90  double y_0;
91  int length;
92  DistType diameter;
93  float *Dij = NULL;
94  float constant_term;
95  int count;
96  double degree;
97  int step;
98  float val;
99  double old_stress, new_stress = 0;
100  boolean converged;
101  int len;
102  double nsizeScale = 0;
103  float maxEdgeLen = 0;
104  double max = 1;
105 
106  initLayout(graph, n, dim, d_coords, nodes);
107  if (n == 1)
108  return 0;
109 
110  for (i = 0; i < n; i++) {
111  for (j = 1; j < graph[i].nedges; j++) {
112  maxEdgeLen = MAX(graph[i].ewgts[j], maxEdgeLen);
113  }
114  }
115 
116  /****************************************************
117  ** Compute the all-pairs-shortest-distances matrix **
118  ****************************************************/
119 
120  if (maxi == 0)
121  return iterations;
122 
123  if (Verbose)
124  start_timer();
125 
126  if (model == MODEL_SUBSET) {
127  /* weight graph to separate high-degree nodes */
128  /* and perform slower Dijkstra-based computation */
129  if (Verbose)
130  fprintf(stderr, "Calculating subset model");
132  } else if (model == MODEL_CIRCUIT) {
133  Dij = circuitModel(graph, n);
134  if (!Dij) {
135  agerr(AGWARN,
136  "graph is disconnected. Hence, the circuit model\n");
137  agerr(AGPREV,
138  "is undefined. Reverting to the shortest path model.\n");
139  }
140  } else if (model == MODEL_MDS) {
141  if (Verbose)
142  fprintf(stderr, "Calculating MDS model");
143  Dij = mdsModel(graph, n);
144  }
145  if (!Dij) {
146  if (Verbose)
147  fprintf(stderr, "Calculating shortest paths");
148  Dij = compute_apsp_packed(graph, n);
149  }
150  if (Verbose) {
151  fprintf(stderr, ": %.2f sec\n", elapsed_sec());
152  fprintf(stderr, "Setting initial positions");
153  start_timer();
154  }
155 
156  diameter = -1;
157  length = n + n * (n - 1) / 2;
158  for (i = 0; i < length; i++) {
159  if (Dij[i] > diameter) {
160  diameter = (int) Dij[i];
161  }
162  }
163 
164  /* for numerical stability, scale down layout */
165  /* No Jiggling, might conflict with constraints */
166  for (i = 0; i < dim; i++) {
167  for (j = 0; j < n; j++) {
168  if (fabs(d_coords[i][j]) > max) {
169  max = fabs(d_coords[i][j]);
170  }
171  }
172  }
173  for (i = 0; i < dim; i++) {
174  for (j = 0; j < n; j++) {
175  d_coords[i][j] *= 10 / max;
176  }
177  }
178 
179  /**************************
180  ** Layout initialization **
181  **************************/
182 
183  for (i = 0; i < dim; i++) {
184  orthog1(n, d_coords[i]);
185  }
186 
187  /* for the y-coords, don't center them, but translate them so y[0]=0 */
188  y_0 = d_coords[1][0];
189  for (i = 0; i < n; i++) {
190  d_coords[1][i] -= y_0;
191  }
192  if (Verbose)
193  fprintf(stderr, ": %.2f sec", elapsed_sec());
194 
195  /**************************
196  ** Laplacian computation **
197  **************************/
198 
199  lap2 = Dij;
200  lap_length = n + n * (n - 1) / 2;
201  square_vec(lap_length, lap2);
202  /* compute off-diagonal entries */
203  invert_vec(lap_length, lap2);
204 
205  if (opt->clusters->nclusters > 0) {
206  int nn = n + opt->clusters->nclusters * 2;
207  int clap_length = nn + nn * (nn - 1) / 2;
208  float *clap = N_GNEW(clap_length, float);
209  int c0, c1;
210  float v;
211  c0 = c1 = 0;
212  for (i = 0; i < nn; i++) {
213  for (j = 0; j < nn - i; j++) {
214  if (i < n && j < n - i) {
215  v = lap2[c0++];
216  } else {
217  /* v=j==1?i%2:0; */
218  if (j == 1 && i % 2 == 1) {
219  v = maxEdgeLen;
220  v *= v;
221  if (v > 0.01) {
222  v = 1.0 / v;
223  }
224  } else
225  v = 0;
226  }
227  clap[c1++] = v;
228  }
229  }
230  free(lap2);
231  lap2 = clap;
232  n = nn;
233  lap_length = clap_length;
234  }
235  /* compute diagonal entries */
236  count = 0;
237  degrees = N_GNEW(n, double);
238  set_vector_val(n, 0, degrees);
239  for (i = 0; i < n - 1; i++) {
240  degree = 0;
241  count++; /* skip main diag entry */
242  for (j = 1; j < n - i; j++, count++) {
243  val = lap2[count];
244  degree += val;
245  degrees[i + j] -= val;
246  }
247  degrees[i] -= degree;
248  }
249  for (step = n, count = 0, i = 0; i < n; i++, count += step, step--) {
250  lap2[count] = (float) degrees[i];
251  }
252 
253  coords = N_GNEW(dim, float *);
254  f_storage = N_GNEW(dim * n, float);
255  for (i = 0; i < dim; i++) {
256  coords[i] = f_storage + i * n;
257  for (j = 0; j < n; j++) {
258  coords[i][j] = j < orig_n ? (float) (d_coords[i][j]) : 0;
259  }
260  }
261 
262  /* compute constant term in stress sum
263  * which is \sum_{i<j} w_{ij}d_{ij}^2
264  */
265  constant_term = (float) (n * (n - 1) / 2);
266 
267  /*************************
268  ** Layout optimization **
269  *************************/
270 
271  b = N_GNEW(dim, float *);
272  b[0] = N_GNEW(dim * n, float);
273  for (k = 1; k < dim; k++) {
274  b[k] = b[0] + k * n;
275  }
276 
277  tmp_coords = N_GNEW(n, float);
278  dist_accumulator = N_GNEW(n, float);
279 
280  old_stress = DBL_MAX; /* at least one iteration */
281 
282  if ((cMajEnvHor = initCMajVPSC(n, lap2, graph, opt, 0)) == NULL) {
283  iterations = -1;
284  goto finish;
285  }
286  if ((cMajEnvVrt = initCMajVPSC(n, lap2, graph, opt, opt->diredges)) == NULL) {
287  iterations = -1;
288  goto finish;
289  }
290 
291  lap1 = N_GNEW(lap_length, float);
292 
293  for (converged = FALSE, iterations = 0;
294  iterations < maxi && !converged; iterations++) {
295 
296  /* First, construct Laplacian of 1/(d_ij*|p_i-p_j|) */
297  set_vector_val(n, 0, degrees);
298  sqrt_vecf(lap_length, lap2, lap1);
299  for (count = 0, i = 0; i < n - 1; i++) {
300  len = n - i - 1;
301  /* init 'dist_accumulator' with zeros */
302  set_vector_valf(n, 0, dist_accumulator);
303 
304  /* put into 'dist_accumulator' all squared distances
305  * between 'i' and 'i'+1,...,'n'-1
306  */
307  for (k = 0; k < dim; k++) {
308  set_vector_valf(len, coords[k][i], tmp_coords);
309  vectors_mult_additionf(len, tmp_coords, -1,
310  coords[k] + i + 1);
311  square_vec(len, tmp_coords);
312  vectors_additionf(len, tmp_coords, dist_accumulator,
313  dist_accumulator);
314  }
315 
316  /* convert to 1/d_{ij} */
317  invert_sqrt_vec(len, dist_accumulator);
318  /* detect overflows */
319  for (j = 0; j < len; j++) {
320  if (dist_accumulator[j] >= FLT_MAX
321  || dist_accumulator[j] < 0) {
322  dist_accumulator[j] = 0;
323  }
324  }
325 
326  count++; /* save place for the main diagonal entry */
327  degree = 0;
328  for (j = 0; j < len; j++, count++) {
329  val = lap1[count] *= dist_accumulator[j];
330  degree += val;
331  degrees[i + j + 1] -= val;
332  }
333  degrees[i] -= degree;
334  }
335  for (step = n, count = 0, i = 0; i < n; i++, count += step, step--) {
336  lap1[count] = (float) degrees[i];
337  }
338 
339  /* Now compute b[] (L^(X(t))*X(t)) */
340  for (k = 0; k < dim; k++) {
341  /* b[k] := lap1*coords[k] */
342  right_mult_with_vector_ff(lap1, n, coords[k], b[k]);
343  }
344 
345  /* compute new stress
346  * remember that the Laplacians are negated, so we subtract
347  * instead of add and vice versa
348  */
349  new_stress = 0;
350  for (k = 0; k < dim; k++) {
351  new_stress += vectors_inner_productf(n, coords[k], b[k]);
352  }
353  new_stress *= 2;
354  new_stress += constant_term; /* only after mult by 2 */
355  for (k = 0; k < dim; k++) {
356  right_mult_with_vector_ff(lap2, n, coords[k], tmp_coords);
357  new_stress -= vectors_inner_productf(n, coords[k], tmp_coords);
358  }
359 
360 #ifdef ALTERNATIVE_STRESS_CALC
361  {
362  double mat_stress = new_stress;
363  double compute_stress(float **coords, float *lap, int dim,
364  int n);
365  new_stress = compute_stress(coords, lap2, dim, n);
366  if (fabs(mat_stress - new_stress) /
367  min(mat_stress, new_stress) > 0.001) {
368  fprintf(stderr,
369  "Diff stress vals: %lf %lf (iteration #%d)\n",
370  mat_stress, new_stress, iterations);
371  }
372  }
373 #endif
374  /* check for convergence */
375  if (Verbose && (iterations % 1 == 0)) {
376  fprintf(stderr, "%.3f ", new_stress);
377  if (iterations % 10 == 0)
378  fprintf(stderr, "\n");
379  }
380  converged = new_stress < old_stress
381  && fabs(new_stress - old_stress) / fabs(old_stress + 1e-10) <
382  Epsilon;
383  /*converged = converged || (iterations>1 && new_stress>old_stress); */
384  /* in first iteration we allowed stress increase, which
385  * might result ny imposing constraints
386  */
387  old_stress = new_stress;
388 
389  /* in determining non-overlap constraints we gradually scale up the
390  * size of nodes to avoid local minima
391  */
392  if ((iterations >= maxi - 1 || converged) && opt->noverlap == 1
393  && nsizeScale < 0.999) {
394  nsizeScale += 0.1;
395  if (Verbose)
396  fprintf(stderr, "nsizescale=%f,iterations=%d\n",
397  nsizeScale, iterations);
398  iterations = 0;
399  converged = FALSE;
400  }
401 
402 
403  /* now we find the optimizer of trace(X'LX)+X'B by solving 'dim'
404  * system of equations, thereby obtaining the new coordinates.
405  * If we use the constraints (given by the var's: 'ordering',
406  * 'levels' and 'num_levels'), we cannot optimize
407  * trace(X'LX)+X'B by simply solving equations, but we have
408  * to use a quadratic programming solver
409  * note: 'lap2' is a packed symmetric matrix, that is its
410  * upper-triangular part is arranged in a vector row-wise
411  * also note: 'lap2' is really the negated laplacian (the
412  * laplacian is -'lap2')
413  */
414 
415  if (opt->noverlap == 1 && nsizeScale > 0.001) {
416  generateNonoverlapConstraints(cMajEnvHor, nsizeScale, coords,
417  0,
418  nsizeScale < 0.5 ? FALSE : TRUE,
419  opt);
420  }
421  if (cMajEnvHor->m > 0) {
422 #ifdef MOSEK
423  if (opt->mosek) {
424  mosek_quad_solve_sep(cMajEnvHor->mosekEnv, n, b[0],
425  coords[0]);
426  } else
427 #endif /* MOSEK */
428  constrained_majorization_vpsc(cMajEnvHor, b[0], coords[0],
429  localConstrMajorIterations);
430  } else {
431  /* if there are no constraints then use conjugate gradient
432  * optimisation which should be considerably faster
433  */
434  if (conjugate_gradient_mkernel(lap2, coords[0], b[0], n,
435  tolerance_cg, n) < 0) {
436  iterations = -1;
437  goto finish;
438  }
439  }
440  if (opt->noverlap == 1 && nsizeScale > 0.001) {
441  generateNonoverlapConstraints(cMajEnvVrt, nsizeScale, coords,
442  1, FALSE, opt);
443  }
444  if (cMajEnvVrt->m > 0) {
445 #ifdef MOSEK
446  if (opt->mosek) {
447  mosek_quad_solve_sep(cMajEnvVrt->mosekEnv, n, b[1],
448  coords[1]);
449  } else
450 #endif /* MOSEK */
451  if (constrained_majorization_vpsc(cMajEnvVrt, b[1], coords[1],
452  localConstrMajorIterations) < 0) {
453  iterations = -1;
454  goto finish;
455  }
456  } else {
457  conjugate_gradient_mkernel(lap2, coords[1], b[1], n,
458  tolerance_cg, n);
459  }
460  }
461  if (Verbose) {
462  fprintf(stderr, "\nfinal e = %f %d iterations %.2f sec\n",
463  new_stress, iterations, elapsed_sec());
464  }
465  deleteCMajEnvVPSC(cMajEnvHor);
466  deleteCMajEnvVPSC(cMajEnvVrt);
467 
468  if (opt->noverlap == 2) {
469  /* fprintf(stderr, "Removing overlaps as post-process...\n"); */
470  removeoverlaps(orig_n, coords, opt);
471  }
472 
473 finish:
474  if (coords != NULL) {
475  for (i = 0; i < dim; i++) {
476  for (j = 0; j < orig_n; j++) {
477  d_coords[i][j] = coords[i][j];
478  }
479  }
480  free(coords[0]);
481  free(coords);
482  }
483 
484  if (b) {
485  free(b[0]);
486  free(b);
487  }
488  free(tmp_coords);
489  free(dist_accumulator);
490  free(degrees);
491  free(lap2);
492  free(lap1);
493 
494  return iterations;
495 }
496 #endif /* IPSEPCOLA */
void vectors_additionf(int n, float *vector1, float *vector2, float *result)
Definition: matrix_ops.c:629
#define MAX(a, b)
Definition: agerror.c:17
void invert_sqrt_vec(int n, float *vec)
Definition: matrix_ops.c:752
void orthog1(int n, double *vec)
Definition: matrix_ops.c:319
void start_timer(void)
Definition: timing.c:45
EXTERN double Epsilon
Definition: globals.h:77
void invert_vec(int n, float *vec)
Definition: matrix_ops.c:714
int initLayout(vtx_data *graph, int n, int dim, double **coords, node_t **nodes)
Definition: stress.c:159
void set_vector_val(int n, double val, double *result)
Definition: matrix_ops.c:677
#define MODEL_CIRCUIT
Definition: neato.h:21
#define tolerance_cg
Definition: stress.h:24
void set_vector_valf(int n, float val, float *result)
Definition: matrix_ops.c:685
int agerr(agerrlevel_t level, const char *fmt,...)
Definition: agerror.c:141
void square_vec(int n, float *vec)
Definition: matrix_ops.c:705
Definition: cgraph.h:388
void vectors_mult_additionf(int n, float *vector1, float alpha, float *vector2)
Definition: matrix_ops.c:639
double vectors_inner_productf(int n, float *vector1, float *vector2)
Definition: matrix_ops.c:665
int
Definition: grammar.c:1264
#define max(x, y)
Definition: stress.c:794
#define MODEL_SUBSET
Definition: neato.h:22
int nedges
Definition: sparsegraph.h:80
void sqrt_vecf(int n, float *source, float *target)
Definition: matrix_ops.c:737
float * circuitModel(vtx_data *graph, int nG)
Definition: stress.c:199
float * compute_apsp_packed(vtx_data *graph, int n)
Definition: stress.c:772
#define MODEL_MDS
Definition: neato.h:23
Agraph_t * graph(char *name)
Definition: gv.cpp:38
int DistType
Definition: sparsegraph.h:92
#define NULL
Definition: logic.h:39
double elapsed_sec(void)
Definition: timing.c:50
EXTERN unsigned char Verbose
Definition: globals.h:64
float * compute_apsp_artifical_weights_packed(vtx_data *graph, int n)
Definition: stress.c:796
float * mdsModel(vtx_data *graph, int nG)
Definition: stress.c:739
Definition: cgraph.h:388
void right_mult_with_vector_ff(float *packed_matrix, int n, float *vector, float *result)
Definition: matrix_ops.c:594
int conjugate_gradient_mkernel(float *A, float *x, float *b, int n, double tol, int max_iterations)
Definition: conjgrad.c:168
#define N_GNEW(n, t)
Definition: agxbuf.c:20
#define FALSE
Definition: cgraph.h:35
#define TRUE
Definition: cgraph.h:38