33#define node_degree(i) (ia[(i)+1] - ia[(i)])
39 int *ia, *ja, i, j, k, l, nz;
41 double len, di, sum, sumd;
56 for (i = 0; i <
D->m; i++) mask[i] = -1;
58 for (i = 0; i <
D->m; i++){
61 for (j = ia[i]; j < ia[i+1]; j++){
62 if (i == ja[j])
continue;
65 for (j = ia[i]; j < ia[i+1]; j++){
69 for (l = ia[k]; l < ia[k+1]; l++){
70 if (mask[ja[l]] == i)
len--;
80 for (i = 0; i <
D->m; i++){
81 for (j = ia[i]; j < ia[i+1]; j++){
82 if (i == ja[j])
continue;
88 sum /= nz; sumd /= nz;
91 for (i = 0; i <
D->m; i++){
92 for (j = ia[i]; j < ia[i+1]; j++){
93 if (i == ja[j])
continue;
104 int ideal_dist_scheme){
108 int i, j, k, l, m =
A->m, *ia =
A->ia, *ja =
A->ja, *iw, *jw, *
id, *jd;
110 double *d, *w, *lambda;
111 double diag_d, diag_w,
dist,
s = 0, stop = 0, sbot = 0;
126 for (i = 0; i < m; i++) sm->
lambda[i] = lambda0;
129 double *avg_dist =
gv_calloc(m,
sizeof(
double));
131 for (i = 0; i < m ;i++){
134 for (j = ia[i]; j < ia[i+1]; j++){
135 if (i == ja[j])
continue;
144 for (i = 0; i < m; i++) mask[i] = -1;
147 for (i = 0; i < m; i++){
149 for (j = ia[i]; j < ia[i+1]; j++){
156 for (j = ia[i]; j < ia[i+1]; j++){
158 for (l = ia[k]; l < ia[k+1]; l++){
159 if (mask[ja[l]] != i){
181 for (i = 0; i < m; i++){
184 for (j = ia[i]; j < ia[i+1]; j++){
193 dist = (avg_dist[i] + avg_dist[k])*0.5;
197 fprintf(stderr,
"ideal_dist_scheme value wrong");
226 for (j = ia[i]; j < ia[i+1]; j++){
228 for (l = ia[k]; l < ia[k+1]; l++){
229 if (mask[ja[l]] != i+m){
235 dist = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
239 fprintf(stderr,
"ideal_dist_scheme value wrong");
270 lambda[i] *= (-diag_w);
272 w[nz] = -diag_w + lambda[i];
281 for (i = 0; i < nz; i++) d[i] *=
s;
298 double stop = 0, sbot = 0;
303 bool has_nonzero =
false;
304 for (
int i = 0; i < m *
dim; i++) {
311 for (
int i = 0; i < m*
dim; i++) x[i] = 72*
drand();
314 const int *
const ia =
A->ia;
315 const int *
const ja =
A->ja;
316 const double *
const a =
A->a;
335 int *
const iw = sm->
Lw->
ia;
336 int *
const jw = sm->
Lw->
ja;
337 int *
const id = sm->
Lwd->
ia;
338 int *
const jd = sm->
Lwd->
ja;
339 double *
const w = sm->
Lw->
a;
340 double *
const d = sm->
Lwd->
a;
344 for (
int i = 0; i < m; i++){
347 for (
int j = ia[i]; j < ia[i+1]; j++){
352 const double dist = a[j];
367 lambda[i] *= (-diag_w);
368 w[nz] = -diag_w + lambda[i];
377 const double s = stop / sbot;
382 for (
int i = 0; i < nz; i++) d[i] *=
s;
393 double total = 0,
dist = 0;
396 for (i = 0; i < m; i++){
398 for (j = 0; j <
dim; j++){
419 int edge_labeling_scheme =
data->edge_labeling_scheme;
420 int n_constr_nodes =
data->n_constr_nodes;
421 int *constr_nodes =
data->constr_nodes;
423 int *ia = A_constr->
ia, *ja = A_constr->
ja, ii, jj, nz, l, ll, i, j;
424 int *irn =
data->irn, *jcn =
data->jcn;
425 double *val =
data->val,
dist, kk, k;
428 double constr_penalty =
data->constr_penalty;
447 assert((!jcn) && (!val));
449 for (i = 0; i < n_constr_nodes; i++){
450 ii = constr_nodes[i];
451 k = ia[ii+1] - ia[ii];
452 nz += (int)((k+1)*(k+1));
460 for (i = 0; i < n_constr_nodes; i++){
461 ii = constr_nodes[i];
462 jj = ja[ia[ii]]; ll = ja[ia[ii] + 1];
463 if (jj == ll)
continue;
467 k = ia[ii+1] - ia[ii];
469 irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(
dist);
470 k = constr_penalty/(k*
dist); kk = constr_penalty/(kk*
dist);
471 for (j = ia[ii]; j < ia[ii+1]; j++){
472 irn[nz] = ii; jcn[nz] = ja[j]; val[nz++] = -k;
474 for (j = ia[ii]; j < ia[ii+1]; j++){
476 irn[nz] = jj; jcn[nz] = ii; val[nz++] = -k;
477 for (l = ia[ii]; l < ia[ii+1]; l++){
479 irn[nz] = jj; jcn[nz] = ll; val[nz++] = kk;
489 assert((!jcn) && (!val));
497 for (i = 0; i < n_constr_nodes; i++){
498 ii = constr_nodes[i];
499 jj = ja[ia[ii]]; ll = ja[ia[ii] + 1];
501 irn[nz] = ii; jcn[nz] = ii; val[nz++] = constr_penalty/(
dist);
502 for (j = ia[ii]; j < ia[ii+1]; j++){
504 for (l = 0; l <
dim; l++){
505 x00[ii*
dim+l] += x[jj*
dim+l];
508 for (l = 0; l <
dim; l++) {
509 x00[ii*
dim+l] *= constr_penalty/(
dist)/(ia[ii+1] - ia[ii]);
520 double *d,
double *x,
double scaling) {
522 double res = 0.,
dist;
525 for (i = 0; i < m; i++){
526 for (j = iw[i]; j < iw[i+1]; j++){
534 return 0.5*res/scaling/scaling;
540 int i, j, k, m, *
id, *jd, *iw, *jw, idiag, iter = 0;
541 double *w, *dd, *d, *y =
NULL, *x0 =
NULL, *x00 =
NULL, diag, diff = 1, *lambda = sm->
lambda;
545 const double tol = 0.001;
550 x0 = calloc(
dim*m,
sizeof(
double));
553 memcpy(x0, x,
sizeof(
double)*
dim*m);
554 y = calloc(
dim*m,
sizeof(
double));
557 id = Lwd->ia; jd = Lwd->ja;
561 iw = Lw->
ia; jw = Lw->
ja;
575 while (iter++ < maxit_sm && diff >
tol){
577 for (i = 0; i < m; i++){
580 for (j =
id[i]; j <
id[i+1]; j++){
593 for (k = 0; k <
dim; k++) x[jd[j]*
dim+k] += 0.0001*(
drand()+.0001)*dij;
608 for (i = 0; i < m; i++){
609 for (j = 0; j <
dim; j++){
610 y[i*
dim+j] += lambda[i]*x0[i*
dim+j];
618 for (i = 0; i < m; i++){
619 for (j = 0; j <
dim; j++){
620 y[i*
dim+j] += x00[i*
dim+j];
643 fprintf(stderr,
"Outer iter = %d, cg res = %g, ||x_{k+1}-x_k||/||x_k|| = %g\n",iter, res, diff);
648 memcpy(x, y,
sizeof(
double)*m*
dim);
652 _statistics[1] += iter-1;
656 if (
Verbose) fprintf(stderr,
"iter = %d, final stress = %f\n", iter,
get_stress(m,
dim, iw, jw, w, d, x, sm->
scaling));
683 bool use_triangularization) {
684 int i, j, k, m =
A->m, *ia =
A->ia, *ja =
A->ja, *iw, *jw, jdiag, nz;
686 double *d, *w, diag_d, diag_w,
dist;
687 double s = 0, stop = 0, sbot = 0;
691 double *avg_dist =
gv_calloc(m,
sizeof(
double));
693 for (i = 0; i < m ;i++){
696 for (j = ia[i]; j < ia[i+1]; j++){
697 if (i == ja[j])
continue;
715 if (use_triangularization){
730 if (!sm->
Lw || !sm->
Lwd) {
740 for (i = 0; i < m; i++){
743 for (j = iw[i]; j < iw[i+1]; j++){
765 lambda[i] *= (-diag_w);
768 w[jdiag] = -diag_w + lambda[i];
773 for (i = 0; i < iw[m]; i++) d[i] *=
s;
793 int i, j, k, l, m =
A->m, *ia =
A->ia, *ja =
A->ja, *
id, *jd;
806 double *avg_dist =
gv_calloc(m,
sizeof(
double));
808 for (i = 0; i < m ;i++){
811 for (j = ia[i]; j < ia[i+1]; j++){
812 if (i == ja[j])
continue;
821 for (i = 0; i < m; i++) mask[i] = -1;
824 for (i = 0; i < m; i++){
826 for (j = ia[i]; j < ia[i+1]; j++){
833 for (j = ia[i]; j < ia[i+1]; j++){
835 for (l = ia[k]; l < ia[k+1]; l++){
836 if (mask[ja[l]] != i){
845 assert(sm->
D !=
NULL);
847 id = sm->
D->
ia; jd = sm->
D->
ja;
852 for (i = 0; i < m; i++){
854 for (j = ia[i]; j < ia[i+1]; j++){
859 d[nz] = (avg_dist[i] + avg_dist[k])*0.5;
865 for (j = ia[i]; j < ia[i+1]; j++){
867 for (l = ia[k]; l < ia[k+1]; l++){
868 if (mask[ja[l]] != i+m){
871 d[nz] = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
965 if (
Verbose) fprintf(stderr,
"post processing %f\n",((
double) (clock() - cpu)) / CLOCKS_PER_SEC);
void SparseMatrix_multiply_dense(SparseMatrix A, const double *v, double *res, int dim)
bool SparseMatrix_is_symmetric(SparseMatrix A, bool test_pattern_symmetry_only)
SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, int *irn, int *jcn, void *val0, int type, size_t sz)
void SparseMatrix_delete(SparseMatrix A)
SparseMatrix SparseMatrix_copy(SparseMatrix A)
SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B)
SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format)
Memory allocation wrappers that exit on failure.
static void * gv_calloc(size_t nmemb, size_t size)
static void * gv_alloc(size_t size)
SparseMatrix call_tri2(int n, int dim, double *xx)
SparseMatrix call_tri(int n, double *x)
static NORETURN void graphviz_exit(int status)
static double dist(int dim, double *x, double *y)
double distance(double *x, int dim, int i, int j)
double vector_product(int n, double *x, double *y)
double distance_cropped(double *x, int dim, int i, int j)
static double len(glCompPoint p)
Arithmetic helper functions.
static bool is_exactly_equal(double a, double b)
are two values precisely the same?
@ ELSCHEME_STRAIGHTLINE_PENALTY2
@ ELSCHEME_STRAIGHTLINE_PENALTY
static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, double *x, SparseMatrix *LL, double **rhs)
StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, double *x)
static double total_distance(int m, int dim, double *x, double *y)
static SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, double *x)
void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm)
TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, double *x, bool use_triangularization)
double StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, double *x, int maxit_sm)
void SpringSmoother_delete(SpringSmoother sm)
double SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, double *x, int maxit_sm)
void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, double *x)
void TriangleSmoother_smooth(TriangleSmoother sm, int dim, double *x)
void StressMajorizationSmoother_delete(StressMajorizationSmoother sm)
SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, double *x)
StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, double lambda0, double *x, int ideal_dist_scheme)
static UNUSED double get_stress(int m, int dim, int *iw, int *jw, double *w, double *d, double *x, double scaling)
void TriangleSmoother_delete(TriangleSmoother sm)
void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, int dim, double *x)
@ SM_SCHEME_NORMAL_ELABEL
#define TriangleSmoother_struct
double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, double maxit)
void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control *ctrl, double *x, int *flag)
@ SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST
@ SMOOTHING_STRESS_MAJORIZATION_POWER_DIST
@ SMOOTHING_STRESS_MAJORIZATION_AVG_DIST
spring_electrical_control ctrl
SparseMatrix Lwd
the laplacian like matrix with offdiag = -scaling × dᵢⱼ ÷ wᵢⱼ. RHS in stress majorization = Lwd....
void(* data_deallocator)(void *)
SparseMatrix Lw
the weighted laplacian. with offdiag = -1 ÷ wᵢⱼ
bool random_start
whether to apply SE from a random layout, or from exisiting layout
graphs, nodes and edges info: Agraphinfo_t, Agnodeinfo_t and Agedgeinfo_t
abstraction for squashing compiler warnings for unused symbols