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 "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
26This is somewhat similar to the binary stress model
27
28*/
29
30UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, real *x, real alpha, real M, int *flag){
31 UniformStressSmoother sm;
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
37 assert(SparseMatrix_is_symmetric(A, FALSE));
38
39 sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct));
40 sm->data = NULL;
41 sm->scheme = SM_SCHEME_UNIFORM_STRESS;
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)) {
58 StressMajorizationSmoother_delete(sm);
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
100void UniformStressSmoother_delete(UniformStressSmoother sm){
101
102 StressMajorizationSmoother_delete(sm);
103
104}
105
106real UniformStressSmoother_smooth(UniformStressSmoother sm, int dim, real *x, int maxit_sm) {
107
108 return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, 0.001);
109
110}
111
112SparseMatrix get_distance_matrix(SparseMatrix A, real scaling){
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){
120 B = SparseMatrix_symmetrize(A, FALSE);
121 } else {
122 B = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
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
129extern void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x);
130
131void uniform_stress(int dim, SparseMatrix A, real *x, int *flag){
132 UniformStressSmoother sm;
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);
164 assert(SparseMatrix_is_symmetric(B, FALSE));
165
166 sm = UniformStressSmoother_new(dim, B, x, 1000000*lambda0, M, flag);
167 UniformStressSmoother_smooth(sm, dim, x, maxit);
168 UniformStressSmoother_delete(sm);
169
170 sm = UniformStressSmoother_new(dim, B, x, 10000*lambda0, M, flag);
171 UniformStressSmoother_smooth(sm, dim, x, maxit);
172 UniformStressSmoother_delete(sm);
173
174 sm = UniformStressSmoother_new(dim, B, x, 100*lambda0, M, flag);
175 UniformStressSmoother_smooth(sm, dim, x, maxit);
176 UniformStressSmoother_delete(sm);
177
178 sm = UniformStressSmoother_new(dim, B, x, lambda0, M, flag);
179 UniformStressSmoother_smooth(sm, dim, x, maxit);
180 UniformStressSmoother_delete(sm);
181
182 scale_to_box(0,0,7*70,10*70,A->m,dim,x);;
183
184 SparseMatrix_delete(B);
185
186}
187