1#include "general.h"
2#include "SparseMatrix.h"
3#include "spring_electrical.h"
4#include "post_process.h"
5#include "stress_model.h"
6
7void stress_model_core(int dim, SparseMatrix B, real **x, int edge_len_weighted, int maxit_sm, real tol, int *flag){
8 int m;
9 SparseStressMajorizationSmoother sm;
10 real lambda = 0;
11 /*int maxit_sm = 1000, i; tol = 0.001*/
12 int i;
13 SparseMatrix A = B;
14
15 if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){
16 if (A->type == MATRIX_TYPE_REAL){
17 A = SparseMatrix_symmetrize(A, FALSE);
18 A = SparseMatrix_remove_diagonal(A);
19 } else {
20 A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
21 }
22 }
23 A = SparseMatrix_remove_diagonal(A);
24
25 *flag = 0;
26 m = A->m;
27 if (!x) {
28 *x = MALLOC(sizeof(real)*m*dim);
29 srand(123);
30 for (i = 0; i < dim*m; i++) (*x)[i] = drand();
31 }
32
33 if (edge_len_weighted){
34 sm = SparseStressMajorizationSmoother_new(A, dim, lambda, *x, WEIGHTING_SCHEME_SQR_DIST, TRUE);/* do not under weight the long distances */
35 //sm = SparseStressMajorizationSmoother_new(A, dim, lambda, *x, WEIGHTING_SCHEME_INV_DIST, TRUE);/* do not under weight the long distances */
36 } else {
37 sm = SparseStressMajorizationSmoother_new(A, dim, lambda, *x, WEIGHTING_SCHEME_NONE, TRUE);/* weight the long distances */
38 }
39
40 if (!sm) {
41 *flag = -1;
42 goto RETURN;
43 }
44
45
46 sm->tol_cg = 0.1; /* we found that there is no need to solve the Laplacian accurately */
47 sm->scheme = SM_SCHEME_STRESS;
48 SparseStressMajorizationSmoother_smooth(sm, dim, *x, maxit_sm, tol);
49 for (i = 0; i < dim*m; i++) {
50 (*x)[i] /= sm->scaling;
51 }
52 SparseStressMajorizationSmoother_delete(sm);
53
54 RETURN:
55 if (A != B) SparseMatrix_delete(A);
56
57}
58
59
60#ifdef GVIEWER
61#include "gviewer.h"
62#include "get_ps.h"
63struct stress_model_data {
64 int dim;
65 SparseMatrix D;
66 real **x;
67 int edge_len_weighted;
68 int maxit_sm;
69 real tol;
70 int *flag;
71};
72
73void stress_model_gv(void* data){
74 struct stress_model_data* d;
75
76 d = (struct stress_model_data*) data;
77 return stress_model_core(d->dim, d->D, d->x, d->edge_len_weighted, d->maxit_sm, d->tol, d->flag);
78}
79void stress_model(int dim, SparseMatrix A, SparseMatrix D, real **x, int edge_len_weighted, int maxit_sm, real tol, int *flag){
80 struct stress_model_data data = {dim, D, x, edge_len_weighted, maxit_sm, tol, flag};
81
82 int argcc = 1;
83 char **argvv;
84
85 if (!Gviewer) return stress_model_core(dim, D, x, edge_len_weighted, maxit_sm, tol, flag);
86 argcc = 1;
87 argvv = malloc(sizeof(char*)*argcc);
88 argvv[0] = malloc(sizeof(char));
89 argvv[0][0] = '1';
90 // gviewer_init(&argcc, argvv, 0.1, 20, 60, 2*1010, 2*770, A, dim, *x, &(data), stress_model_gv);
91 gviewer_set_edge_color_scheme(COLOR_SCHEME_NO);
92 gviewer_toggle_bgcolor();
93 gviewer_init(&argcc, argvv, 0.1, 20, 60, 720, 720, A, dim, *x, &(data), stress_model_gv);
94 free(argvv);
95
96}
97#else
98void stress_model(int dim, SparseMatrix A, SparseMatrix D, real **x, int edge_len_weighted, int maxit_sm, real tol, int *flag){
99 stress_model_core(dim, D, x, edge_len_weighted, maxit_sm, tol, flag);
100}
101#endif
102
103