1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14#include "config.h"
15
16#include <time.h>
17#include <string.h>
18#include <math.h>
19
20#include "types.h"
21#include "memory.h"
22#include "globals.h"
23
24#include "sparse_solve.h"
25#include "post_process.h"
26#include "overlap.h"
27#include "spring_electrical.h"
28#include "call_tri.h"
29#include "sfdpinternal.h"
30
31#define node_degree(i) (ia[(i)+1] - ia[(i)])
32
33#ifdef UNUSED
34static void get_neighborhood_precision_recall(char *outfile, SparseMatrix A0, real *ideal_dist_matrix, real *dist_matrix){
35 SparseMatrix A = A0;
36 int i, j, k, n = A->m;
37 // int *ia, *ja;
38 int *g_order = NULL, *p_order = NULL;/* ordering using graph/physical distance */
39 real *gdist, *pdist, radius;
40 int np_neighbors;
41 int ng_neighbors; /*number of (graph theoretical) neighbors */
42 real node_dist;/* distance of a node to the center node */
43 real true_positive;
44 real recall;
45 FILE *fp;
46
47 fp = fopen(outfile,"w");
48
49 if (!SparseMatrix_is_symmetric(A, FALSE)){
50 A = SparseMatrix_symmetrize(A, FALSE);
51 }
52 // ia = A->ia;
53 // ja = A->ja;
54
55 for (k = 5; k <= 50; k+= 5){
56 recall = 0;
57 for (i = 0; i < n; i++){
58 gdist = &(ideal_dist_matrix[i*n]);
59 vector_ordering(n, gdist, &g_order, TRUE);
60 pdist = &(dist_matrix[i*n]);
61 vector_ordering(n, pdist, &p_order, TRUE);
62 ng_neighbors = MIN(n-1, k); /* set the number of closest neighbor in the graph space to consider, excluding the node itself */
63 np_neighbors = ng_neighbors;/* set the number of closest neighbor in the embedding to consider, excluding the node itself */
64 radius = pdist[p_order[np_neighbors]];
65 true_positive = 0;
66 for (j = 1; j <= ng_neighbors; j++){
67 node_dist = pdist[g_order[j]];/* the phisical distance for j-th closest node (in graph space) */
68 if (node_dist <= radius) true_positive++;
69 }
70 recall += true_positive/np_neighbors;
71 }
72 recall /= n;
73
74 fprintf(fp,"%d %f\n", k, recall);
75 }
76 fprintf(stderr,"wrote precision/recall in file %s\n", outfile);
77 fclose(fp);
78
79 if (A != A0) SparseMatrix_delete(A);
80 FREE(g_order); FREE(p_order);
81}
82
83
84
85void dump_distance_edge_length(char *outfile, SparseMatrix A, int dim, real *x){
86 int weighted = TRUE;
87 int n, i, j, nzz;
88 real *dist_matrix = NULL;
89 int flag;
90 real *dij, *xij, wij, top = 0, bot = 0, t;
91
92 int *p = NULL;
93 real dmin, dmax, xmax, xmin, bsta, bwidth, *xbin, x25, x75, median;
94 int nbins = 30, nsta, nz = 0;
95 FILE *fp;
96
97 fp = fopen(outfile,"w");
98
99 flag = SparseMatrix_distance_matrix(A, weighted, &dist_matrix);
100 assert(!flag);
101
102 n = A->m;
103 dij = MALLOC(sizeof(real)*(n*(n-1)/2));
104 xij = MALLOC(sizeof(real)*(n*(n-1)/2));
105 for (i = 0; i < n; i++){
106 for (j = i+1; j < n; j++){
107 dij[nz] = dist_matrix[i*n+j];
108 xij[nz] = distance(x, dim, i, j);
109 if (dij[nz] > 0){
110 wij = 1/(dij[nz]*dij[nz]);
111 } else {
112 wij = 1;
113 }
114 top += wij * dij[nz] * xij[nz];
115 bot += wij*xij[nz]*xij[nz];
116 nz++;
117 }
118 }
119 if (bot > 0){
120 t = top/bot;
121 } else {
122 t = 1;
123 }
124 fprintf(stderr,"scaling factor = %f\n", t);
125
126 for (i = 0; i < nz; i++) xij[i] *= t;
127
128 vector_ordering(nz, dij, &p, TRUE);
129 dmin = dij[p[0]];
130 dmax = dij[p[nz-1]];
131 nbins = MIN(nbins, dmax/MAX(1,dmin));/* scale by dmin since edge length of 1 is treated as 72 points in stress/maxent/full_stress */
132 bwidth = (dmax - dmin)/nbins;
133
134 nsta = 0; bsta = dmin;
135 xbin = MALLOC(sizeof(real)*nz);
136 nzz = nz;
137
138 for (i = 0; i < nbins; i++){
139 /* the bin is [dmin + i*(dmax-dmin)/nbins, dmin + (i+1)*(dmax-dmin)/nbins] */
140 nz = 0; xmin = xmax = xij[p[nsta]];
141 while (nsta < nzz && dij[p[nsta]] >= bsta && dij[p[nsta]] <= bsta + bwidth){
142 xbin[nz++] = xij[p[nsta]];
143 xmin = MIN(xmin, xij[p[nsta]]);
144 xmax = MAX(xmax, xij[p[nsta]]);
145 nsta++;
146 }
147 /* find the median, and 25/75% */
148 if (nz > 0){
149 median = vector_median(nz, xbin);
150 x25 = vector_percentile(nz, xbin, 0.25);
151 x75 = vector_percentile(nz, xbin, 0.75);
152 fprintf(fp, "%d %f %f %f %f %f %f %f\n", nz, bsta, bsta + bwidth, xmin, x25, median, x75, xmax);
153 } else {
154 xmin = xmax = median = x25 = x75 = (bsta + 0.5*bwidth);
155 }
156 bsta += bwidth;
157 }
158
159 FREE(dij); FREE(xij); FREE(xbin); FREE(p);
160 FREE(dist_matrix);
161
162}
163
164real get_full_stress(SparseMatrix A, int dim, real *x, int weighting_scheme){
165 /* get the full weighted stress, \sum_ij w_ij (t ||x_i-x_j||^2-d_ij)^2, where t
166 is the optimal scaling factor t = \sum w_ij ||x_i-x_j||^2/(\sum w_ij d_ij ||x_i - x_j||),
167
168 Matrix A is assume to be sparse and a full distance matrix will be generated.
169 We assume undirected graph and only check an undirected edge i--j, not i->j and j->i.
170 */
171 int weighted = TRUE;
172 int n, i, j;
173 real *ideal_dist_matrix = NULL, *dist_matrix;
174 int flag;
175 real t, top = 0, bot = 0, lstress, stress = 0, dij, wij, xij;
176 real sne_disimilarity = 0;
177
178 flag = SparseMatrix_distance_matrix(A, weighted, &ideal_dist_matrix);
179 assert(!flag);
180
181 n = A->m;
182 dist_matrix = MALLOC(sizeof(real)*n*n);
183
184 for (i = 0; i < n; i++){
185 for (j = i+1; j < n; j++){
186 dij = ideal_dist_matrix[i*n+j];
187 assert(dij >= 0);
188 xij = distance(x, dim, i, j);
189 if (dij > 0){
190 switch (weighting_scheme){
191 case WEIGHTING_SCHEME_SQR_DIST:
192 wij = 1/(dij*dij);
193 break;
194 case WEIGHTING_SCHEME_INV_DIST:
195 wij = 1/(dij);
196 break;
197 break;
198 case WEIGHTING_SCHEME_NONE:
199 wij = 1;
200 default:
201 assert(0);
202 }
203 } else {
204 wij = 1;
205 }
206 top += wij * dij * xij;
207 bot += wij*xij*xij;
208 }
209 }
210 if (bot > 0){
211 t = top/bot;
212 } else {
213 t = 1;
214 }
215
216 for (i = 0; i < n; i++){
217 dist_matrix[i*n+i] = 0.;
218 for (j = i+1; j < n; j++){
219 dij = ideal_dist_matrix[i*n+j];
220 assert(dij >= 0);
221 xij = distance(x, dim, i, j);
222 dist_matrix[i*n+j] = xij*t;
223 dist_matrix[j*n+i] = xij*t;
224 if (dij > 0){
225 wij = 1/(dij*dij);
226 } else {
227 wij = 1;
228 }
229 lstress = (xij*t - dij);
230 stress += wij*lstress*lstress;
231 }
232 }
233
234 {int K = 20;
235 sne_disimilarity = get_sne_disimilarity(n, ideal_dist_matrix, dist_matrix, K);
236 }
237
238 get_neighborhood_precision_recall("/tmp/recall.txt", A, ideal_dist_matrix, dist_matrix);
239 get_neighborhood_precision_recall("/tmp/precision.txt", A, dist_matrix, ideal_dist_matrix);
240
241 fprintf(stderr,"sne_disimilarity = %f\n",sne_disimilarity);
242 if (n > 0) fprintf(stderr,"sstress per edge = %f\n",stress/n/(n-1)*2);
243
244 FREE(dist_matrix);
245 FREE(ideal_dist_matrix);
246 return stress;
247}
248#endif
249
250
251SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real *x){
252 /* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]|
253 */
254 SparseMatrix D;
255 int *ia, *ja, i, j, k, l, nz;
256 real *d;
257 int *mask = NULL;
258 real len, di, sum, sumd;
259
260 assert(SparseMatrix_is_symmetric(A, FALSE));
261
262 D = SparseMatrix_copy(A);
263 ia = D->ia;
264 ja = D->ja;
265 if (D->type != MATRIX_TYPE_REAL){
266 FREE(D->a);
267 D->type = MATRIX_TYPE_REAL;
268 D->a = N_GNEW(D->nz,real);
269 }
270 d = (real*) D->a;
271
272 mask = N_GNEW(D->m,int);
273 for (i = 0; i < D->m; i++) mask[i] = -1;
274
275 for (i = 0; i < D->m; i++){
276 di = node_degree(i);
277 mask[i] = i;
278 for (j = ia[i]; j < ia[i+1]; j++){
279 if (i == ja[j]) continue;
280 mask[ja[j]] = i;
281 }
282 for (j = ia[i]; j < ia[i+1]; j++){
283 k = ja[j];
284 if (i == k) continue;
285 len = di + node_degree(k);
286 for (l = ia[k]; l < ia[k+1]; l++){
287 if (mask[ja[l]] == i) len--;
288 }
289 d[j] = len;
290 assert(len > 0);
291 }
292
293 }
294
295 sum = 0; sumd = 0;
296 nz = 0;
297 for (i = 0; i < D->m; i++){
298 for (j = ia[i]; j < ia[i+1]; j++){
299 if (i == ja[j]) continue;
300 nz++;
301 sum += distance(x, dim, i, ja[j]);
302 sumd += d[j];
303 }
304 }
305 sum /= nz; sumd /= nz;
306 sum = sum/sumd;
307
308 for (i = 0; i < D->m; i++){
309 for (j = ia[i]; j < ia[i+1]; j++){
310 if (i == ja[j]) continue;
311 d[j] = sum*d[j];
312 }
313 }
314
315
316 return D;
317}
318
319
320StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, real lambda0, real *x,
321 int ideal_dist_scheme){
322 /* use up to dist 2 neighbor */
323 /* use up to dist 2 neighbor. This is used in overcoming pherical effect with ideal distance of
324 2-neighbors equal graph distance etc.
325 */
326 StressMajorizationSmoother sm;
327 int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
328 int *mask, nz;
329 real *d, *w, *lambda;
330 real *avg_dist, diag_d, diag_w, dist, s = 0, stop = 0, sbot = 0;
331 SparseMatrix ID;
332
333 assert(SparseMatrix_is_symmetric(A, FALSE));
334
335 ID = ideal_distance_matrix(A, dim, x);
336
337 sm = GNEW(struct StressMajorizationSmoother_struct);
338 sm->scaling = 1.;
339 sm->data = NULL;
340 sm->scheme = SM_SCHEME_NORMAL;
341 sm->tol_cg = 0.01;
342 sm->maxit_cg = (int)sqrt((double) A->m);
343
344 lambda = sm->lambda = N_GNEW(m,real);
345 for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
346 mask = N_GNEW(m,int);
347
348 avg_dist = N_GNEW(m,real);
349
350 for (i = 0; i < m ;i++){
351 avg_dist[i] = 0;
352 nz = 0;
353 for (j = ia[i]; j < ia[i+1]; j++){
354 if (i == ja[j]) continue;
355 avg_dist[i] += distance(x, dim, i, ja[j]);
356 nz++;
357 }
358 assert(nz > 0);
359 avg_dist[i] /= nz;
360 }
361
362
363 for (i = 0; i < m; i++) mask[i] = -1;
364
365 nz = 0;
366 for (i = 0; i < m; i++){
367 mask[i] = i;
368 for (j = ia[i]; j < ia[i+1]; j++){
369 k = ja[j];
370 if (mask[k] != i){
371 mask[k] = i;
372 nz++;
373 }
374 }
375 for (j = ia[i]; j < ia[i+1]; j++){
376 k = ja[j];
377 for (l = ia[k]; l < ia[k+1]; l++){
378 if (mask[ja[l]] != i){
379 mask[ja[l]] = i;
380 nz++;
381 }
382 }
383 }
384 }
385
386 sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
387 sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
388 if (!(sm->Lw) || !(sm->Lwd)) {
389 StressMajorizationSmoother_delete(sm);
390 return NULL;
391 }
392
393 iw = sm->Lw->ia; jw = sm->Lw->ja;
394
395 w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
396
397 id = sm->Lwd->ia; jd = sm->Lwd->ja;
398 iw[0] = id[0] = 0;
399
400 nz = 0;
401 for (i = 0; i < m; i++){
402 mask[i] = i+m;
403 diag_d = diag_w = 0;
404 for (j = ia[i]; j < ia[i+1]; j++){
405 k = ja[j];
406 if (mask[k] != i+m){
407 mask[k] = i+m;
408
409 jw[nz] = k;
410 if (ideal_dist_scheme == IDEAL_GRAPH_DIST){
411 dist = 1;
412 } else if (ideal_dist_scheme == IDEAL_AVG_DIST){
413 dist = (avg_dist[i] + avg_dist[k])*0.5;
414 } else if (ideal_dist_scheme == IDEAL_POWER_DIST){
415 dist = pow(distance_cropped(x,dim,i,k),.4);
416 } else {
417 fprintf(stderr,"ideal_dist_scheme value wrong");
418 assert(0);
419 exit(1);
420 }
421
422 /*
423 w[nz] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
424 w[nz] = -2./(avg_dist[i]+avg_dist[k]);*/
425 /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
426 w[nz] = -1/(dist*dist);
427
428 diag_w += w[nz];
429
430 jd[nz] = k;
431 /*
432 d[nz] = w[nz]*distance(x,dim,i,k);
433 d[nz] = w[nz]*dd[j];
434 d[nz] = w[nz]*(avg_dist[i] + avg_dist[k])*0.5;
435 */
436 d[nz] = w[nz]*dist;
437 stop += d[nz]*distance(x,dim,i,k);
438 sbot += d[nz]*dist;
439 diag_d += d[nz];
440
441 nz++;
442 }
443 }
444
445 /* distance 2 neighbors */
446 for (j = ia[i]; j < ia[i+1]; j++){
447 k = ja[j];
448 for (l = ia[k]; l < ia[k+1]; l++){
449 if (mask[ja[l]] != i+m){
450 mask[ja[l]] = i+m;
451
452 if (ideal_dist_scheme == IDEAL_GRAPH_DIST){
453 dist = 2;
454 } else if (ideal_dist_scheme == IDEAL_AVG_DIST){
455 dist = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
456 } else if (ideal_dist_scheme == IDEAL_POWER_DIST){
457 dist = pow(distance_cropped(x,dim,i,ja[l]),.4);
458 } else {
459 fprintf(stderr,"ideal_dist_scheme value wrong");
460 assert(0);
461 exit(1);
462 }
463
464 jw[nz] = ja[l];
465 /*
466 w[nz] = -1/(ia[i+1]-ia[i]+ia[ja[l]+1]-ia[ja[l]]);
467 w[nz] = -2/(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]]);*/
468 /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
469
470 w[nz] = -1/(dist*dist);
471
472 diag_w += w[nz];
473
474 jd[nz] = ja[l];
475 /*
476 d[nz] = w[nz]*(distance(x,dim,i,k)+distance(x,dim,k,ja[l]));
477 d[nz] = w[nz]*(dd[j]+dd[l]);
478 d[nz] = w[nz]*(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
479 */
480 d[nz] = w[nz]*dist;
481 stop += d[nz]*distance(x,dim,ja[l],k);
482 sbot += d[nz]*dist;
483 diag_d += d[nz];
484
485 nz++;
486 }
487 }
488 }
489 jw[nz] = i;
490 lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
491
492 w[nz] = -diag_w + lambda[i];
493 jd[nz] = i;
494 d[nz] = -diag_d;
495 nz++;
496
497 iw[i+1] = nz;
498 id[i+1] = nz;
499 }
500 s = stop/sbot;
501 for (i = 0; i < nz; i++) d[i] *= s;
502
503 sm->scaling = s;
504 sm->Lw->nz = nz;
505 sm->Lwd->nz = nz;
506
507 FREE(mask);
508 FREE(avg_dist);
509 SparseMatrix_delete(ID);
510 return sm;
511}
512
513StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x,
514 int weighting_scheme, int scale_initial_coord){
515 /* solve a stress model to achieve the ideal distance among a sparse set of edges recorded in A.
516 A must be a real matrix.
517 */
518 StressMajorizationSmoother sm;
519 int i, j, k, m = A->m, *ia, *ja, *iw, *jw, *id, *jd;
520 int nz;
521 real *d, *w, *lambda;
522 real diag_d, diag_w, *a, dist, s = 0, stop = 0, sbot = 0;
523 real xdot = 0;
524
525 assert(SparseMatrix_is_symmetric(A, FALSE) && A->type == MATRIX_TYPE_REAL);
526
527 /* if x is all zero, make it random */
528 for (i = 0; i < m*dim; i++) xdot += x[i]*x[i];
529 if (xdot == 0){
530 for (i = 0; i < m*dim; i++) x[i] = 72*drand();
531 }
532
533 ia = A->ia;
534 ja = A->ja;
535 a = (real*) A->a;
536
537
538 sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct));
539 sm->scaling = 1.;
540 sm->data = NULL;
541 sm->scheme = SM_SCHEME_NORMAL;
542 sm->D = A;
543 sm->tol_cg = 0.01;
544 sm->maxit_cg = (int)sqrt((double) A->m);
545
546 lambda = sm->lambda = MALLOC(sizeof(real)*m);
547 for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
548
549 nz = A->nz;
550
551 sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
552 sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
553 if (!(sm->Lw) || !(sm->Lwd)) {
554 StressMajorizationSmoother_delete(sm);
555 return NULL;
556 }
557
558 iw = sm->Lw->ia; jw = sm->Lw->ja;
559 id = sm->Lwd->ia; jd = sm->Lwd->ja;
560 w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
561 iw[0] = id[0] = 0;
562
563 nz = 0;
564 for (i = 0; i < m; i++){
565 diag_d = diag_w = 0;
566 for (j = ia[i]; j < ia[i+1]; j++){
567 k = ja[j];
568 if (k != i){
569
570 jw[nz] = k;
571 dist = a[j];
572 switch (weighting_scheme){
573 case WEIGHTING_SCHEME_SQR_DIST:
574 if (dist*dist == 0){
575 w[nz] = -100000;
576 } else {
577 w[nz] = -1/(dist*dist);
578 }
579 break;
580 case WEIGHTING_SCHEME_INV_DIST:
581 if (dist*dist == 0){
582 w[nz] = -100000;
583 } else {
584 w[nz] = -1/(dist);
585 }
586 break;
587 case WEIGHTING_SCHEME_NONE:
588 w[nz] = -1;
589 break;
590 default:
591 assert(0);
592 return NULL;
593 }
594 diag_w += w[nz];
595 jd[nz] = k;
596 d[nz] = w[nz]*dist;
597
598 stop += d[nz]*distance(x,dim,i,k);
599 sbot += d[nz]*dist;
600 diag_d += d[nz];
601
602 nz++;
603 }
604 }
605
606 jw[nz] = i;
607 lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
608 w[nz] = -diag_w + lambda[i];
609
610 jd[nz] = i;
611 d[nz] = -diag_d;
612 nz++;
613
614 iw[i+1] = nz;
615 id[i+1] = nz;
616 }
617 if (scale_initial_coord){
618 s = stop/sbot;
619 } else {
620 s = 1.;
621 }
622 if (s == 0) {
623 return NULL;
624 }
625 for (i = 0; i < nz; i++) d[i] *= s;
626
627
628 sm->scaling = s;
629 sm->Lw->nz = nz;
630 sm->Lwd->nz = nz;
631
632 return sm;
633}
634
635static real total_distance(int m, int dim, real* x, real* y){
636 real total = 0, dist = 0;
637 int i, j;
638
639 for (i = 0; i < m; i++){
640 dist = 0.;
641 for (j = 0; j < dim; j++){
642 dist += (y[i*dim+j] - x[i*dim+j])*(y[i*dim+j] - x[i*dim+j]);
643 }
644 total += sqrt(dist);
645 }
646 return total;
647
648}
649
650
651
652void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm){
653 StressMajorizationSmoother_delete(sm);
654}
655
656
657real SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol){
658
659 return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, tol);
660
661
662}
663static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, real *x, SparseMatrix *LL, real **rhs){
664 int edge_labeling_scheme = data->edge_labeling_scheme;
665 int n_constr_nodes = data->n_constr_nodes;
666 int *constr_nodes = data->constr_nodes;
667 SparseMatrix A_constr = data->A_constr;
668 int *ia = A_constr->ia, *ja = A_constr->ja, ii, jj, nz, l, ll, i, j;
669 int *irn = data->irn, *jcn = data->jcn;
670 real *val = data->val, dist, kk, k;
671 real *x00 = NULL;
672 SparseMatrix Lc = NULL;
673 real constr_penalty = data->constr_penalty;
674
675 if (edge_labeling_scheme == ELSCHEME_PENALTY || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY){
676 /* 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
677 . i j k
678 i 1 -0.5 -0.5
679 j -0.5 0.25 0.25
680 k -0.5 0.25 0.25
681 in general, if i has m neighbors j, k, ..., then
682 p_ii = 1
683 p_jj = 1/m^2
684 p_ij = -1/m
685 p jk = 1/m^2
686 . i j k
687 i 1 -1/m -1/m ...
688 j -1/m 1/m^2 1/m^2 ...
689 k -1/m 1/m^2 1/m^2 ...
690 */
691 if (!irn){
692 assert((!jcn) && (!val));
693 nz = 0;
694 for (i = 0; i < n_constr_nodes; i++){
695 ii = constr_nodes[i];
696 k = ia[ii+1] - ia[ii];/*usually k = 2 */
697 nz += (int)((k+1)*(k+1));
698
699 }
700 irn = data->irn = MALLOC(sizeof(int)*nz);
701 jcn = data->jcn = MALLOC(sizeof(int)*nz);
702 val = data->val = MALLOC(sizeof(double)*nz);
703 }
704 nz = 0;
705 for (i = 0; i < n_constr_nodes; i++){
706 ii = constr_nodes[i];
707 jj = ja[ia[ii]]; ll = ja[ia[ii] + 1];
708 if (jj == ll) continue; /* do not do loops */
709 dist = distance_cropped(x, dim, jj, ll);
710 dist *= dist;
711
712 k = ia[ii+1] - ia[ii];/* usually k = 2 */
713 kk = k*k;
714 irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(dist);
715 k = constr_penalty/(k*dist); kk = constr_penalty/(kk*dist);
716 for (j = ia[ii]; j < ia[ii+1]; j++){
717 irn[nz] = ii; jcn[nz] = ja[j]; val[nz++] = -k;
718 }
719 for (j = ia[ii]; j < ia[ii+1]; j++){
720 jj = ja[j];
721 irn[nz] = jj; jcn[nz] = ii; val[nz++] = -k;
722 for (l = ia[ii]; l < ia[ii+1]; l++){
723 ll = ja[l];
724 irn[nz] = jj; jcn[nz] = ll; val[nz++] = kk;
725 }
726 }
727 }
728 Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(real));
729 } else if (edge_labeling_scheme == ELSCHEME_PENALTY2 || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY2){
730 /* 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
731 1/d_jk, and to the right hand side: {0,...,average_position_of_i's neighbor if i is an edge node,...}
732 */
733 if (!irn){
734 assert((!jcn) && (!val));
735 nz = n_constr_nodes;
736 irn = data->irn = MALLOC(sizeof(int)*nz);
737 jcn = data->jcn = MALLOC(sizeof(int)*nz);
738 val = data->val = MALLOC(sizeof(double)*nz);
739 }
740 x00 = MALLOC(sizeof(real)*m*dim);
741 for (i = 0; i < m*dim; i++) x00[i] = 0;
742 nz = 0;
743 for (i = 0; i < n_constr_nodes; i++){
744 ii = constr_nodes[i];
745 jj = ja[ia[ii]]; ll = ja[ia[ii] + 1];
746 dist = distance_cropped(x, dim, jj, ll);
747 irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(dist);
748 for (j = ia[ii]; j < ia[ii+1]; j++){
749 jj = ja[j];
750 for (l = 0; l < dim; l++){
751 x00[ii*dim+l] += x[jj*dim+l];
752 }
753 }
754 for (l = 0; l < dim; l++) {
755 x00[ii*dim+l] *= constr_penalty/(dist)/(ia[ii+1] - ia[ii]);
756 }
757 }
758 Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(real));
759
760 }
761 *LL = Lc;
762 *rhs = x00;
763}
764
765real get_stress(int m, int dim, int *iw, int *jw, real *w, real *d, real *x, real scaling, void *data, int weighted){
766 int i, j;
767 real res = 0., dist;
768 /* we use the fact that d_ij = w_ij*graph_dist(i,j). Also, d_ij and x are scalinged by *scaling, so divide by it to get actual unscaled streee. */
769 for (i = 0; i < m; i++){
770 for (j = iw[i]; j < iw[i+1]; j++){
771 if (i == jw[j]) {
772 continue;
773 }
774 dist = d[j]/w[j];/* both negative*/
775 if (weighted){
776 res += -w[j]*(dist - distance(x, dim, i, jw[j]))*(dist - distance(x, dim, i, jw[j]));
777 } else {
778 res += (dist - distance(x, dim, i, jw[j]))*(dist - distance(x, dim, i, jw[j]));
779 }
780 }
781 }
782 return 0.5*res/scaling/scaling;
783
784}
785
786static void uniform_stress_augment_rhs(int m, int dim, real *x, real *y, real alpha, real M){
787 int i, j, k;
788 real dist, distij;
789 for (i = 0; i < m; i++){
790 for (j = i+1; j < m; j++){
791 dist = distance_cropped(x, dim, i, j);
792 for (k = 0; k < dim; k++){
793 distij = (x[i*dim+k] - x[j*dim+k])/dist;
794 y[i*dim+k] += alpha*M*distij;
795 y[j*dim+k] += alpha*M*(-distij);
796 }
797 }
798 }
799}
800
801static real uniform_stress_solve(SparseMatrix Lw, real alpha, int dim, real *x0, real *rhs, real tol, int maxit, int *flag){
802 Operator Ax;
803 Operator Precon;
804
805 Ax = Operator_uniform_stress_matmul(Lw, alpha);
806 Precon = Operator_uniform_stress_diag_precon_new(Lw, alpha);
807
808 return cg(Ax, Precon, Lw->m, dim, x0, rhs, tol, maxit, flag);
809
810}
811
812real StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol) {
813 SparseMatrix Lw = sm->Lw, Lwd = sm->Lwd, Lwdd = NULL;
814 int i, j, k, m, *id, *jd, *iw, *jw, idiag, flag = 0, iter = 0;
815 real *w, *dd, *d, *y = NULL, *x0 = NULL, *x00 = NULL, diag, diff = 1, *lambda = sm->lambda, res, alpha = 0., M = 0.;
816 SparseMatrix Lc = NULL;
817 real dij, dist;
818
819
820 Lwdd = SparseMatrix_copy(Lwd);
821 m = Lw->m;
822 x0 = N_GNEW(dim*m,real);
823 if (!x0) goto RETURN;
824
825 x0 = MEMCPY(x0, x, sizeof(real)*dim*m);
826 y = N_GNEW(dim*m,real);
827 if (!y) goto RETURN;
828
829 id = Lwd->ia; jd = Lwd->ja;
830 d = (real*) Lwd->a;
831 dd = (real*) Lwdd->a;
832 w = (real*) Lw->a;
833 iw = Lw->ia; jw = Lw->ja;
834
835#ifdef DEBUG_PRINT
836 if (Verbose) fprintf(stderr, "initial stress = %f\n", get_stress(m, dim, iw, jw, w, d, x, sm->scaling, sm->data, 1));
837#endif
838 /* for the additional matrix L due to the position constraints */
839 if (sm->scheme == SM_SCHEME_NORMAL_ELABEL){
840 get_edge_label_matrix(sm->data, m, dim, x, &Lc, &x00);
841 if (Lc) Lw = SparseMatrix_add(Lw, Lc);
842 } else if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){
843 alpha = ((real*) (sm->data))[0];
844 M = ((real*) (sm->data))[1];
845 }
846
847 while (iter++ < maxit_sm && diff > tol){
848#ifdef GVIEWER
849 if (Gviewer) {
850 drawScene();
851 if (iter%2 == 0) gviewer_dump_current_frame();
852 }
853#endif
854
855 if (sm->scheme != SM_SCHEME_STRESS_APPROX){
856 for (i = 0; i < m; i++){
857 idiag = -1;
858 diag = 0.;
859 for (j = id[i]; j < id[i+1]; j++){
860 if (i == jd[j]) {
861 idiag = j;
862 continue;
863 }
864
865 dist = distance(x, dim, i, jd[j]);
866 //if (d[j] >= -0.0001*dist){
867 // /* sometimes d[j] = 0 and ||x_i-x_j||=0*/
868 // dd[j] = d[j]/MAX(0.0001, dist);
869 if (d[j] == 0){
870 dd[j] = 0;
871 } else {
872 if (dist == 0){
873 dij = d[j]/w[j];/* the ideal distance */
874 /* perturb so points do not sit at the same place */
875 for (k = 0; k < dim; k++) x[jd[j]*dim+k] += 0.0001*(drand()+.0001)*dij;
876 dist = distance(x, dim, i, jd[j]);
877 }
878 dd[j] = d[j]/dist;
879
880#if 0
881 /* if two points are at the same place, we do not want a huge entry,
882 as this cause problem with CG./ In any case,
883 at thw limit d[j] == ||x[i] - x[jd[j]]||,
884 or close. */
885 if (dist < -d[j]*0.0000001){
886 dd[j] = -10000.;
887 } else {
888 dd[j] = d[j]/dist;
889 }
890#endif
891
892 }
893 diag += dd[j];
894 }
895 assert(idiag >= 0);
896 dd[idiag] = -diag;
897 }
898 /* solve (Lw+lambda*I) x = Lwdd y + lambda x0 */
899
900 SparseMatrix_multiply_dense(Lwdd, FALSE, x, FALSE, &y, FALSE, dim);
901 } else {
902 for (i = 0; i < m; i++){
903 for (j = 0; j < dim; j++){
904 y[i*dim+j] = 0;/* for stress_approx scheme, the whole rhs is calculated in stress_maxent_augment_rhs */
905 }
906 }
907 }
908
909 if (lambda){/* is there a penalty term? */
910 for (i = 0; i < m; i++){
911 for (j = 0; j < dim; j++){
912 y[i*dim+j] += lambda[i]*x0[i*dim+j];
913 }
914 }
915 }
916
917 /* additional term added to the rhs */
918 switch (sm->scheme){
919 case SM_SCHEME_NORMAL_ELABEL: {
920 for (i = 0; i < m; i++){
921 for (j = 0; j < dim; j++){
922 y[i*dim+j] += x00[i*dim+j];
923 }
924 }
925 break;
926 }
927 case SM_SCHEME_UNIFORM_STRESS:{/* this part can be done more efficiently using octree approximation */
928 uniform_stress_augment_rhs(m, dim, x, y, alpha, M);
929 break;
930 }
931#if UNIMPEMENTED
932 case SM_SCHEME_MAXENT:{
933#ifdef GVIEWER
934 if (Gviewer){
935 char *lab;
936 lab = MALLOC(sizeof(char)*100);
937 sprintf(lab,"maxent. alpha=%10.2g, iter=%d",stress_maxent_data_get_alpha(sm->data), iter);
938 gviewer_set_label(lab);
939 FREE(lab);
940 }
941#endif
942 stress_maxent_augment_rhs_fast(sm, dim, x, y, &flag);
943 if (flag) goto RETURN;
944 break;
945 }
946 case SM_SCHEME_STRESS_APPROX:{
947 stress_approx_augment_rhs(sm, dim, x, y, &flag);
948 if (flag) goto RETURN;
949 break;
950 }
951 case SM_SCHEME_STRESS:{
952#ifdef GVIEWER
953 if (Gviewer){
954 char *lab;
955 lab = MALLOC(sizeof(char)*100);
956 sprintf(lab,"pmds(k), iter=%d", iter);
957 gviewer_set_label(lab);
958 FREE(lab);
959 }
960#endif
961 }
962#endif /* UNIMPEMENTED */
963 default:
964 break;
965 }
966
967#ifdef DEBUG_PRINT
968 if (Verbose) {
969 fprintf(stderr, "stress1 = %g\n",get_stress(m, dim, iw, jw, w, d, x, sm->scaling, sm->data, 1));
970 }
971#endif
972
973 if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){
974 res = uniform_stress_solve(Lw, alpha, dim, x, y, sm->tol_cg, sm->maxit_cg, &flag);
975 } else {
976 res = SparseMatrix_solve(Lw, dim, x, y, sm->tol_cg, sm->maxit_cg, SOLVE_METHOD_CG, &flag);
977 //res = SparseMatrix_solve(Lw, dim, x, y, sm->tol_cg, 1, SOLVE_METHOD_JACOBI, &flag);
978 }
979
980 if (flag) goto RETURN;
981#ifdef DEBUG_PRINT
982 if (Verbose) fprintf(stderr, "stress2 = %g\n",get_stress(m, dim, iw, jw, w, d, y, sm->scaling, sm->data, 1));
983#endif
984 diff = total_distance(m, dim, x, y)/sqrt(vector_product(m*dim, x, x));
985#ifdef DEBUG_PRINT
986 if (Verbose){
987 fprintf(stderr, "Outer iter = %d, cg res = %g, ||x_{k+1}-x_k||/||x_k|| = %g\n",iter, res, diff);
988 }
989#endif
990
991
992 MEMCPY(x, y, sizeof(real)*m*dim);
993 }
994
995#ifdef DEBUG
996 _statistics[1] += iter-1;
997#endif
998
999#ifdef DEBUG_PRINT
1000 if (Verbose) fprintf(stderr, "iter = %d, final stress = %f\n", iter, get_stress(m, dim, iw, jw, w, d, x, sm->scaling, sm->data, 1));
1001#endif
1002
1003 RETURN:
1004 SparseMatrix_delete(Lwdd);
1005 if (Lc) {
1006 SparseMatrix_delete(Lc);
1007 SparseMatrix_delete(Lw);
1008 }
1009
1010 if (x0) FREE(x0);
1011 if (y) FREE(y);
1012 if (x00) FREE(x00);
1013 return diff;
1014
1015}
1016
1017void StressMajorizationSmoother_delete(StressMajorizationSmoother sm){
1018 if (!sm) return;
1019 if (sm->Lw) SparseMatrix_delete(sm->Lw);
1020 if (sm->Lwd) SparseMatrix_delete(sm->Lwd);
1021 if (sm->lambda) FREE(sm->lambda);
1022 if (sm->data) sm->data_deallocator(sm->data);
1023 FREE(sm);
1024}
1025
1026
1027TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int use_triangularization){
1028 TriangleSmoother sm;
1029 int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, jdiag, nz;
1030 SparseMatrix B;
1031 real *avg_dist, *lambda, *d, *w, diag_d, diag_w, dist;
1032 real s = 0, stop = 0, sbot = 0;
1033
1034 assert(SparseMatrix_is_symmetric(A, FALSE));
1035
1036 avg_dist = N_GNEW(m,real);
1037
1038 for (i = 0; i < m ;i++){
1039 avg_dist[i] = 0;
1040 nz = 0;
1041 for (j = ia[i]; j < ia[i+1]; j++){
1042 if (i == ja[j]) continue;
1043 avg_dist[i] += distance(x, dim, i, ja[j]);
1044 nz++;
1045 }
1046 assert(nz > 0);
1047 avg_dist[i] /= nz;
1048 }
1049
1050 sm = N_GNEW(1,struct TriangleSmoother_struct);
1051 sm->scaling = 1;
1052 sm->data = NULL;
1053 sm->scheme = SM_SCHEME_NORMAL;
1054 sm->tol_cg = 0.01;
1055 sm->maxit_cg = (int)sqrt((double) A->m);
1056
1057 lambda = sm->lambda = N_GNEW(m,real);
1058 for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
1059
1060 if (m > 2){
1061 if (use_triangularization){
1062 B= call_tri(m, dim, x);
1063 } else {
1064 B= call_tri2(m, dim, x);
1065 }
1066 } else {
1067 B = SparseMatrix_copy(A);
1068 }
1069
1070
1071
1072 sm->Lw = SparseMatrix_add(A, B);
1073
1074 SparseMatrix_delete(B);
1075 sm->Lwd = SparseMatrix_copy(sm->Lw);
1076 if (!(sm->Lw) || !(sm->Lwd)) {
1077 TriangleSmoother_delete(sm);
1078 return NULL;
1079 }
1080
1081 iw = sm->Lw->ia; jw = sm->Lw->ja;
1082
1083 w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
1084
1085 for (i = 0; i < m; i++){
1086 diag_d = diag_w = 0;
1087 jdiag = -1;
1088 for (j = iw[i]; j < iw[i+1]; j++){
1089 k = jw[j];
1090 if (k == i){
1091 jdiag = j;
1092 continue;
1093 }
1094 /* w[j] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
1095 w[j] = -2./(avg_dist[i]+avg_dist[k]);
1096 w[j] = -1.*/;/* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
1097 dist = pow(distance_cropped(x,dim,i,k),0.6);
1098 w[j] = 1/(dist*dist);
1099 diag_w += w[j];
1100
1101 /* d[j] = w[j]*distance(x,dim,i,k);
1102 d[j] = w[j]*(avg_dist[i] + avg_dist[k])*0.5;*/
1103 d[j] = w[j]*dist;
1104 stop += d[j]*distance(x,dim,i,k);
1105 sbot += d[j]*dist;
1106 diag_d += d[j];
1107
1108 }
1109
1110 lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
1111
1112 assert(jdiag >= 0);
1113 w[jdiag] = -diag_w + lambda[i];
1114 d[jdiag] = -diag_d;
1115 }
1116
1117 s = stop/sbot;
1118 for (i = 0; i < iw[m]; i++) d[i] *= s;
1119 sm->scaling = s;
1120
1121 FREE(avg_dist);
1122
1123 return sm;
1124}
1125
1126void TriangleSmoother_delete(TriangleSmoother sm){
1127
1128 StressMajorizationSmoother_delete(sm);
1129
1130}
1131
1132void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x){
1133
1134 StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001);
1135}
1136
1137
1138
1139
1140/* ================================ spring and spring-electrical based smoother ================ */
1141SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x){
1142 SpringSmoother sm;
1143 int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *id, *jd;
1144 int *mask, nz;
1145 real *d, *dd;
1146 real *avg_dist;
1147 SparseMatrix ID = NULL;
1148
1149 assert(SparseMatrix_is_symmetric(A, FALSE));
1150
1151 ID = ideal_distance_matrix(A, dim, x);
1152 dd = (real*) ID->a;
1153
1154 sm = N_GNEW(1,struct SpringSmoother_struct);
1155 mask = N_GNEW(m,int);
1156
1157 avg_dist = N_GNEW(m,real);
1158
1159 for (i = 0; i < m ;i++){
1160 avg_dist[i] = 0;
1161 nz = 0;
1162 for (j = ia[i]; j < ia[i+1]; j++){
1163 if (i == ja[j]) continue;
1164 avg_dist[i] += distance(x, dim, i, ja[j]);
1165 nz++;
1166 }
1167 assert(nz > 0);
1168 avg_dist[i] /= nz;
1169 }
1170
1171
1172 for (i = 0; i < m; i++) mask[i] = -1;
1173
1174 nz = 0;
1175 for (i = 0; i < m; i++){
1176 mask[i] = i;
1177 for (j = ia[i]; j < ia[i+1]; j++){
1178 k = ja[j];
1179 if (mask[k] != i){
1180 mask[k] = i;
1181 nz++;
1182 }
1183 }
1184 for (j = ia[i]; j < ia[i+1]; j++){
1185 k = ja[j];
1186 for (l = ia[k]; l < ia[k+1]; l++){
1187 if (mask[ja[l]] != i){
1188 mask[ja[l]] = i;
1189 nz++;
1190 }
1191 }
1192 }
1193 }
1194
1195 sm->D = SparseMatrix_new(m, m, nz, MATRIX_TYPE_REAL, FORMAT_CSR);
1196 if (!(sm->D)){
1197 SpringSmoother_delete(sm);
1198 return NULL;
1199 }
1200
1201 id = sm->D->ia; jd = sm->D->ja;
1202 d = (real*) sm->D->a;
1203 id[0] = 0;
1204
1205 nz = 0;
1206 for (i = 0; i < m; i++){
1207 mask[i] = i+m;
1208 for (j = ia[i]; j < ia[i+1]; j++){
1209 k = ja[j];
1210 if (mask[k] != i+m){
1211 mask[k] = i+m;
1212 jd[nz] = k;
1213 d[nz] = (avg_dist[i] + avg_dist[k])*0.5;
1214 d[nz] = dd[j];
1215 nz++;
1216 }
1217 }
1218
1219 for (j = ia[i]; j < ia[i+1]; j++){
1220 k = ja[j];
1221 for (l = ia[k]; l < ia[k+1]; l++){
1222 if (mask[ja[l]] != i+m){
1223 mask[ja[l]] = i+m;
1224 jd[nz] = ja[l];
1225 d[nz] = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
1226 d[nz] = dd[j]+dd[l];
1227 nz++;
1228 }
1229 }
1230 }
1231 id[i+1] = nz;
1232 }
1233 sm->D->nz = nz;
1234 sm->ctrl = spring_electrical_control_new();
1235 *(sm->ctrl) = *ctrl;
1236 sm->ctrl->random_start = FALSE;
1237 sm->ctrl->multilevels = 1;
1238 sm->ctrl->step /= 2;
1239 sm->ctrl->maxiter = 20;
1240
1241 FREE(mask);
1242 FREE(avg_dist);
1243 SparseMatrix_delete(ID);
1244
1245 return sm;
1246}
1247
1248
1249void SpringSmoother_delete(SpringSmoother sm){
1250 if (!sm) return;
1251 if (sm->D) SparseMatrix_delete(sm->D);
1252 if (sm->ctrl) spring_electrical_control_delete(sm->ctrl);
1253}
1254
1255
1256
1257
1258void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x){
1259 int flag = 0;
1260
1261 spring_electrical_spring_embedding(dim, A, sm->D, sm->ctrl, node_weights, x, &flag);
1262 assert(!flag);
1263
1264}
1265
1266/*=============================== end of spring and spring-electrical based smoother =========== */
1267
1268void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
1269#ifdef TIME
1270 clock_t cpu;
1271#endif
1272
1273#ifdef TIME
1274 cpu = clock();
1275#endif
1276 *flag = 0;
1277
1278 switch (ctrl->smoothing){
1279 case SMOOTHING_RNG:
1280 case SMOOTHING_TRIANGLE:{
1281 TriangleSmoother sm;
1282
1283 if (A->m > 2) { /* triangles need at least 3 nodes */
1284 if (ctrl->smoothing == SMOOTHING_RNG){
1285 sm = TriangleSmoother_new(A, dim, 0, x, FALSE);
1286 } else {
1287 sm = TriangleSmoother_new(A, dim, 0, x, TRUE);
1288 }
1289 TriangleSmoother_smooth(sm, dim, x);
1290 TriangleSmoother_delete(sm);
1291 }
1292 break;
1293 }
1294 case SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST:
1295 case SMOOTHING_STRESS_MAJORIZATION_POWER_DIST:
1296 case SMOOTHING_STRESS_MAJORIZATION_AVG_DIST:
1297 {
1298 StressMajorizationSmoother sm;
1299 int k, dist_scheme = IDEAL_AVG_DIST;
1300
1301 if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST){
1302 dist_scheme = IDEAL_GRAPH_DIST;
1303 } else if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_AVG_DIST){
1304 dist_scheme = IDEAL_AVG_DIST;
1305 } else if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_POWER_DIST){
1306 dist_scheme = IDEAL_POWER_DIST;
1307 }
1308
1309 for (k = 0; k < 1; k++){
1310 sm = StressMajorizationSmoother2_new(A, dim, 0.05, x, dist_scheme);
1311 StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001);
1312 StressMajorizationSmoother_delete(sm);
1313 }
1314 break;
1315 }
1316 case SMOOTHING_SPRING:{
1317 SpringSmoother sm;
1318 int k;
1319
1320 for (k = 0; k < 1; k++){
1321 sm = SpringSmoother_new(A, dim, ctrl, x);
1322 SpringSmoother_smooth(sm, A, node_weights, dim, x);
1323 SpringSmoother_delete(sm);
1324 }
1325
1326 break;
1327 }
1328
1329 }/* end switch between smoothing methods */
1330
1331#ifdef TIME
1332 if (Verbose) fprintf(stderr, "post processing %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
1333#endif
1334}
1335