Graphviz  2.41.20170921.2350
post_process.c
Go to the documentation of this file.
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
34 static 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 
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 
85 void 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 
164 real 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){
192  wij = 1/(dij*dij);
193  break;
195  wij = 1/(dij);
196  break;
197  break;
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 
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 
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 
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  */
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;
332 
334 
335  ID = ideal_distance_matrix(A, dim, x);
336 
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)) {
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);
510  return sm;
511 }
512 
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  */
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 
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)) {
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){
574  if (dist*dist == 0){
575  w[nz] = -100000;
576  } else {
577  w[nz] = -1/(dist*dist);
578  }
579  break;
581  if (dist*dist == 0){
582  w[nz] = -100000;
583  } else {
584  w[nz] = -1/(dist);
585  }
586  break;
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 
635 static 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 
654 }
655 
656 
658 
659  return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, tol);
660 
661 
662 }
663 static 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 
765 real 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 
786 static 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 
801 static 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 
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){
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  }
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 
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 
1027 TriangleSmoother 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 
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 
1075  sm->Lwd = SparseMatrix_copy(sm->Lw);
1076  if (!(sm->Lw) || !(sm->Lwd)) {
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 
1127 
1129 
1130 }
1131 
1133 
1134  StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001);
1135 }
1136 
1137 
1138 
1139 
1140 /* ================================ spring and spring-electrical based smoother ================ */
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 
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)){
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;
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 
1250  if (!sm) return;
1251  if (sm->D) SparseMatrix_delete(sm->D);
1253 }
1254 
1255 
1256 
1257 
1258 void 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 
1268 void 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);
1291  }
1292  break;
1293  }
1297  {
1299  int k, dist_scheme = IDEAL_AVG_DIST;
1300 
1302  dist_scheme = IDEAL_GRAPH_DIST;
1304  dist_scheme = IDEAL_AVG_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);
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);
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 }
#define MAX(a, b)
Definition: agerror.c:17
void StressMajorizationSmoother_delete(StressMajorizationSmoother sm)
Definition: post_process.h:40
SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real *x)
Definition: post_process.c:251
spring_electrical_control spring_electrical_control_new()
double xmax
Definition: geometry.c:20
real drand()
Definition: general.c:52
void spring_electrical_control_delete(spring_electrical_control ctrl)
Definition: post_process.h:82
#define MIN(a, b)
Definition: arith.h:35
SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x)
Definition: post_process.h:19
SparseMatrix SparseMatrix_copy(SparseMatrix A)
#define FREE
Definition: PriorityQueue.c:23
SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format)
Definition: SparseMatrix.c:384
#define assert(x)
Definition: cghdr.h:47
void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x)
double xmin
Definition: geometry.c:20
Definition: post_process.h:40
real distance(real *x, int dim, int i, int j)
Definition: post_process.h:40
void SparseMatrix_multiply_dense(SparseMatrix A, int ATransposed, real *v, int vTransposed, real **res, int res_transposed, int dim)
StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int weighting_scheme, int scale_initial_coord)
Definition: post_process.c:513
#define node_degree(i)
Definition: post_process.c:31
int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only)
Definition: SparseMatrix.c:178
#define TriangleSmoother_struct
Definition: post_process.h:47
Definition: post_process.h:19
Definition: xdot.h:148
int
Definition: grammar.c:1264
void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
void SpringSmoother_delete(SpringSmoother sm)
real get_stress(int m, int dim, int *iw, int *jw, real *w, real *d, real *x, real scaling, void *data, int weighted)
Definition: post_process.c:765
void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm)
Definition: post_process.c:652
StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, real lambda0, real *x, int ideal_dist_scheme)
Definition: post_process.c:320
int SparseMatrix_distance_matrix(SparseMatrix D0, int weighted, real **dist0)
real StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol)
Definition: post_process.c:812
#define MALLOC
Definition: PriorityQueue.c:21
real vector_product(int n, real *x, real *y)
Definition: general.c:101
Definition: post_process.h:82
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only)
Definition: SparseMatrix.c:151
Definition: post_process.h:19
TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int use_triangularization)
Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha)
Definition: sparse_solve.c:100
spring_electrical_control ctrl
Definition: post_process.h:61
void vector_ordering(int n, real *v, int **p, int ascending)
Definition: general.c:210
void SparseMatrix_delete(SparseMatrix A)
Definition: SparseMatrix.c:411
Definition: grammar.c:79
#define MEMCPY
Definition: PriorityQueue.c:24
real SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol)
Definition: post_process.c:657
if(aagss+aagstacksize-1<=aagssp)
Definition: grammar.c:1332
SparseMatrix call_tri(int n, int dim, real *x)
Definition: call_tri.c:21
void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag)
void TriangleSmoother_delete(TriangleSmoother sm)
#define NULL
Definition: logic.h:39
Definition: post_process.h:19
real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit, int *flag)
Definition: sparse_solve.c:227
#define GNEW(t)
Definition: memory.h:37
real get_full_stress(SparseMatrix A, int dim, real *x, int weighting_scheme)
Definition: post_process.h:19
void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x)
SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz)
Definition: SparseMatrix.c:957
real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag)
Definition: sparse_solve.c:296
EXTERN unsigned char Verbose
Definition: globals.h:64
Definition: post_process.h:19
#define top(sp)
Definition: stack.h:35
for(;;)
Definition: grammar.c:1846
Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha)
Definition: sparse_solve.c:56
real vector_median(int n, real *x)
Definition: general.c:21
#define alpha
Definition: shapes.c:3902
real vector_percentile(int n, real *x, real y)
Definition: general.c:35
#define ID
void dump_distance_edge_length(char *outfile, SparseMatrix A, int dim, real *x)
SparseMatrix call_tri2(int n, int dim, real *xx)
Definition: call_tri.c:69
double dist(Site *s, Site *t)
Definition: site.c:41
real distance_cropped(real *x, int dim, int i, int j)
SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B)
Definition: SparseMatrix.c:966
#define N_GNEW(n, t)
Definition: agxbuf.c:20
#define FALSE
Definition: cgraph.h:35
Definition: legal.c:60
Definition: post_process.h:82
#define TRUE
Definition: cgraph.h:38
#define real
Definition: general.h:34