Graphviz 13.0.0~dev.20250607.1528
Loading...
Searching...
No Matches
post_process.c
Go to the documentation of this file.
1/*************************************************************************
2 * Copyright (c) 2011 AT&T Intellectual Property
3 * All rights reserved. This program and the accompanying materials
4 * are made available under the terms of the Eclipse Public License v1.0
5 * which accompanies this distribution, and is available at
6 * https://www.eclipse.org/legal/epl-v10.html
7 *
8 * Contributors: Details at https://graphviz.org
9 *************************************************************************/
10
11#include "config.h"
12
13#include <assert.h>
14#include <time.h>
15#include <string.h>
16#include <math.h>
17#include <stdbool.h>
18#include <stdlib.h>
19#include <common/types.h>
20#include <common/globals.h>
21
24#include <neatogen/overlap.h>
26#include <neatogen/call_tri.h>
27#include <sfdpgen/sfdp.h>
28#include <util/alloc.h>
29#include <util/exit.h>
30#include <util/gv_math.h>
31#include <util/unused.h>
32
33#define node_degree(i) (ia[(i)+1] - ia[(i)])
34
36 /* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]|
37 */
39 int *ia, *ja, i, j, k, l, nz;
40 double *d;
41 double len, di, sum, sumd;
42
43 assert(SparseMatrix_is_symmetric(A, false));
44
46 ia = D->ia;
47 ja = D->ja;
48 if (D->type != MATRIX_TYPE_REAL){
49 free(D->a);
50 D->type = MATRIX_TYPE_REAL;
51 D->a = gv_calloc(D->nz, sizeof(double));
52 }
53 d = D->a;
54
55 int *mask = gv_calloc(D->m, sizeof(int));
56 for (i = 0; i < D->m; i++) mask[i] = -1;
57
58 for (i = 0; i < D->m; i++){
59 di = node_degree(i);
60 mask[i] = i;
61 for (j = ia[i]; j < ia[i+1]; j++){
62 if (i == ja[j]) continue;
63 mask[ja[j]] = i;
64 }
65 for (j = ia[i]; j < ia[i+1]; j++){
66 k = ja[j];
67 if (i == k) continue;
68 len = di + node_degree(k);
69 for (l = ia[k]; l < ia[k+1]; l++){
70 if (mask[ja[l]] == i) len--;
71 }
72 d[j] = len;
73 assert(len > 0);
74 }
75
76 }
77
78 sum = 0; sumd = 0;
79 nz = 0;
80 for (i = 0; i < D->m; i++){
81 for (j = ia[i]; j < ia[i+1]; j++){
82 if (i == ja[j]) continue;
83 nz++;
84 sum += distance(x, dim, i, ja[j]);
85 sumd += d[j];
86 }
87 }
88 sum /= nz; sumd /= nz;
89 sum = sum/sumd;
90
91 for (i = 0; i < D->m; i++){
92 for (j = ia[i]; j < ia[i+1]; j++){
93 if (i == ja[j]) continue;
94 d[j] = sum*d[j];
95 }
96 }
97
98 free(mask);
99 return D;
100}
101
102
104 int ideal_dist_scheme){
105 /* use up to dist 2 neighbor. This is used in overcoming pherical effect with ideal distance of
106 2-neighbors equal graph distance etc.
107 */
108 int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
109 int nz;
110 double *d, *w, *lambda;
111 double diag_d, diag_w, dist, s = 0, stop = 0, sbot = 0;
113
114 assert(SparseMatrix_is_symmetric(A, false));
115
117
119 sm->scaling = 1.;
120 sm->data = NULL;
122 sm->tol_cg = 0.01;
123 sm->maxit_cg = floor(sqrt(A->m));
124
125 lambda = sm->lambda = gv_calloc(m, sizeof(double));
126 for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
127 int *mask = gv_calloc(m, sizeof(int));
128
129 double *avg_dist = gv_calloc(m, sizeof(double));
130
131 for (i = 0; i < m ;i++){
132 avg_dist[i] = 0;
133 nz = 0;
134 for (j = ia[i]; j < ia[i+1]; j++){
135 if (i == ja[j]) continue;
136 avg_dist[i] += distance(x, dim, i, ja[j]);
137 nz++;
138 }
139 assert(nz > 0);
140 avg_dist[i] /= nz;
141 }
142
143
144 for (i = 0; i < m; i++) mask[i] = -1;
145
146 nz = 0;
147 for (i = 0; i < m; i++){
148 mask[i] = i;
149 for (j = ia[i]; j < ia[i+1]; j++){
150 k = ja[j];
151 if (mask[k] != i){
152 mask[k] = i;
153 nz++;
154 }
155 }
156 for (j = ia[i]; j < ia[i+1]; j++){
157 k = ja[j];
158 for (l = ia[k]; l < ia[k+1]; l++){
159 if (mask[ja[l]] != i){
160 mask[ja[l]] = i;
161 nz++;
162 }
163 }
164 }
165 }
166
167 sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
168 assert(sm->Lw != NULL);
169 sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
170 assert(sm->Lwd != NULL);
171
172 iw = sm->Lw->ia; jw = sm->Lw->ja;
173
174 w = sm->Lw->a;
175 d = sm->Lwd->a;
176
177 id = sm->Lwd->ia; jd = sm->Lwd->ja;
178 iw[0] = id[0] = 0;
179
180 nz = 0;
181 for (i = 0; i < m; i++){
182 mask[i] = i+m;
183 diag_d = diag_w = 0;
184 for (j = ia[i]; j < ia[i+1]; j++){
185 k = ja[j];
186 if (mask[k] != i+m){
187 mask[k] = i+m;
188
189 jw[nz] = k;
190 if (ideal_dist_scheme == IDEAL_GRAPH_DIST){
191 dist = 1;
192 } else if (ideal_dist_scheme == IDEAL_AVG_DIST){
193 dist = (avg_dist[i] + avg_dist[k])*0.5;
194 } else if (ideal_dist_scheme == IDEAL_POWER_DIST){
195 dist = pow(distance_cropped(x,dim,i,k),.4);
196 } else {
197 fprintf(stderr,"ideal_dist_scheme value wrong");
198 assert(0);
199 graphviz_exit(1);
200 }
201
202 /*
203 w[nz] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
204 w[nz] = -2./(avg_dist[i]+avg_dist[k]);*/
205 /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
206 w[nz] = -1/(dist*dist);
207
208 diag_w += w[nz];
209
210 jd[nz] = k;
211 /*
212 d[nz] = w[nz]*distance(x,dim,i,k);
213 d[nz] = w[nz]*dd[j];
214 d[nz] = w[nz]*(avg_dist[i] + avg_dist[k])*0.5;
215 */
216 d[nz] = w[nz]*dist;
217 stop += d[nz]*distance(x,dim,i,k);
218 sbot += d[nz]*dist;
219 diag_d += d[nz];
220
221 nz++;
222 }
223 }
224
225 /* distance 2 neighbors */
226 for (j = ia[i]; j < ia[i+1]; j++){
227 k = ja[j];
228 for (l = ia[k]; l < ia[k+1]; l++){
229 if (mask[ja[l]] != i+m){
230 mask[ja[l]] = i+m;
231
232 if (ideal_dist_scheme == IDEAL_GRAPH_DIST){
233 dist = 2;
234 } else if (ideal_dist_scheme == IDEAL_AVG_DIST){
235 dist = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
236 } else if (ideal_dist_scheme == IDEAL_POWER_DIST){
237 dist = pow(distance_cropped(x,dim,i,ja[l]),.4);
238 } else {
239 fprintf(stderr,"ideal_dist_scheme value wrong");
240 assert(0);
241 graphviz_exit(1);
242 }
243
244 jw[nz] = ja[l];
245 /*
246 w[nz] = -1/(ia[i+1]-ia[i]+ia[ja[l]+1]-ia[ja[l]]);
247 w[nz] = -2/(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]]);*/
248 /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
249
250 w[nz] = -1/(dist*dist);
251
252 diag_w += w[nz];
253
254 jd[nz] = ja[l];
255 /*
256 d[nz] = w[nz]*(distance(x,dim,i,k)+distance(x,dim,k,ja[l]));
257 d[nz] = w[nz]*(dd[j]+dd[l]);
258 d[nz] = w[nz]*(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
259 */
260 d[nz] = w[nz]*dist;
261 stop += d[nz]*distance(x,dim,ja[l],k);
262 sbot += d[nz]*dist;
263 diag_d += d[nz];
264
265 nz++;
266 }
267 }
268 }
269 jw[nz] = i;
270 lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
271
272 w[nz] = -diag_w + lambda[i];
273 jd[nz] = i;
274 d[nz] = -diag_d;
275 nz++;
276
277 iw[i+1] = nz;
278 id[i+1] = nz;
279 }
280 s = stop/sbot;
281 for (i = 0; i < nz; i++) d[i] *= s;
282
283 sm->scaling = s;
284 sm->Lw->nz = nz;
285 sm->Lwd->nz = nz;
286
287 free(mask);
288 free(avg_dist);
290 return sm;
291}
292
294 /* solve a stress model to achieve the ideal distance among a sparse set of edges recorded in A.
295 A must be a real matrix.
296 */
297 int m = A->m;
298 double stop = 0, sbot = 0;
299
300 assert(SparseMatrix_is_symmetric(A, false) && A->type == MATRIX_TYPE_REAL);
301
302 /* if x is all zero, make it random */
303 bool has_nonzero = false;
304 for (int i = 0; i < m * dim; i++) {
305 if (!is_exactly_equal(x[i], 0)) {
306 has_nonzero = true;
307 break;
308 }
309 }
310 if (!has_nonzero) {
311 for (int i = 0; i < m*dim; i++) x[i] = 72*drand();
312 }
313
314 const int *const ia = A->ia;
315 const int *const ja = A->ja;
316 const double *const a = A->a;
317
319 sm->scaling = 1.;
320 sm->data = NULL;
322 sm->D = A;
323 sm->tol_cg = 0.01;
324 sm->maxit_cg = floor(sqrt(A->m));
325
326 double *const lambda = sm->lambda = gv_calloc(m, sizeof(double));
327
328 int nz = A->nz;
329
330 sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
331 assert(sm->Lw != NULL);
332 sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
333 assert(sm->Lwd != NULL);
334
335 int *const iw = sm->Lw->ia;
336 int *const jw = sm->Lw->ja;
337 int *const id = sm->Lwd->ia;
338 int *const jd = sm->Lwd->ja;
339 double *const w = sm->Lw->a;
340 double *const d = sm->Lwd->a;
341 iw[0] = id[0] = 0;
342
343 nz = 0;
344 for (int i = 0; i < m; i++){
345 double diag_d = 0;
346 double diag_w = 0;
347 for (int j = ia[i]; j < ia[i+1]; j++){
348 const int k = ja[j];
349 if (k != i){
350
351 jw[nz] = k;
352 const double dist = a[j];
353 w[nz] = -1;
354 diag_w += w[nz];
355 jd[nz] = k;
356 d[nz] = w[nz]*dist;
357
358 stop += d[nz]*distance(x,dim,i,k);
359 sbot += d[nz]*dist;
360 diag_d += d[nz];
361
362 nz++;
363 }
364 }
365
366 jw[nz] = i;
367 lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
368 w[nz] = -diag_w + lambda[i];
369
370 jd[nz] = i;
371 d[nz] = -diag_d;
372 nz++;
373
374 iw[i+1] = nz;
375 id[i+1] = nz;
376 }
377 const double s = stop / sbot;
378 if (s == 0) {
380 return NULL;
381 }
382 for (int i = 0; i < nz; i++) d[i] *= s;
383
384
385 sm->scaling = s;
386 sm->Lw->nz = nz;
387 sm->Lwd->nz = nz;
388
389 return sm;
390}
391
392static double total_distance(int m, int dim, double* x, double* y){
393 double total = 0, dist = 0;
394 int i, j;
395
396 for (i = 0; i < m; i++){
397 dist = 0.;
398 for (j = 0; j < dim; j++){
399 dist += (y[i*dim+j] - x[i*dim+j])*(y[i*dim+j] - x[i*dim+j]);
400 }
401 total += sqrt(dist);
402 }
403 return total;
404
405}
406
407
408
412
413
415 return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm);
416}
417
418static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, double *x, SparseMatrix *LL, double **rhs){
419 int edge_labeling_scheme = data->edge_labeling_scheme;
420 int n_constr_nodes = data->n_constr_nodes;
421 int *constr_nodes = data->constr_nodes;
422 SparseMatrix A_constr = data->A_constr;
423 int *ia = A_constr->ia, *ja = A_constr->ja, ii, jj, nz, l, ll, i, j;
424 int *irn = data->irn, *jcn = data->jcn;
425 double *val = data->val, dist, kk, k;
426 double *x00 = NULL;
427 SparseMatrix Lc = NULL;
428 double constr_penalty = data->constr_penalty;
429
430 if (edge_labeling_scheme == ELSCHEME_PENALTY || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY){
431 /* for an node with two neighbors j--i--k, and assume i needs to be between j and k, then the contribution to P is
432 . i j k
433 i 1 -0.5 -0.5
434 j -0.5 0.25 0.25
435 k -0.5 0.25 0.25
436 in general, if i has m neighbors j, k, ..., then
437 pᵢᵢ = 1
438 pⱼⱼ = 1 ÷ m²
439 pᵢⱼ = -1 ÷ m
440 pⱼₖ = 1 ÷ m²
441 . i j k
442 i 1 -1 ÷ m -1 ÷ m ...
443 j -1 ÷ m 1 ÷ m² 1 ÷ m² ...
444 k -1 ÷ m 1 ÷ m² 1 ÷ m² ...
445 */
446 if (!irn){
447 assert((!jcn) && (!val));
448 nz = 0;
449 for (i = 0; i < n_constr_nodes; i++){
450 ii = constr_nodes[i];
451 k = ia[ii+1] - ia[ii];/*usually k = 2 */
452 nz += (int)((k+1)*(k+1));
453
454 }
455 irn = data->irn = gv_calloc(nz, sizeof(int));
456 jcn = data->jcn = gv_calloc(nz, sizeof(int));
457 val = data->val = gv_calloc(nz, sizeof(double));
458 }
459 nz = 0;
460 for (i = 0; i < n_constr_nodes; i++){
461 ii = constr_nodes[i];
462 jj = ja[ia[ii]]; ll = ja[ia[ii] + 1];
463 if (jj == ll) continue; /* do not do loops */
464 dist = distance_cropped(x, dim, jj, ll);
465 dist *= dist;
466
467 k = ia[ii+1] - ia[ii];/* usually k = 2 */
468 kk = k*k;
469 irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(dist);
470 k = constr_penalty/(k*dist); kk = constr_penalty/(kk*dist);
471 for (j = ia[ii]; j < ia[ii+1]; j++){
472 irn[nz] = ii; jcn[nz] = ja[j]; val[nz++] = -k;
473 }
474 for (j = ia[ii]; j < ia[ii+1]; j++){
475 jj = ja[j];
476 irn[nz] = jj; jcn[nz] = ii; val[nz++] = -k;
477 for (l = ia[ii]; l < ia[ii+1]; l++){
478 ll = ja[l];
479 irn[nz] = jj; jcn[nz] = ll; val[nz++] = kk;
480 }
481 }
482 }
483 Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(double));
484 } else if (edge_labeling_scheme == ELSCHEME_PENALTY2 || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY2){
485 /* for an node with two neighbors j--i--k, and assume i needs to be between the old position of j and k, then the contribution to P is
486 1/d_jk, and to the right hand side: {0,...,average_position_of_i's neighbor if i is an edge node,...}
487 */
488 if (!irn){
489 assert((!jcn) && (!val));
490 nz = n_constr_nodes;
491 irn = data->irn = gv_calloc(nz, sizeof(int));
492 jcn = data->jcn = gv_calloc(nz, sizeof(int));
493 val = data->val = gv_calloc(nz, sizeof(double));
494 }
495 x00 = gv_calloc(m * dim, sizeof(double));
496 nz = 0;
497 for (i = 0; i < n_constr_nodes; i++){
498 ii = constr_nodes[i];
499 jj = ja[ia[ii]]; ll = ja[ia[ii] + 1];
500 dist = distance_cropped(x, dim, jj, ll);
501 irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(dist);
502 for (j = ia[ii]; j < ia[ii+1]; j++){
503 jj = ja[j];
504 for (l = 0; l < dim; l++){
505 x00[ii*dim+l] += x[jj*dim+l];
506 }
507 }
508 for (l = 0; l < dim; l++) {
509 x00[ii*dim+l] *= constr_penalty/(dist)/(ia[ii+1] - ia[ii]);
510 }
511 }
512 Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(double));
513
514 }
515 *LL = Lc;
516 *rhs = x00;
517}
518
519static UNUSED double get_stress(int m, int dim, int *iw, int *jw, double *w,
520 double *d, double *x, double scaling) {
521 int i, j;
522 double res = 0., dist;
523 // We use the fact that dᵢⱼ = wᵢⱼ × graph_dist(i, j). Also, dᵢⱼ and x are
524 // scaled by ×scaling, so divide by it to get actual unscaled stress.
525 for (i = 0; i < m; i++){
526 for (j = iw[i]; j < iw[i+1]; j++){
527 if (i == jw[j]) {
528 continue;
529 }
530 dist = d[j]/w[j];/* both negative*/
531 res += -w[j]*(dist - distance(x, dim, i, jw[j]))*(dist - distance(x, dim, i, jw[j]));
532 }
533 }
534 return 0.5*res/scaling/scaling;
535
536}
537
538double StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, double *x, int maxit_sm) {
539 SparseMatrix Lw = sm->Lw, Lwd = sm->Lwd, Lwdd = NULL;
540 int i, j, k, m, *id, *jd, *iw, *jw, idiag, iter = 0;
541 double *w, *dd, *d, *y = NULL, *x0 = NULL, *x00 = NULL, diag, diff = 1, *lambda = sm->lambda;
542 SparseMatrix Lc = NULL;
543 double dij, dist;
544
545 const double tol = 0.001;
546
547
548 Lwdd = SparseMatrix_copy(Lwd);
549 m = Lw->m;
550 x0 = calloc(dim*m, sizeof(double));
551 if (!x0) goto RETURN;
552
553 memcpy(x0, x, sizeof(double)*dim*m);
554 y = calloc(dim*m, sizeof(double));
555 if (!y) goto RETURN;
556
557 id = Lwd->ia; jd = Lwd->ja;
558 d = Lwd->a;
559 dd = Lwdd->a;
560 w = Lw->a;
561 iw = Lw->ia; jw = Lw->ja;
562
563#ifdef DEBUG_PRINT
564 if (Verbose) fprintf(stderr, "initial stress = %f\n", get_stress(m, dim, iw, jw, w, d, x, sm->scaling));
565#else
566 (void)iw;
567 (void)jw;
568#endif
569 /* for the additional matrix L due to the position constraints */
571 get_edge_label_matrix(sm->data, m, dim, x, &Lc, &x00);
572 if (Lc) Lw = SparseMatrix_add(Lw, Lc);
573 }
574
575 while (iter++ < maxit_sm && diff > tol){
576
577 for (i = 0; i < m; i++){
578 idiag = -1;
579 diag = 0.;
580 for (j = id[i]; j < id[i+1]; j++){
581 if (i == jd[j]) {
582 idiag = j;
583 continue;
584 }
585
586 dist = distance(x, dim, i, jd[j]);
587 if (d[j] == 0){
588 dd[j] = 0;
589 } else {
590 if (dist == 0){
591 dij = d[j]/w[j];/* the ideal distance */
592 /* perturb so points do not sit at the same place */
593 for (k = 0; k < dim; k++) x[jd[j]*dim+k] += 0.0001*(drand()+.0001)*dij;
594 dist = distance(x, dim, i, jd[j]);
595 }
596 dd[j] = d[j]/dist;
597 }
598 diag += dd[j];
599 }
600 assert(idiag >= 0);
601 dd[idiag] = -diag;
602 }
603 /* solve (Lw+lambda*I) x = Lwdd y + lambda x0 */
604
605 SparseMatrix_multiply_dense(Lwdd, x, y, dim);
606
607 if (lambda){/* is there a penalty term? */
608 for (i = 0; i < m; i++){
609 for (j = 0; j < dim; j++){
610 y[i*dim+j] += lambda[i]*x0[i*dim+j];
611 }
612 }
613 }
614
615 /* additional term added to the rhs */
616 switch (sm->scheme){
618 for (i = 0; i < m; i++){
619 for (j = 0; j < dim; j++){
620 y[i*dim+j] += x00[i*dim+j];
621 }
622 }
623 break;
624 }
625 default:
626 break;
627 }
628
629#ifdef DEBUG_PRINT
630 if (Verbose) {
631 fprintf(stderr, "stress1 = %g\n",get_stress(m, dim, iw, jw, w, d, x, sm->scaling));
632 }
633#endif
634
635 SparseMatrix_solve(Lw, dim, x, y, sm->tol_cg, sm->maxit_cg);
636
637#ifdef DEBUG_PRINT
638 if (Verbose) fprintf(stderr, "stress2 = %g\n",get_stress(m, dim, iw, jw, w, d, y, sm->scaling));
639#endif
640 diff = total_distance(m, dim, x, y)/sqrt(vector_product(m*dim, x, x));
641#ifdef DEBUG_PRINT
642 if (Verbose){
643 fprintf(stderr, "Outer iter = %d, cg res = %g, ||x_{k+1}-x_k||/||x_k|| = %g\n",iter, res, diff);
644 }
645#endif
646
647
648 memcpy(x, y, sizeof(double)*m*dim);
649 }
650
651#ifdef DEBUG
652 _statistics[1] += iter-1;
653#endif
654
655#ifdef DEBUG_PRINT
656 if (Verbose) fprintf(stderr, "iter = %d, final stress = %f\n", iter, get_stress(m, dim, iw, jw, w, d, x, sm->scaling));
657#endif
658
659 RETURN:
661 if (Lc) {
664 }
665
666 free(x0);
667 free(y);
668 free(x00);
669 return diff;
670
671}
672
674 if (!sm) return;
675 if (sm->Lw) SparseMatrix_delete(sm->Lw);
676 if (sm->Lwd) SparseMatrix_delete(sm->Lwd);
677 free(sm->lambda);
678 if (sm->data) sm->data_deallocator(sm->data);
679 free(sm);
680}
681
683 bool use_triangularization) {
684 int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, jdiag, nz;
686 double *d, *w, diag_d, diag_w, dist;
687 double s = 0, stop = 0, sbot = 0;
688
689 assert(SparseMatrix_is_symmetric(A, false));
690
691 double *avg_dist = gv_calloc(m, sizeof(double));
692
693 for (i = 0; i < m ;i++){
694 avg_dist[i] = 0;
695 nz = 0;
696 for (j = ia[i]; j < ia[i+1]; j++){
697 if (i == ja[j]) continue;
698 avg_dist[i] += distance(x, dim, i, ja[j]);
699 nz++;
700 }
701 assert(nz > 0);
702 avg_dist[i] /= nz;
703 }
704
706 sm->scaling = 1;
707 sm->data = NULL;
709 sm->tol_cg = 0.01;
710 sm->maxit_cg = floor(sqrt(A->m));
711
712 double *lambda = sm->lambda = gv_calloc(m, sizeof(double));
713
714 if (m > 2){
715 if (use_triangularization){
716 B= call_tri(m, x);
717 } else {
718 B= call_tri2(m, dim, x);
719 }
720 } else {
722 }
723
724
725
726 sm->Lw = SparseMatrix_add(A, B);
727
729 sm->Lwd = SparseMatrix_copy(sm->Lw);
730 if (!sm->Lw || !sm->Lwd) {
732 return NULL;
733 }
734
735 iw = sm->Lw->ia; jw = sm->Lw->ja;
736
737 w = sm->Lw->a;
738 d = sm->Lwd->a;
739
740 for (i = 0; i < m; i++){
741 diag_d = diag_w = 0;
742 jdiag = -1;
743 for (j = iw[i]; j < iw[i+1]; j++){
744 k = jw[j];
745 if (k == i){
746 jdiag = j;
747 continue;
748 }
749 /* w[j] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
750 w[j] = -2./(avg_dist[i]+avg_dist[k]);
751 w[j] = -1.*/;/* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
752 dist = pow(distance_cropped(x,dim,i,k),0.6);
753 w[j] = 1/(dist*dist);
754 diag_w += w[j];
755
756 /* d[j] = w[j]*distance(x,dim,i,k);
757 d[j] = w[j]*(avg_dist[i] + avg_dist[k])*0.5;*/
758 d[j] = w[j]*dist;
759 stop += d[j]*distance(x,dim,i,k);
760 sbot += d[j]*dist;
761 diag_d += d[j];
762
763 }
764
765 lambda[i] *= (-diag_w);
766
767 assert(jdiag >= 0);
768 w[jdiag] = -diag_w + lambda[i];
769 d[jdiag] = -diag_d;
770 }
771
772 s = stop/sbot;
773 for (i = 0; i < iw[m]; i++) d[i] *= s;
774 sm->scaling = s;
775
776 free(avg_dist);
777
778 return sm;
779}
780
786
790
791/* ================================ spring and spring-electrical based smoother ================ */
793 int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *id, *jd;
794 int nz;
795 double *d, *dd;
797
798 assert(SparseMatrix_is_symmetric(A, false));
799
801 dd = ID->a;
802
803 SpringSmoother sm = gv_alloc(sizeof(struct SpringSmoother_struct));
804 int *mask = gv_calloc(m, sizeof(int));
805
806 double *avg_dist = gv_calloc(m, sizeof(double));
807
808 for (i = 0; i < m ;i++){
809 avg_dist[i] = 0;
810 nz = 0;
811 for (j = ia[i]; j < ia[i+1]; j++){
812 if (i == ja[j]) continue;
813 avg_dist[i] += distance(x, dim, i, ja[j]);
814 nz++;
815 }
816 assert(nz > 0);
817 avg_dist[i] /= nz;
818 }
819
820
821 for (i = 0; i < m; i++) mask[i] = -1;
822
823 nz = 0;
824 for (i = 0; i < m; i++){
825 mask[i] = i;
826 for (j = ia[i]; j < ia[i+1]; j++){
827 k = ja[j];
828 if (mask[k] != i){
829 mask[k] = i;
830 nz++;
831 }
832 }
833 for (j = ia[i]; j < ia[i+1]; j++){
834 k = ja[j];
835 for (l = ia[k]; l < ia[k+1]; l++){
836 if (mask[ja[l]] != i){
837 mask[ja[l]] = i;
838 nz++;
839 }
840 }
841 }
842 }
843
845 assert(sm->D != NULL);
846
847 id = sm->D->ia; jd = sm->D->ja;
848 d = sm->D->a;
849 id[0] = 0;
850
851 nz = 0;
852 for (i = 0; i < m; i++){
853 mask[i] = i+m;
854 for (j = ia[i]; j < ia[i+1]; j++){
855 k = ja[j];
856 if (mask[k] != i+m){
857 mask[k] = i+m;
858 jd[nz] = k;
859 d[nz] = (avg_dist[i] + avg_dist[k])*0.5;
860 d[nz] = dd[j];
861 nz++;
862 }
863 }
864
865 for (j = ia[i]; j < ia[i+1]; j++){
866 k = ja[j];
867 for (l = ia[k]; l < ia[k+1]; l++){
868 if (mask[ja[l]] != i+m){
869 mask[ja[l]] = i+m;
870 jd[nz] = ja[l];
871 d[nz] = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
872 d[nz] = dd[j]+dd[l];
873 nz++;
874 }
875 }
876 }
877 id[i+1] = nz;
878 }
879 sm->D->nz = nz;
880 sm->ctrl = ctrl;
881 sm->ctrl.random_start = false;
882 sm->ctrl.multilevels = 1;
883 sm->ctrl.step /= 2;
884 sm->ctrl.maxiter = 20;
885
886 free(mask);
887 free(avg_dist);
889
890 return sm;
891}
892
893
895 if (!sm) return;
896 if (sm->D) SparseMatrix_delete(sm->D);
897}
898
900 int flag = 0;
901
902 spring_electrical_spring_embedding(dim, A, sm->D, &sm->ctrl, x, &flag);
903 assert(!flag);
904
905}
906
907/*=============================== end of spring and spring-electrical based smoother =========== */
908
910#ifdef TIME
911 clock_t cpu;
912#endif
913
914#ifdef TIME
915 cpu = clock();
916#endif
917
918 switch (ctrl.smoothing){
919 case SMOOTHING_RNG:
920 case SMOOTHING_TRIANGLE:{
922
923 if (A->m > 2) { /* triangles need at least 3 nodes */
924 if (ctrl.smoothing == SMOOTHING_RNG){
925 sm = TriangleSmoother_new(A, dim, x, false);
926 } else {
927 sm = TriangleSmoother_new(A, dim, x, true);
928 }
931 }
932 break;
933 }
937 {
939 int dist_scheme = IDEAL_AVG_DIST;
940
942 dist_scheme = IDEAL_GRAPH_DIST;
944 dist_scheme = IDEAL_AVG_DIST;
946 dist_scheme = IDEAL_POWER_DIST;
947 }
948
949 sm = StressMajorizationSmoother2_new(A, dim, 0.05, x, dist_scheme);
952 break;
953 }
954 case SMOOTHING_SPRING:{
955 SpringSmoother sm = SpringSmoother_new(A, dim, ctrl, x);
956 SpringSmoother_smooth(sm, A, dim, x);
958
959 break;
960 }
961
962 }/* end switch between smoothing methods */
963
964#ifdef TIME
965 if (Verbose) fprintf(stderr, "post processing %f\n",((double) (clock() - cpu)) / CLOCKS_PER_SEC);
966#endif
967}
void SparseMatrix_multiply_dense(SparseMatrix A, const double *v, double *res, int dim)
bool SparseMatrix_is_symmetric(SparseMatrix A, bool test_pattern_symmetry_only)
SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz)
void SparseMatrix_delete(SparseMatrix A)
SparseMatrix SparseMatrix_copy(SparseMatrix A)
SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B)
SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format)
@ FORMAT_CSR
@ MATRIX_TYPE_REAL
Memory allocation wrappers that exit on failure.
static void * gv_calloc(size_t nmemb, size_t size)
Definition alloc.h:26
static void * gv_alloc(size_t size)
Definition alloc.h:47
SparseMatrix call_tri2(int n, int dim, double *xx)
Definition call_tri.c:64
SparseMatrix call_tri(int n, double *x)
Definition call_tri.c:19
static NORETURN void graphviz_exit(int status)
Definition exit.h:23
#define A(n, t)
Definition expr.h:76
static double dist(int dim, double *x, double *y)
double drand(void)
Definition general.c:22
double distance(double *x, int dim, int i, int j)
Definition general.c:121
double vector_product(int n, double *x, double *y)
Definition general.c:32
double distance_cropped(double *x, int dim, int i, int j)
Definition general.c:116
static double len(glCompPoint p)
Definition glutils.c:150
static bool Verbose
Definition gml2gv.c:23
void free(void *)
#define ID
Definition gmlparse.h:134
node NULL
Definition grammar.y:180
static uint64_t id
Definition gv2gml.c:40
Arithmetic helper functions.
static bool is_exactly_equal(double a, double b)
are two values precisely the same?
Definition gv_math.h:45
#define B
Definition hierarchy.c:118
#define D
Definition hierarchy.c:120
static const int dim
@ ELSCHEME_PENALTY2
Definition overlap.h:29
@ ELSCHEME_STRAIGHTLINE_PENALTY2
Definition overlap.h:29
@ ELSCHEME_STRAIGHTLINE_PENALTY
Definition overlap.h:29
@ ELSCHEME_PENALTY
Definition overlap.h:29
static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, double *x, SparseMatrix *LL, double **rhs)
StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, double *x)
static double total_distance(int m, int dim, double *x, double *y)
static SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, double *x)
void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm)
TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, double *x, bool use_triangularization)
double StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, double *x, int maxit_sm)
void SpringSmoother_delete(SpringSmoother sm)
double SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, double *x, int maxit_sm)
void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, double *x)
void TriangleSmoother_smooth(TriangleSmoother sm, int dim, double *x)
void StressMajorizationSmoother_delete(StressMajorizationSmoother sm)
SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, double *x)
StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, double lambda0, double *x, int ideal_dist_scheme)
static UNUSED double get_stress(int m, int dim, int *iw, int *jw, double *w, double *d, double *x, double scaling)
void TriangleSmoother_delete(TriangleSmoother sm)
void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, int dim, double *x)
#define node_degree(i)
@ IDEAL_GRAPH_DIST
@ IDEAL_AVG_DIST
@ IDEAL_POWER_DIST
@ SM_SCHEME_NORMAL_ELABEL
@ SM_SCHEME_NORMAL
#define TriangleSmoother_struct
double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, double maxit)
static const double tol
void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control *ctrl, double *x, int *flag)
@ SMOOTHING_SPRING
@ SMOOTHING_RNG
@ SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST
@ SMOOTHING_STRESS_MAJORIZATION_POWER_DIST
@ SMOOTHING_TRIANGLE
@ SMOOTHING_STRESS_MAJORIZATION_AVG_DIST
#define RETURN(v)
Definition strmatch.c:144
spring_electrical_control ctrl
SparseMatrix Lwd
the laplacian like matrix with offdiag = -scaling × dᵢⱼ ÷ wᵢⱼ. RHS in stress majorization = Lwd....
SparseMatrix Lw
the weighted laplacian. with offdiag = -1 ÷ wᵢⱼ
Definition legal.c:50
bool random_start
whether to apply SE from a random layout, or from exisiting layout
graphs, nodes and edges info: Agraphinfo_t, Agnodeinfo_t and Agedgeinfo_t
Definition grammar.c:89
abstraction for squashing compiler warnings for unused symbols
#define UNUSED
Definition unused.h:25