Graphviz  2.41.20170921.2350
uniform_stress.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
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 "general.h"
15 #include "SparseMatrix.h"
16 #include "spring_electrical.h"
17 #include "post_process.h"
18 #include "sparse_solve.h"
19 #include <time.h>
20 #include "uniform_stress.h"
21
22 /* uniform stress solves the model:
23
24 \Sum_{i<->j} (||x_i-x_j||-1)^2 + alpha * \Sum_{i!=j} (||x_i-x_j||^2-M)^2
25
26 This is somewhat similar to the binary stress model
27
28 */
29
32  int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
33  int nz;
34  real *d, *w, *a = (real*) A->a;
35  real diag_d, diag_w, dist, epsilon = 0.01;
36
38
39  sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct));
40  sm->data = NULL;
42  sm->lambda = NULL;
43  sm->data = MALLOC(sizeof(real)*2);
44  ((real*) sm->data)[0] = alpha;
45  ((real*) sm->data)[1] = M;
46  sm->data_deallocator = FREE;
47  sm->tol_cg = 0.01;
48  sm->maxit_cg = (int)sqrt((double) A->m);
49
50  /* Lw and Lwd have diagonals */
51  sm->Lw = SparseMatrix_new(m, m, A->nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
52  sm->Lwd = SparseMatrix_new(m, m, A->nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
53  iw = sm->Lw->ia; jw = sm->Lw->ja;
54  id = sm->Lwd->ia; jd = sm->Lwd->ja;
55  w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
56
57  if (!(sm->Lw) || !(sm->Lwd)) {
59  return NULL;
60  }
61
62  iw = sm->Lw->ia; jw = sm->Lw->ja;
63  id = sm->Lwd->ia; jd = sm->Lwd->ja;
64  w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
65  iw[0] = id[0] = 0;
66
67  nz = 0;
68  for (i = 0; i < m; i++){
69  diag_d = diag_w = 0;
70  for (j = ia[i]; j < ia[i+1]; j++){
71  k = ja[j];
72  if (k != i){
73  dist = MAX(ABS(a[j]), epsilon);
74  jd[nz] = jw[nz] = k;
75  w[nz] = -1/(dist*dist);
76  w[nz] = -1.;
77  d[nz] = w[nz]*dist;
78  diag_w += w[nz];
79  diag_d += d[nz];
80  nz++;
81  }
82  }
83  jd[nz] = jw[nz] = i;
84  w[nz] = -diag_w;
85  d[nz] = -diag_d;
86  nz++;
87
88  iw[i+1] = nz;
89  id[i+1] = nz;
90
91  }
92
93  sm->Lw->nz = nz;
94  sm->Lwd->nz = nz;
95
96  return sm;
97 }
98
99
101
103
104 }
105
107
108  return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, 0.001);
109
110 }
111
113  /* get a distance matrix from a graph, at the moment we just symmetrize the matrix. At the moment if the matrix is not real,
114  we just assume distance of 1 among edges. Then we apply scaling to the entire matrix */
115  SparseMatrix B;
116  real *val;
117  int i;
118
119  if (A->type == MATRIX_TYPE_REAL){
121  } else {
123  }
124  val = (real*) B->a;
125  if (scaling != 1) for (i = 0; i < B->nz; i++) val[i] *= scaling;
126  return B;
127 }
128
129 extern void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x);
130
131 void uniform_stress(int dim, SparseMatrix A, real *x, int *flag){
133  real lambda0 = 10.1, M = 100, scaling = 1.;
134  int maxit = 300, samepoint = TRUE, i, k, n = A->m;
135  SparseMatrix B = NULL;
136
137  *flag = 0;
138
139  /* just set random initial for now */
140  for (i = 0; i < dim*n; i++) {
141  x[i] = M*drand();
142  }
143
144  /* make sure x is not all at the same point */
145  for (i = 1; i < n; i++){
146  for (k = 0; k < dim; k++) {
147  if (ABS(x[0*dim+k] - x[i*dim+k]) > MACHINEACC){
148  samepoint = FALSE;
149  i = n;
150  break;
151  }
152  }
153  }
154
155  if (samepoint){
156  srand(1);
157 #ifdef DEBUG_PRINT
158  fprintf(stderr,"input coordinates to uniform_stress are the same, use random coordinates as initial input");
159 #endif
160  for (i = 0; i < dim*n; i++) x[i] = M*drand();
161  }
162
163  B = get_distance_matrix(A, scaling);
165
166  sm = UniformStressSmoother_new(dim, B, x, 1000000*lambda0, M, flag);
167  UniformStressSmoother_smooth(sm, dim, x, maxit);
169
170  sm = UniformStressSmoother_new(dim, B, x, 10000*lambda0, M, flag);
171  UniformStressSmoother_smooth(sm, dim, x, maxit);
173
174  sm = UniformStressSmoother_new(dim, B, x, 100*lambda0, M, flag);
175  UniformStressSmoother_smooth(sm, dim, x, maxit);
177
178  sm = UniformStressSmoother_new(dim, B, x, lambda0, M, flag);
179  UniformStressSmoother_smooth(sm, dim, x, maxit);
181
182  scale_to_box(0,0,7*70,10*70,A->m,dim,x);;
183
185
186 }
void uniform_stress(int dim, SparseMatrix A, real *x, int *flag)
#define MAX(a, b)
Definition: agerror.c:17
void StressMajorizationSmoother_delete(StressMajorizationSmoother sm)
#define ABS(a)
Definition: arith.h:45
double xmax
Definition: geometry.c:20
real drand()
Definition: general.c:52
void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x)
Definition: general.c:312
#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
double xmin
Definition: geometry.c:20
SparseMatrix get_distance_matrix(SparseMatrix A, real scaling)
double ymax
Definition: geometry.c:20
double ymin
Definition: geometry.c:20
int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only)
Definition: SparseMatrix.c:178
void UniformStressSmoother_delete(UniformStressSmoother sm)
Definition: post_process.h:19
int
Definition: grammar.c:1264
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
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only)
Definition: SparseMatrix.c:151
UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, real *x, real alpha, real M, int *flag)
void SparseMatrix_delete(SparseMatrix A)
Definition: SparseMatrix.c:411
if(aagss+aagstacksize-1<=aagssp)
Definition: grammar.c:1332
#define NULL
Definition: logic.h:39
#define alpha
Definition: shapes.c:3902
real UniformStressSmoother_smooth(UniformStressSmoother sm, int dim, real *x, int maxit_sm)
double dist(Site *s, Site *t)
Definition: site.c:41
#define MACHINEACC
Definition: general.h:120
#define FALSE
Definition: cgraph.h:35
#define TRUE
Definition: cgraph.h:38
#define real
Definition: general.h:34