#define VERSION "VVVV"
#define BUILDDATE "DDDD"
+#define DIGCOLA 1
+#define IPSEPCOLA 1
+/* #define MOSEK 1 */
+
#ifndef DEFAULT_FONTPATH
# ifdef _UWIN
# define DEFAULT_FONTPATH "/win/fonts"
+# $Id$ $Revision$
## Process this file with automake to produce Makefile.in
AM_CPPFLAGS = \
-I$(top_srcdir)/lib/pack \
-I$(top_srcdir)/lib/pathplan \
-I$(top_srcdir)/lib/graph \
- -I$(top_srcdir)/lib/cdt \
- -I$(top_srcdir)/lib/libvpsc
+ -I$(top_srcdir)/lib/cdt
noinst_HEADERS = adjust.h edges.h geometry.h heap.h hedges.h info.h mem.h \
neato.h poly.h neatoprocs.h simple.h site.h voronoi.h \
bfs.h closest.h conjgrad.h defs.h dijkstra.h embed_graph.h kkutils.h \
- matrix_ops.h pca.h stress.h quad_prog_solver.h digcola.h \
- quad_prog_vpsc.h mosek_quad_solve.h
+ matrix_ops.h pca.h stress.h quad_prog_solver.h digcola.h
noinst_LTLIBRARIES = libneatogen.la
libneatogen_la_SOURCES = adjust.c circuit.c edges.c find_ints.c geometry.c \
voronoi.c stress.c kkutils.c matrix_ops.c embed_graph.c dijkstra.c \
conjgrad.c pca.c closest.c bfs.c constraint.c quad_prog_solve.c \
smart_ini_x.c constrained_majorization.c opt_arrangement.c \
- compute_hierarchy.c \
- quad_prog_vpsc.c constrained_majorization_ipsep.c \
- mosek_quad_solve.c
+ compute_hierarchy.c
EXTRA_DIST = Makefile.old
INCS = -I. -I$(ROOT) \
-I../common \
-I../pack \
+ -I../vpsc \
-I../gvc \
-I../cdt \
-I../pathplan \
DEFINES = -DHAVE_CONFIG_H
HDRS = adjust.h edges.h geometry.h heap.h hedges.h info.h mem.h \
- neato.h neatoprocs.h poly.h simple.h site.h voronoi.h
+ neato.h neatoprocs.h poly.h simple.h site.h voronoi.h \
+ quad_prog_vpsc.h mosek_quad_solve.h
NOBJS = stress.o kkutils.o pca.o matrix_ops.o embed_graph.o dijkstra.o \
conjgrad.o closest.o bfs.o compute_hierarchy.o opt_arrangement.o \
smart_ini_x.o constrained_majorization.o quad_prog_solve.o
+POBJS = quad_prog_vpsc.o constrained_majorization_ipsep.o mosek_quad_solve.o
OBJS = adjust.o circuit.o edges.o find_ints.o geometry.o heap.o hedges.o \
info.o neatoinit.o intersect.o legal.o lu.o matinv.o memory.o poly.o \
printvis.o site.o solve.o neatosplines.o stuff.o voronoi.o \
- $(NOBJS) constraint.o
+ $(NOBJS) $(POBJS) constraint.o
libneatogen.a : $(OBJS)
$(RM) libneatogen.a
double conj_tol=tolerance_cg; /* tolerance of Conjugate Gradient */
float ** unpackedLap = NULL;
CMajEnv *cMajEnv = NULL;
- clock_t start_time;
double y_0;
int length;
DistType diameter;
/* Additionally, we never consider gaps smaller than 'abs_tol'*<avg_gap> */
double relative_tol=levels_sep_tol;
int *ordering=NULL, *levels=NULL;
- double hierarchy_spread;
float constant_term;
int count;
double degree;
x = d_coords[0]; y = d_coords[1];
compute_y_coords(graph, n, y, n);
- hierarchy_spread = compute_hierarchy(graph, n, abs_tol, relative_tol, y, &ordering, &levels, &num_levels);
+ compute_hierarchy(graph, n, abs_tol, relative_tol, y, &ordering, &levels, &num_levels);
if (num_levels<=1) {
/* no hierarchy found, use faster algorithm */
return stress_majorization_kD_mkernel(graph, n, nedges_graph, d_coords, dim, smart_ini, model, maxi);
}
else {
initLayout(graph, n, dim, d_coords);
- hierarchy_spread = compute_hierarchy(graph, n, abs_tol, relative_tol, NULL, &ordering, &levels, &num_levels);
+ compute_hierarchy(graph, n, abs_tol, relative_tol, NULL, &ordering, &levels, &num_levels);
}
if (n == 1) return 0;
old_stress=DBL_MAX; /* at least one iteration */
- start_time = clock();
-
unpackedLap = unpackMatrix(lap2, n);
cMajEnv=initConstrainedMajorization(lap2, n, ordering, levels, num_levels);
+/* vim:set shiftwidth=4 ts=8: */
+
+/**
+ *
+ * Authors:
+ * Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * This version is released under the CPL (Common Public License) with
+ * the Graphviz distribution.
+ * A version is also available under the LGPL as part of the Adaptagrams
+ * project: http://sourceforge.net/projects/adaptagrams.
+ * If you make improvements or bug fixes to this code it would be much
+ * appreciated if you could also contribute those changes back to the
+ * Adaptagrams repository.
+ */
+
/**********************************************************
* Based on constrained_majorization.c
*
**********************************************************/
#include "digcola.h"
-#ifdef DIGCOLA
+#ifdef IPSEPCOLA
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define localConstrMajorIterations 1000
-int
-stress_majorization_cola(
- vtx_data* graph, /* Input graph in sparse representation */
+int stress_majorization_cola(
+ vtx_data * graph, /* Input graph in sparse representation */
int n, /* Number of nodes */
- int nedges_graph, /* Number of edges */
- double** d_coords, /* Coordinates of nodes (output layout) */
+ int nedges_graph, /* Number of edges */
+ double **d_coords, /* Coordinates of nodes (output layout) */
int dim, /* Dimemsionality of layout */
- int model, /* difference model */
- int maxi, /* max iterations */
- ipsep_options* opt
-)
+ int model, /* difference model */
+ int maxi, /* max iterations */
+ ipsep_options * opt)
{
- int iterations = 0; /* Output: number of iteration of the process */
+ int iterations = 0; /* Output: number of iteration of the process */
/*************************************************
** Computation of full, dense, unrestricted k-D **
** This function imposes HIERARCHY CONSTRAINTS **
*************************************************/
- int i,j,k;
- float * lap1 = NULL;
- float * dist_accumulator = NULL;
- float * tmp_coords = NULL;
- float ** b = NULL;
- double * degrees = NULL;
- float * lap2=NULL;
- int lap_length;
- float * f_storage=NULL;
- float ** coords=NULL;
- int orig_n=n;
-
- //double conj_tol=tolerance_cg; /* tolerance of Conjugate Gradient */
- CMajEnvVPSC *cMajEnvHor = NULL;
- CMajEnvVPSC *cMajEnvVrt = NULL;
- clock_t start_time;
- double y_0;
- int length;
- DistType diameter;
- float * Dij=NULL;
- float constant_term;
- int count;
- double degree;
- int step;
- float val;
- double old_stress, new_stress=0;
- bool converged;
- int len;
- double nsizeScale=0;
- float maxEdgeLen=0;
+ int i, j, k;
+ float *lap1 = NULL;
+ float *dist_accumulator = NULL;
+ float *tmp_coords = NULL;
+ float **b = NULL;
+ double *degrees = NULL;
+ float *lap2 = NULL;
+ int lap_length;
+ float *f_storage = NULL;
+ float **coords = NULL;
+ int orig_n = n;
+
+ /*double conj_tol=tolerance_cg; *//* tolerance of Conjugate Gradient */
+ CMajEnvVPSC *cMajEnvHor = NULL;
+ CMajEnvVPSC *cMajEnvVrt = NULL;
+ double y_0;
+ int length;
+ DistType diameter;
+ float *Dij = NULL;
+ float constant_term;
+ int count;
+ double degree;
+ int step;
+ float val;
+ double old_stress, new_stress = 0;
+ bool converged;
+ int len;
+ double nsizeScale = 0;
+ float maxEdgeLen = 0;
+ double max = 1;
initLayout(graph, n, dim, d_coords);
- if (n == 1) return 0;
+ if (n == 1)
+ return 0;
- for(i=0;i<n;i++) {
- for(j=1;j<graph[i].nedges;j++) {
- maxEdgeLen=MAX(graph[i].ewgts[j],maxEdgeLen);
- }
+ for (i = 0; i < n; i++) {
+ for (j = 1; j < graph[i].nedges; j++) {
+ maxEdgeLen = MAX(graph[i].ewgts[j], maxEdgeLen);
+ }
}
/****************************************************
** Compute the all-pairs-shortest-distances matrix **
****************************************************/
- if (maxi==0)
- return iterations;
+ if (maxi == 0)
+ return iterations;
if (model == MODEL_SUBSET) {
- /* weight graph to separate high-degree nodes */
- /* and perform slower Dijkstra-based computation */
- if (Verbose)
- fprintf(stderr, "Calculating subset model");
- Dij = compute_apsp_artifical_weights_packed(graph, n);
+ /* weight graph to separate high-degree nodes */
+ /* and perform slower Dijkstra-based computation */
+ if (Verbose)
+ fprintf(stderr, "Calculating subset model");
+ Dij = compute_apsp_artifical_weights_packed(graph, n);
} else if (model == MODEL_CIRCUIT) {
- Dij = circuitModel(graph, n);
- if (!Dij) {
- agerr(AGWARN,
- "graph is disconnected. Hence, the circuit model\n");
- agerr(AGPREV,
- "is undefined. Reverting to the shortest path model.\n");
- }
+ Dij = circuitModel(graph, n);
+ if (!Dij) {
+ agerr(AGWARN,
+ "graph is disconnected. Hence, the circuit model\n");
+ agerr(AGPREV,
+ "is undefined. Reverting to the shortest path model.\n");
+ }
}
if (!Dij) {
- if (Verbose)
- fprintf(stderr, "Calculating shortest paths\n ");
- Dij = compute_apsp_packed(graph, n);
+ if (Verbose)
+ fprintf(stderr, "Calculating shortest paths\n ");
+ Dij = compute_apsp_packed(graph, n);
}
- diameter=-1;
- length = n+n*(n-1)/2;
- for (i=0; i<length; i++) {
- if (Dij[i]>diameter) {
- diameter = (int)Dij[i];
- }
+ diameter = -1;
+ length = n + n * (n - 1) / 2;
+ for (i = 0; i < length; i++) {
+ if (Dij[i] > diameter) {
+ diameter = (int) Dij[i];
}
+ }
- /* for numerical stability, scale down layout */
- /* No Jiggling, might conflict with constraints */
- double max=1;
- for (i=0; i<dim; i++) {
- for (j=0; j<n; j++) {
- if (fabs(d_coords[i][j])>max) {
- max=fabs(d_coords[i][j]);
- }
- }
+ /* for numerical stability, scale down layout */
+ /* No Jiggling, might conflict with constraints */
+ for (i = 0; i < dim; i++) {
+ for (j = 0; j < n; j++) {
+ if (fabs(d_coords[i][j]) > max) {
+ max = fabs(d_coords[i][j]);
+ }
+ }
}
- for (i=0; i<dim; i++) {
- for (j=0; j<n; j++) {
- d_coords[i][j]*=10/max;
- }
+ for (i = 0; i < dim; i++) {
+ for (j = 0; j < n; j++) {
+ d_coords[i][j] *= 10 / max;
+ }
}
/**************************
** Layout initialization **
**************************/
- for (i=0; i<dim; i++) {
- orthog1(n, d_coords[i]);
- }
+ for (i = 0; i < dim; i++) {
+ orthog1(n, d_coords[i]);
+ }
- /* for the y-coords, don't center them, but translate them so y[0]=0 */
- y_0 = d_coords[1][0];
- for (i=0; i<n; i++) {
- d_coords[1][i] -= y_0;
- }
+ /* for the y-coords, don't center them, but translate them so y[0]=0 */
+ y_0 = d_coords[1][0];
+ for (i = 0; i < n; i++) {
+ d_coords[1][i] -= y_0;
+ }
/**************************
** Laplacian computation **
**************************/
-
+
lap2 = Dij;
- lap_length = n+n*(n-1)/2;
+ lap_length = n + n * (n - 1) / 2;
square_vec(lap_length, lap2);
/* compute off-diagonal entries */
invert_vec(lap_length, lap2);
-
- if(opt->clusters->nclusters>0) {
- int nn = n+opt->clusters->nclusters*2;
- int clap_length = nn+nn*(nn-1)/2;
- float *clap = N_GNEW(clap_length, float);
- int c0,c1;
- float v;
- c0=c1=0;
- for(i=0;i<nn;i++) {
- for(j=0;j<nn-i;j++) {
- if(i<n && j<n-i) {
- v=lap2[c0++];
- } else {
- //v=j==1?i%2:0;
- if(j==1&&i%2==1) {
- v=maxEdgeLen;
- v*=v;
- if(v>0.01) {
- v=1.0/v;
- }
- }
- else v=0;
- }
- clap[c1++]=v;
- }
- }
- free(lap2);
- lap2=clap;
- n=nn;
- lap_length=clap_length;
+
+ if (opt->clusters->nclusters > 0) {
+ int nn = n + opt->clusters->nclusters * 2;
+ int clap_length = nn + nn * (nn - 1) / 2;
+ float *clap = N_GNEW(clap_length, float);
+ int c0, c1;
+ float v;
+ c0 = c1 = 0;
+ for (i = 0; i < nn; i++) {
+ for (j = 0; j < nn - i; j++) {
+ if (i < n && j < n - i) {
+ v = lap2[c0++];
+ } else {
+ /* v=j==1?i%2:0; */
+ if (j == 1 && i % 2 == 1) {
+ v = maxEdgeLen;
+ v *= v;
+ if (v > 0.01) {
+ v = 1.0 / v;
+ }
+ } else
+ v = 0;
+ }
+ clap[c1++] = v;
+ }
+ }
+ free(lap2);
+ lap2 = clap;
+ n = nn;
+ lap_length = clap_length;
}
/* compute diagonal entries */
- count=0;
+ count = 0;
degrees = N_GNEW(n, double);
set_vector_val(n, 0, degrees);
- for (i=0; i<n-1; i++) {
- degree=0;
- count++; // skip main diag entry
- for (j=1; j<n-i; j++,count++) {
- val = lap2[count];
- degree+=val; degrees[i+j]-=val;
- }
- degrees[i]-=degree;
+ for (i = 0; i < n - 1; i++) {
+ degree = 0;
+ count++; /* skip main diag entry */
+ for (j = 1; j < n - i; j++, count++) {
+ val = lap2[count];
+ degree += val;
+ degrees[i + j] -= val;
+ }
+ degrees[i] -= degree;
}
- for (step=n,count=0,i=0; i<n; i++,count+=step,step--) {
- lap2[count]=(float)degrees[i];
+ for (step = n, count = 0, i = 0; i < n; i++, count += step, step--) {
+ lap2[count] = (float) degrees[i];
}
- coords = N_GNEW(dim, float*);
- f_storage = N_GNEW(dim*n, float);
- for (i=0; i<dim; i++) {
- coords[i] = f_storage+i*n;
- for (j=0; j<n; j++) {
- coords[i][j] = j<orig_n?(float)(d_coords[i][j]):0;
- }
+ coords = N_GNEW(dim, float *);
+ f_storage = N_GNEW(dim * n, float);
+ for (i = 0; i < dim; i++) {
+ coords[i] = f_storage + i * n;
+ for (j = 0; j < n; j++) {
+ coords[i][j] = j < orig_n ? (float) (d_coords[i][j]) : 0;
}
+ }
- /* compute constant term in stress sum
- * which is \sum_{i<j} w_{ij}d_{ij}^2
+ /* compute constant term in stress sum
+ * which is \sum_{i<j} w_{ij}d_{ij}^2
*/
- constant_term=(float)(n*(n-1)/2);
-
+ constant_term = (float) (n * (n - 1) / 2);
+
/*************************
** Layout optimization **
*************************/
-
- b = N_GNEW (dim, float*);
- b[0] = N_GNEW (dim*n, float);
- for (k=1; k<dim; k++) {
- b[k] = b[0]+k*n;
- }
- tmp_coords = N_GNEW(n, float);
- dist_accumulator = N_GNEW(n, float);
-
- old_stress=DBL_MAX; /* at least one iteration */
-
- start_time = clock();
-
- cMajEnvHor=initCMajVPSC(n,lap2,graph,opt,0);
- cMajEnvVrt=initCMajVPSC(n,lap2,graph,opt,opt->diredges);
-
- lap1 = N_GNEW(lap_length, float);
-
- for (converged=false,iterations=0; iterations<maxi && !converged; iterations++) {
-
- /* First, construct Laplacian of 1/(d_ij*|p_i-p_j|) */
- set_vector_val(n, 0, degrees);
- sqrt_vecf(lap_length, lap2, lap1);
- for (count=0,i=0; i<n-1; i++) {
- len = n-i-1;
- /* init 'dist_accumulator' with zeros */
- set_vector_valf(n, 0, dist_accumulator);
-
- /* put into 'dist_accumulator' all squared distances
- * between 'i' and 'i'+1,...,'n'-1
- */
- for (k=0; k<dim; k++) {
- set_vector_valf(len, coords[k][i], tmp_coords);
- vectors_mult_additionf(len, tmp_coords, -1, coords[k]+i+1);
- square_vec(len, tmp_coords);
- vectors_additionf(len, tmp_coords, dist_accumulator, dist_accumulator);
- }
-
- /* convert to 1/d_{ij} */
- invert_sqrt_vec(len, dist_accumulator);
- /* detect overflows */
- for (j=0; j<len; j++) {
- if (dist_accumulator[j]>=FLT_MAX || dist_accumulator[j]<0) {
- dist_accumulator[j]=0;
- }
- }
-
- count++; /* save place for the main diagonal entry */
- degree=0;
- for (j=0; j<len; j++,count++) {
- val=lap1[count]*=dist_accumulator[j];
- degree+=val; degrees[i+j+1]-=val;
- }
- degrees[i]-=degree;
- }
- for (step=n,count=0,i=0; i<n; i++,count+=step,step--) {
- lap1[count]=(float)degrees[i];
- }
+ b = N_GNEW(dim, float *);
+ b[0] = N_GNEW(dim * n, float);
+ for (k = 1; k < dim; k++) {
+ b[k] = b[0] + k * n;
+ }
- /* Now compute b[] (L^(X(t))*X(t)) */
- for (k=0; k<dim; k++) {
- /* b[k] := lap1*coords[k] */
- right_mult_with_vector_ff(lap1, n, coords[k], b[k]);
- }
-
- /* compute new stress
- * remember that the Laplacians are negated, so we subtract
- * instead of add and vice versa
- */
- new_stress=0;
- for (k=0; k<dim; k++) {
- new_stress+=vectors_inner_productf(n, coords[k], b[k]);
- }
- new_stress*=2;
- new_stress+=constant_term; // only after mult by 2
- for (k=0; k<dim; k++) {
- right_mult_with_vector_ff(lap2, n, coords[k], tmp_coords);
- new_stress-=vectors_inner_productf(n, coords[k], tmp_coords);
+ tmp_coords = N_GNEW(n, float);
+ dist_accumulator = N_GNEW(n, float);
+
+ old_stress = DBL_MAX; /* at least one iteration */
+
+ cMajEnvHor = initCMajVPSC(n, lap2, graph, opt, 0);
+ cMajEnvVrt = initCMajVPSC(n, lap2, graph, opt, opt->diredges);
+
+ lap1 = N_GNEW(lap_length, float);
+
+ for (converged = false, iterations = 0;
+ iterations < maxi && !converged; iterations++) {
+
+ /* First, construct Laplacian of 1/(d_ij*|p_i-p_j|) */
+ set_vector_val(n, 0, degrees);
+ sqrt_vecf(lap_length, lap2, lap1);
+ for (count = 0, i = 0; i < n - 1; i++) {
+ len = n - i - 1;
+ /* init 'dist_accumulator' with zeros */
+ set_vector_valf(n, 0, dist_accumulator);
+
+ /* put into 'dist_accumulator' all squared distances
+ * between 'i' and 'i'+1,...,'n'-1
+ */
+ for (k = 0; k < dim; k++) {
+ set_vector_valf(len, coords[k][i], tmp_coords);
+ vectors_mult_additionf(len, tmp_coords, -1,
+ coords[k] + i + 1);
+ square_vec(len, tmp_coords);
+ vectors_additionf(len, tmp_coords, dist_accumulator,
+ dist_accumulator);
+ }
+
+ /* convert to 1/d_{ij} */
+ invert_sqrt_vec(len, dist_accumulator);
+ /* detect overflows */
+ for (j = 0; j < len; j++) {
+ if (dist_accumulator[j] >= FLT_MAX
+ || dist_accumulator[j] < 0) {
+ dist_accumulator[j] = 0;
}
+ }
+
+ count++; /* save place for the main diagonal entry */
+ degree = 0;
+ for (j = 0; j < len; j++, count++) {
+ val = lap1[count] *= dist_accumulator[j];
+ degree += val;
+ degrees[i + j + 1] -= val;
+ }
+ degrees[i] -= degree;
+ }
+ for (step = n, count = 0, i = 0; i < n; i++, count += step, step--) {
+ lap1[count] = (float) degrees[i];
+ }
+
+ /* Now compute b[] (L^(X(t))*X(t)) */
+ for (k = 0; k < dim; k++) {
+ /* b[k] := lap1*coords[k] */
+ right_mult_with_vector_ff(lap1, n, coords[k], b[k]);
+ }
+
+ /* compute new stress
+ * remember that the Laplacians are negated, so we subtract
+ * instead of add and vice versa
+ */
+ new_stress = 0;
+ for (k = 0; k < dim; k++) {
+ new_stress += vectors_inner_productf(n, coords[k], b[k]);
+ }
+ new_stress *= 2;
+ new_stress += constant_term; /* only after mult by 2 */
+ for (k = 0; k < dim; k++) {
+ right_mult_with_vector_ff(lap2, n, coords[k], tmp_coords);
+ new_stress -= vectors_inner_productf(n, coords[k], tmp_coords);
+ }
#ifdef ALTERNATIVE_STRESS_CALC
- {
- double mat_stress=new_stress;
- double compute_stress(float ** coords, float * lap, int dim, int n);
- new_stress = compute_stress(coords, lap2, dim, n);
- if (fabs(mat_stress-new_stress)/min(mat_stress,new_stress)>0.001) {
- fprintf(stderr,"Diff stress vals: %lf %lf (iteration #%d)\n", mat_stress, new_stress,iterations);
- }
- }
-#endif
- /* check for convergence */
- if (Verbose && (iterations % 1 == 0)) {
- fprintf(stderr, "%.3f ", new_stress);
- if (iterations % 10 == 0)
- fprintf(stderr, "\n");
+ {
+ double mat_stress = new_stress;
+ double compute_stress(float **coords, float *lap, int dim,
+ int n);
+ new_stress = compute_stress(coords, lap2, dim, n);
+ if (fabs(mat_stress - new_stress) /
+ min(mat_stress, new_stress) > 0.001) {
+ fprintf(stderr,
+ "Diff stress vals: %lf %lf (iteration #%d)\n",
+ mat_stress, new_stress, iterations);
}
- converged = new_stress < old_stress && fabs(new_stress-old_stress)/fabs(old_stress+1e-10) < Epsilon;
- //converged = converged || (iterations>1 && new_stress>old_stress);
- /* in first iteration we allowed stress increase, which
- * might result ny imposing constraints
- */
- old_stress = new_stress;
-
- // in determining non-overlap constraints we gradually scale up the
- // size of nodes to avoid local minima
- if((iterations>=maxi-1||converged)&&opt->noverlap==1&&nsizeScale<0.999) {
- nsizeScale+=0.1;
- if(Verbose)
- fprintf(stderr,"nsizescale=%f,iterations=%d\n",nsizeScale,iterations);
- iterations=0;
- converged = false;
- }
-
-
- /* now we find the optimizer of trace(X'LX)+X'B by solving 'dim'
- * system of equations, thereby obtaining the new coordinates.
- * If we use the constraints (given by the var's: 'ordering',
- * 'levels' and 'num_levels'), we cannot optimize
- * trace(X'LX)+X'B by simply solving equations, but we have
- * to use a quadratic programming solver
- * note: 'lap2' is a packed symmetric matrix, that is its
- * upper-triangular part is arranged in a vector row-wise
- * also note: 'lap2' is really the negated laplacian (the
- * laplacian is -'lap2')
- */
-
- if(opt->noverlap==1 && nsizeScale > 0.001) {
- generateNonoverlapConstraints(cMajEnvHor,nsizeScale,coords,0,nsizeScale<0.5?false:true,opt);
- }
- if(cMajEnvHor->m > 0) {
+ }
+#endif
+ /* check for convergence */
+ if (Verbose && (iterations % 1 == 0)) {
+ fprintf(stderr, "%.3f ", new_stress);
+ if (iterations % 10 == 0)
+ fprintf(stderr, "\n");
+ }
+ converged = new_stress < old_stress
+ && fabs(new_stress - old_stress) / fabs(old_stress + 1e-10) <
+ Epsilon;
+ /*converged = converged || (iterations>1 && new_stress>old_stress); */
+ /* in first iteration we allowed stress increase, which
+ * might result ny imposing constraints
+ */
+ old_stress = new_stress;
+
+ /* in determining non-overlap constraints we gradually scale up the
+ * size of nodes to avoid local minima
+ */
+ if ((iterations >= maxi - 1 || converged) && opt->noverlap == 1
+ && nsizeScale < 0.999) {
+ nsizeScale += 0.1;
+ if (Verbose)
+ fprintf(stderr, "nsizescale=%f,iterations=%d\n",
+ nsizeScale, iterations);
+ iterations = 0;
+ converged = false;
+ }
+
+
+ /* now we find the optimizer of trace(X'LX)+X'B by solving 'dim'
+ * system of equations, thereby obtaining the new coordinates.
+ * If we use the constraints (given by the var's: 'ordering',
+ * 'levels' and 'num_levels'), we cannot optimize
+ * trace(X'LX)+X'B by simply solving equations, but we have
+ * to use a quadratic programming solver
+ * note: 'lap2' is a packed symmetric matrix, that is its
+ * upper-triangular part is arranged in a vector row-wise
+ * also note: 'lap2' is really the negated laplacian (the
+ * laplacian is -'lap2')
+ */
+
+ if (opt->noverlap == 1 && nsizeScale > 0.001) {
+ generateNonoverlapConstraints(cMajEnvHor, nsizeScale, coords,
+ 0,
+ nsizeScale < 0.5 ? false : true,
+ opt);
+ }
+ if (cMajEnvHor->m > 0) {
#ifdef MOSEK
- if(opt->mosek) {
- mosek_quad_solve_sep(cMajEnvHor->mosekEnv,n,b[0],coords[0]);
- } else
-#endif // MOSEK
- constrained_majorization_vpsc(cMajEnvHor,b[0],coords[0],localConstrMajorIterations);
- } else {
- // if there are no constraints then use conjugate gradient
- // optimisation which should be considerably faster
- conjugate_gradient_mkernel(lap2, coords[0], b[0], n, tolerance_cg, n);
- }
- if(opt->noverlap==1 && nsizeScale > 0.001) {
- generateNonoverlapConstraints(cMajEnvVrt,nsizeScale,coords,1,false,opt);
- }
- if(cMajEnvVrt->m > 0) {
+ if (opt->mosek) {
+ mosek_quad_solve_sep(cMajEnvHor->mosekEnv, n, b[0],
+ coords[0]);
+ } else
+#endif /* MOSEK */
+ constrained_majorization_vpsc(cMajEnvHor, b[0], coords[0],
+ localConstrMajorIterations);
+ } else {
+ /* if there are no constraints then use conjugate gradient
+ * optimisation which should be considerably faster
+ */
+ conjugate_gradient_mkernel(lap2, coords[0], b[0], n,
+ tolerance_cg, n);
+ }
+ if (opt->noverlap == 1 && nsizeScale > 0.001) {
+ generateNonoverlapConstraints(cMajEnvVrt, nsizeScale, coords,
+ 1, false, opt);
+ }
+ if (cMajEnvVrt->m > 0) {
#ifdef MOSEK
- if(opt->mosek) {
- mosek_quad_solve_sep(cMajEnvVrt->mosekEnv,n,b[1],coords[1]);
- } else
-#endif // MOSEK
- constrained_majorization_vpsc(cMajEnvVrt,b[1],coords[1],localConstrMajorIterations);
- } else {
- conjugate_gradient_mkernel(lap2,coords[1],b[1],n,tolerance_cg,n);
- }
+ if (opt->mosek) {
+ mosek_quad_solve_sep(cMajEnvVrt->mosekEnv, n, b[1],
+ coords[1]);
+ } else
+#endif /* MOSEK */
+ constrained_majorization_vpsc(cMajEnvVrt, b[1], coords[1],
+ localConstrMajorIterations);
+ } else {
+ conjugate_gradient_mkernel(lap2, coords[1], b[1], n,
+ tolerance_cg, n);
}
+ }
if (Verbose) {
- fprintf(stderr, "\nfinal e = %f %d iterations %.2f sec\n",
- new_stress, iterations, elapsed_sec());
+ fprintf(stderr, "\nfinal e = %f %d iterations %.2f sec\n",
+ new_stress, iterations, elapsed_sec());
}
- deleteCMajEnvVPSC(cMajEnvHor);
- deleteCMajEnvVPSC(cMajEnvVrt);
+ deleteCMajEnvVPSC(cMajEnvHor);
+ deleteCMajEnvVPSC(cMajEnvVrt);
- if (opt->noverlap==2) {
- fprintf(stderr,"Removing overlaps as post-process...\n");
- removeoverlaps(orig_n,coords,opt);
+ if (opt->noverlap == 2) {
+ fprintf(stderr, "Removing overlaps as post-process...\n");
+ removeoverlaps(orig_n, coords, opt);
}
-
- if (coords!=NULL) {
- for (i=0; i<dim; i++) {
- for (j=0; j<orig_n; j++) {
- d_coords[i][j] = coords[i][j];
- }
- }
- free (coords[0]); free (coords);
- }
-
- if (b) {
- free (b[0]); free (b);
+
+ if (coords != NULL) {
+ for (i = 0; i < dim; i++) {
+ for (j = 0; j < orig_n; j++) {
+ d_coords[i][j] = coords[i][j];
+ }
}
- free (tmp_coords);
- free (dist_accumulator);
- free (degrees);
- free (lap2);
- free (lap1);
+ free(coords[0]);
+ free(coords);
+ }
+
+ if (b) {
+ free(b[0]);
+ free(b);
+ }
+ free(tmp_coords);
+ free(dist_accumulator);
+ free(degrees);
+ free(lap2);
+ free(lap1);
return iterations;
}
-#endif /* DIGCOLA */
-
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
+#endif /* IPSEPCOLA */
float *edists; /* directed dist reflecting the direction of the edge */
};
- struct cluster_data {
- int nvars; /* total count of vars in clusters */
- int nclusters; /* number of clusters */
- int *clustersizes; /* number of vars in each cluster */
- int **clusters; /* list of var indices for constituents of each c */
- int ntoplevel; /* number of nodes not in any cluster */
- int *toplevel; /* array of nodes not in any cluster */
- boxf *bb; /* bounding box of each cluster */
- };
-
typedef int DistType; /* must be signed!! */
inline double max(double x, double y) {
#include <macros.h>
extern void *gmalloc(size_t);
-#define DIGCOLA 1
#ifdef USE_STYLES
typedef enum { regular, invisible } Style;
#endif
} vtx_data;
+#ifdef DIGCOLA
+#ifdef IPSEPCOLA
typedef struct cluster_data {
int nvars; /* total count of vars in clusters */
int nclusters; /* number of clusters */
int *toplevel; /* array of nodes not in any cluster */
boxf *bb; /* bounding box of each cluster */
} cluster_data;
-
+#endif
+#endif
typedef int DistType; /* must be signed!! */
extern void IMDS_given_dim(vtx_data*, int, double*, double*, double);
extern int stress_majorization_with_hierarchy(vtx_data*, int, int, double**,
int, int, int, int, double);
+#ifdef IPSEPCOLA
typedef struct ipsep_options {
int diredges; /* 1=generate directed edge constraints */
/* 2=generate directed hierarchy level constraints (DiG-CoLa) */
/* list of node indices for each cluster */
#ifdef MOSEK
int mosek; /* use Mosek as constraint optimization engine */
-#endif //MOSEK
+#endif /* MOSEK */
} ipsep_options;
-// stress majorization, for Constraint Layout
+ /* stress majorization, for Constraint Layout */
extern int stress_majorization_cola(vtx_data*, int, int, double**, int, int, int, ipsep_options*);
#endif
#endif
+#endif
#ifdef __cplusplus
}
+/* vim:set shiftwidth=4 ts=8: */
+
+/**
+ * Authors:
+ * Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * This version is released under the CPL (Common Public License) with
+ * the Graphviz distribution.
+ * A version is also available under the LGPL as part of the Adaptagrams
+ * project: http://sourceforge.net/projects/adaptagrams.
+ * If you make improvements or bug fixes to this code it would be much
+ * appreciated if you could also contribute those changes back to the
+ * Adaptagrams repository.
+ */
+
/*
* Interface to Mosek (www.mosek.com) quadratic programming solver for solving
* instances of the "Variable Placement with Separation Constraints" problem,
#include "mosek_quad_solve.h"
#include "quad_prog_vpsc.h"
-//#define DUMP_CONSTRAINTS
-//#define EQUAL_WIDTH_LEVELS
+/* #define DUMP_CONSTRAINTS */
+/* #define EQUAL_WIDTH_LEVELS */
static FILE *logfile;
static void MSKAPI printstr(void *handle, char str[])
{
- fprintf(logfile,"%s",str);
-} /* printstr */
+ fprintf(logfile, "%s", str);
+}
+
+#define INIT_sub_val(a,b) \
+ MSKidxt subi[2]; \
+ double vali[2]; \
+ subi[0] = a; \
+ subi[1] = b; \
+ vali[0] = 1.0; \
+ vali[1] = -1.0;
+
+#define INIT_sub_val3(a,b,c) \
+ MSKidxt subi[3]; \
+ double vali[3]; \
+ subi[0] = a; \
+ subi[1] = b; \
+ subi[2] = c; \
+ vali[0] = 1.0; \
+ vali[1] = -2.0; \
+ vali[2] = 1.0;
/**********************
lap: the upper RHS of the symmetric graph laplacian matrix which will be transformed
num_divisions: number of divisions between levels, ie number of levels - 1
separation: the minimum separation between nodes on different levels
***********************/
-MosekEnv* mosek_init_hier(
- float* lap, int n,int *ordering,int *level_indexes,int num_divisions, float separation)
+MosekEnv *mosek_init_hier(float *lap, int n, int *ordering,
+ int *level_indexes, int num_divisions,
+ float separation)
{
- int i,j,num_levels=num_divisions+1;
- MosekEnv *mskEnv = GNEW(MosekEnv);
- // vars for nodes (except x0) + dummy nodes between levels
- // x0 is fixed at 0, and therefore is not included in the opt problem
- // add 2 more vars for top and bottom constraints
- mskEnv->num_variables=n+num_divisions+1;
-
- logfile = fopen("quad_solve_log","w");
- DigColaLevel *levels=assign_digcola_levels(ordering,n,level_indexes,num_divisions);
+ int count = 0;
+ int i, j, num_levels = num_divisions + 1;
+ int num_constraints;
+ MosekEnv *mskEnv = GNEW(MosekEnv);
+ DigColaLevel *levels;
+ int nonzero_lapsize = (n * (n - 1)) / 2;
+ /* vars for nodes (except x0) + dummy nodes between levels
+ * x0 is fixed at 0, and therefore is not included in the opt problem
+ * add 2 more vars for top and bottom constraints
+ */
+ mskEnv->num_variables = n + num_divisions + 1;
+
+ logfile = fopen("quad_solve_log", "w");
+ levels = assign_digcola_levels(ordering, n, level_indexes, num_divisions);
#ifdef DUMP_CONSTRAINTS
- print_digcola_levels(logfile,levels,num_levels);
+ print_digcola_levels(logfile, levels, num_levels);
#endif
- /* nonlinear coefficients matrix of objective function */
- //int lapsize=mskEnv->num_variables+(mskEnv->num_variables*(mskEnv->num_variables-1))/2;
- int nonzero_lapsize=(n*(n-1))/2;
- mskEnv->qval = N_GNEW(nonzero_lapsize,double);
- mskEnv->qsubi = N_GNEW(nonzero_lapsize,int);
- mskEnv->qsubj = N_GNEW(nonzero_lapsize,int);
-
- /* solution vector */
- mskEnv->xx = N_GNEW(mskEnv->num_variables,double);
-
- /* constraint matrix */
- separation/=2.0; // separation between each node and it's adjacent constraint
- int num_constraints = get_num_digcola_constraints(levels, num_levels)+num_divisions+1;
- // constraints of the form x_i - x_j >= sep so 2 non-zero entries per constraint in LHS matrix
- // except x_0 (fixed at 0) constraints which have 1 nz val each.
+ /* nonlinear coefficients matrix of objective function */
+ /* int lapsize=mskEnv->num_variables+(mskEnv->num_variables*(mskEnv->num_variables-1))/2; */
+ mskEnv->qval = N_GNEW(nonzero_lapsize, double);
+ mskEnv->qsubi = N_GNEW(nonzero_lapsize, int);
+ mskEnv->qsubj = N_GNEW(nonzero_lapsize, int);
+
+ /* solution vector */
+ mskEnv->xx = N_GNEW(mskEnv->num_variables, double);
+
+ /* constraint matrix */
+ separation /= 2.0; /* separation between each node and it's adjacent constraint */
+ num_constraints = get_num_digcola_constraints(levels,
+ num_levels) + num_divisions + 1;
+ /* constraints of the form x_i - x_j >= sep so 2 non-zero entries per constraint in LHS matrix
+ * except x_0 (fixed at 0) constraints which have 1 nz val each.
+ */
#ifdef EQUAL_WIDTH_LEVELS
- num_constraints+=num_divisions;
+ num_constraints += num_divisions;
#endif
- // pointer to beginning of nonzero sequence in a column
-
- int count=0;
- for(i=0;i<n-1;i++) {
- for(j=i;j<n-1;j++) {
- mskEnv->qval[count]=-2*lap[count+n];
- assert(mskEnv->qval[count]!=0);
- mskEnv->qsubi[count]=j;
- mskEnv->qsubj[count]=i;
- count++;
- }
+ /* pointer to beginning of nonzero sequence in a column */
+
+ for (i = 0; i < n - 1; i++) {
+ for (j = i; j < n - 1; j++) {
+ mskEnv->qval[count] = -2 * lap[count + n];
+ assert(mskEnv->qval[count] != 0);
+ mskEnv->qsubi[count] = j;
+ mskEnv->qsubj[count] = i;
+ count++;
}
+ }
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"Q=[");
- int lapcntr=n;
- for(i=0;i<mskEnv->num_variables;i++) {
- if(i!=0) fprintf(logfile,";");
- for(j=0;j<mskEnv->num_variables;j++) {
- if(j<i||i>=n-1||j>=n-1) {
- fprintf(logfile,"0 ");
- } else {
- fprintf(logfile,"%f ",-2*lap[lapcntr++]);
- }
- }
+ fprintf(logfile, "Q=[");
+ int lapcntr = n;
+ for (i = 0; i < mskEnv->num_variables; i++) {
+ if (i != 0)
+ fprintf(logfile, ";");
+ for (j = 0; j < mskEnv->num_variables; j++) {
+ if (j < i || i >= n - 1 || j >= n - 1) {
+ fprintf(logfile, "0 ");
+ } else {
+ fprintf(logfile, "%f ", -2 * lap[lapcntr++]);
+ }
}
- fprintf(logfile,"]\nQ=Q-diag(diag(Q))+Q'\n");
+ }
+ fprintf(logfile, "]\nQ=Q-diag(diag(Q))+Q'\n");
#endif
- fprintf(logfile,"\n");
- /* Make the mosek environment. */
- mskEnv->r = MSK_makeenv(&mskEnv->env,NULL,NULL,NULL,NULL);
-
- /* Check whether the return code is ok. */
- if ( mskEnv->r==MSK_RES_OK )
- {
- /* Directs the log stream to the user
- specified procedure 'printstr'. */
- MSK_linkfunctoenvstream(mskEnv->env,MSK_STREAM_LOG,NULL,printstr);
- }
+ fprintf(logfile, "\n");
+ /* Make the mosek environment. */
+ mskEnv->r = MSK_makeenv(&mskEnv->env, NULL, NULL, NULL, NULL);
+
+ /* Check whether the return code is ok. */
+ if (mskEnv->r == MSK_RES_OK) {
+ /* Directs the log stream to the user
+ * specified procedure 'printstr'.
+ */
+ MSK_linkfunctoenvstream(mskEnv->env, MSK_STREAM_LOG, NULL,
+ printstr);
+ }
+
+ /* Initialize the environment. */
+ mskEnv->r = MSK_initenv(mskEnv->env);
+ if (mskEnv->r == MSK_RES_OK) {
+ /* Make the optimization task. */
+ mskEnv->r =
+ MSK_maketask(mskEnv->env, num_constraints,
+ mskEnv->num_variables, &mskEnv->task);
+
+ if (mskEnv->r == MSK_RES_OK) {
+ int c_ind = 0;
+ int c_var = n - 1;
+ mskEnv->r =
+ MSK_linkfunctotaskstream(mskEnv->task, MSK_STREAM_LOG,
+ NULL, printstr);
+ /* Resize the task. */
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r = MSK_resizetask(mskEnv->task, num_constraints, mskEnv->num_variables, 0, /* no cones!! */
+ /* each constraint applies to 2 vars */
+ 2 * num_constraints +
+ num_divisions, nonzero_lapsize);
- /* Initialize the environment. */
- mskEnv->r = MSK_initenv(mskEnv->env);
- if (mskEnv->r==MSK_RES_OK)
- {
- /* Make the optimization task. */
- mskEnv->r = MSK_maketask(mskEnv->env,num_constraints,mskEnv->num_variables,&mskEnv->task);
-
- if ( mskEnv->r==MSK_RES_OK )
- {
- mskEnv->r = MSK_linkfunctotaskstream(mskEnv->task,MSK_STREAM_LOG,NULL,printstr);
- /* Resize the task. */
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_resizetask(mskEnv->task,num_constraints,mskEnv->num_variables,
- 0, // no cones!!
- // each constraint applies to 2 vars
- 2*num_constraints+num_divisions,
- nonzero_lapsize);
-
- /* Append the constraints. */
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_append(mskEnv->task,1,num_constraints);
-
- /* Append the variables. */
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_append(mskEnv->task,0,mskEnv->num_variables);
- /* Put variable bounds. */
- for(j=0; j<mskEnv->num_variables && mskEnv->r==MSK_RES_OK; ++j)
- mskEnv->r = MSK_putbound(mskEnv->task,0,j,MSK_BK_RA,-MSK_INFINITY,MSK_INFINITY);
- int c_ind=0;
- int c_var=n-1;
- for(j=0;j<levels[0].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
- int node = levels[0].nodes[j]-1;
- if(node>=0){
- MSKidxt subi[]={c_var,node};
- double vali[]={1.0,-1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
- } else {
- // constraint for y0 (fixed at 0)
- mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,1.0);
- }
- mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
- c_ind++;
- }
- for(i=0;i<num_divisions && mskEnv->r==MSK_RES_OK;i++) {
- c_var=n+i;
- for(j=0;j<levels[i].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
- // create separation constraint a>=b+separation
- int node = levels[i].nodes[j]-1;
- if(node>=0){ // no constraint for fixed node
- MSKidxt subi[]={node,c_var};
- double vali[]={1.0,-1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
- } else {
- // constraint for y0 (fixed at 0)
- mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,-1.0);
- }
- mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
- c_ind++;
- }
- for(j=0;j<levels[i+1].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
- int node = levels[i+1].nodes[j]-1;
- if(node>=0){
- MSKidxt subi[]={c_var,node};
- double vali[]={1.0,-1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
- } else {
- // constraint for y0 (fixed at 0)
- mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,1.0);
- }
- mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
- c_ind++;
- }
- }
- c_var=n+i;
- for(j=0;j<levels[i].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
- // create separation constraint a>=b+separation
- int node = levels[i].nodes[j]-1;
- if(node>=0){ // no constraint for fixed node
- MSKidxt subi[]={node,c_var};
- double vali[]={1.0,-1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
- } else {
- // constraint for y0 (fixed at 0)
- mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,-1.0);
- }
- mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
- c_ind++;
- }
- // create constraints preserving the order of dummy vars
- for(i=0;i<num_divisions+1 && mskEnv->r==MSK_RES_OK;i++) {
- int c_var=n-1+i, c_var2=c_var+1;
- MSKidxt subi[]={c_var,c_var2};
- double vali[]={1.0,-1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
- mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,0,MSK_INFINITY);
- c_ind++;
- }
+ /* Append the constraints. */
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r = MSK_append(mskEnv->task, 1, num_constraints);
+
+ /* Append the variables. */
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r =
+ MSK_append(mskEnv->task, 0, mskEnv->num_variables);
+ /* Put variable bounds. */
+ for (j = 0;
+ j < mskEnv->num_variables && mskEnv->r == MSK_RES_OK; ++j)
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 0, j, MSK_BK_RA,
+ -MSK_INFINITY, MSK_INFINITY);
+ for (j = 0; j < levels[0].num_nodes && mskEnv->r == MSK_RES_OK;
+ j++) {
+ int node = levels[0].nodes[j] - 1;
+ if (node >= 0) {
+ INIT_sub_val(c_var,node);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, c_ind, 2, subi, vali);
+ } else {
+ /* constraint for y0 (fixed at 0) */
+ mskEnv->r =
+ MSK_putaij(mskEnv->task, c_ind, c_var, 1.0);
+ }
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
+ separation, MSK_INFINITY);
+ c_ind++;
+ }
+ for (i = 0; i < num_divisions && mskEnv->r == MSK_RES_OK; i++) {
+ c_var = n + i;
+ for (j = 0;
+ j < levels[i].num_nodes && mskEnv->r == MSK_RES_OK;
+ j++) {
+ /* create separation constraint a>=b+separation */
+ int node = levels[i].nodes[j] - 1;
+ if (node >= 0) { /* no constraint for fixed node */
+ INIT_sub_val(node,c_var);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, c_ind, 2, subi,
+ vali);
+ } else {
+ /* constraint for y0 (fixed at 0) */
+ mskEnv->r =
+ MSK_putaij(mskEnv->task, c_ind, c_var, -1.0);
+ }
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
+ separation, MSK_INFINITY);
+ c_ind++;
+ }
+ for (j = 0;
+ j < levels[i + 1].num_nodes
+ && mskEnv->r == MSK_RES_OK; j++) {
+ int node = levels[i + 1].nodes[j] - 1;
+ if (node >= 0) {
+ INIT_sub_val(c_var,node);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, c_ind, 2, subi,
+ vali);
+ } else {
+ /* constraint for y0 (fixed at 0) */
+ mskEnv->r =
+ MSK_putaij(mskEnv->task, c_ind, c_var, 1.0);
+ }
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
+ separation, MSK_INFINITY);
+ c_ind++;
+ }
+ }
+ c_var = n + i;
+ for (j = 0; j < levels[i].num_nodes && mskEnv->r == MSK_RES_OK;
+ j++) {
+ /* create separation constraint a>=b+separation */
+ int node = levels[i].nodes[j] - 1;
+ if (node >= 0) { /* no constraint for fixed node */
+ INIT_sub_val(node,c_var);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, c_ind, 2, subi, vali);
+ } else {
+ /* constraint for y0 (fixed at 0) */
+ mskEnv->r =
+ MSK_putaij(mskEnv->task, c_ind, c_var, -1.0);
+ }
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
+ separation, MSK_INFINITY);
+ c_ind++;
+ }
+ /* create constraints preserving the order of dummy vars */
+ for (i = 0; i < num_divisions + 1 && mskEnv->r == MSK_RES_OK;
+ i++) {
+ int c_var = n - 1 + i, c_var2 = c_var + 1;
+ INIT_sub_val(c_var,c_var2);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, c_ind, 2, subi, vali);
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO, 0,
+ MSK_INFINITY);
+ c_ind++;
+ }
#ifdef EQUAL_WIDTH_LEVELS
- for(i=1;i<num_divisions+1 && mskEnv->r==MSK_RES_OK;i++) {
- int c_var=n-1+i, c_var_lo=c_var-1, c_var_hi=c_var+1;
- MSKidxt subi[]={c_var_lo,c_var,c_var_hi};
- double vali[]={1.0,-2.0,1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,3,subi,vali);
- mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_FX,0,0);
- c_ind++;
- }
+ for (i = 1; i < num_divisions + 1 && mskEnv->r == MSK_RES_OK;
+ i++) {
+ int c_var = n - 1 + i, c_var_lo = c_var - 1, c_var_hi =
+ c_var + 1;
+ INIT_sub_val3(c_var_lo, c_var, c_var_h);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, c_ind, 3, subi, vali);
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_FX, 0, 0);
+ c_ind++;
+ }
#endif
- assert(c_ind==num_constraints);
+ assert(c_ind == num_constraints);
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"A=[");
- for(i=0;i<num_constraints;i++) {
- if(i!=0) fprintf(logfile,";");
- for(j=0;j<mskEnv->num_variables;j++) {
- double aij;
- MSK_getaij(mskEnv->task,i,j,&aij);
- fprintf(logfile,"%f ",aij);
- }
- }
- fprintf(logfile,"]\n");
- fprintf(logfile,"b=[");
- for(i=0;i<num_constraints;i++) {
- fprintf(logfile,"%f ",separation);
- }
- fprintf(logfile,"]\n");
-#endif
- if ( mskEnv->r==MSK_RES_OK )
- {
- /*
- * The lower triangular part of the Q
- * matrix in the objective is specified.
- */
- mskEnv->r = MSK_putqobj(mskEnv->task,nonzero_lapsize,mskEnv->qsubi,mskEnv->qsubj,mskEnv->qval);
- }
+ fprintf(logfile, "A=[");
+ for (i = 0; i < num_constraints; i++) {
+ if (i != 0)
+ fprintf(logfile, ";");
+ for (j = 0; j < mskEnv->num_variables; j++) {
+ double aij;
+ MSK_getaij(mskEnv->task, i, j, &aij);
+ fprintf(logfile, "%f ", aij);
}
+ }
+ fprintf(logfile, "]\n");
+ fprintf(logfile, "b=[");
+ for (i = 0; i < num_constraints; i++) {
+ fprintf(logfile, "%f ", separation);
+ }
+ fprintf(logfile, "]\n");
+#endif
+ if (mskEnv->r == MSK_RES_OK) {
+ /*
+ * The lower triangular part of the Q
+ * matrix in the objective is specified.
+ */
+ mskEnv->r =
+ MSK_putqobj(mskEnv->task, nonzero_lapsize,
+ mskEnv->qsubi, mskEnv->qsubj,
+ mskEnv->qval);
+ }
}
- delete_digcola_levels(levels,num_levels);
- return mskEnv;
+ }
+ delete_digcola_levels(levels, num_levels);
+ return mskEnv;
}
+
/*
b: coefficients of linear part of optimisation function
n: number of nodes
hierarchy_boundaries: y coord of boundaries between levels
(ie, solution values for the dummy variables used in constraints)
*/
-void mosek_quad_solve_hier(MosekEnv *mskEnv, float *b,int n,float* coords, float *hierarchy_boundaries) {
- int i,j;
- for(i=1;i<n && mskEnv->r==MSK_RES_OK;i++) {
- mskEnv->r = MSK_putcj(mskEnv->task,i-1,-2*b[i]);
- }
+void mosek_quad_solve_hier(MosekEnv * mskEnv, float *b, int n,
+ float *coords, float *hierarchy_boundaries)
+{
+ int i, j;
+ for (i = 1; i < n && mskEnv->r == MSK_RES_OK; i++) {
+ mskEnv->r = MSK_putcj(mskEnv->task, i - 1, -2 * b[i]);
+ }
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"x0=[");
- for(j=0;j<mskEnv->num_variables;j++) {
- fprintf(logfile,"%f ",j<n?b[j]:0);
- }
- fprintf(logfile,"]\n");
- fprintf(logfile,"f=[");
- double *c=N_GNEW(mskEnv->num_variables,double);
- MSK_getc(mskEnv->task,c);
- for(j=0;j<mskEnv->num_variables;j++) {
- fprintf(logfile,"%f ",c[j]);
- }
- free(c);
- fprintf(logfile,"]\n");
+ fprintf(logfile, "x0=[");
+ for (j = 0; j < mskEnv->num_variables; j++) {
+ fprintf(logfile, "%f ", j < n ? b[j] : 0);
+ }
+ fprintf(logfile, "]\n");
+ fprintf(logfile, "f=[");
+ double *c = N_GNEW(mskEnv->num_variables, double);
+ MSK_getc(mskEnv->task, c);
+ for (j = 0; j < mskEnv->num_variables; j++) {
+ fprintf(logfile, "%f ", c[j]);
+ }
+ free(c);
+ fprintf(logfile, "]\n");
#endif
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_optimize(mskEnv->task);
-
- if ( mskEnv->r==MSK_RES_OK )
- {
- MSK_getsolutionslice(mskEnv->task,
- MSK_SOL_ITR,
- MSK_SOL_ITEM_XX,
- 0,
- mskEnv->num_variables,
- mskEnv->xx);
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r = MSK_optimize(mskEnv->task);
+
+ if (mskEnv->r == MSK_RES_OK) {
+ MSK_getsolutionslice(mskEnv->task,
+ MSK_SOL_ITR,
+ MSK_SOL_ITEM_XX,
+ 0, mskEnv->num_variables, mskEnv->xx);
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"Primal solution\n");
+ fprintf(logfile, "Primal solution\n");
#endif
- coords[0]=0;
- for(j=0; j<mskEnv->num_variables; ++j) {
+ coords[0] = 0;
+ for (j = 0; j < mskEnv->num_variables; ++j) {
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"x[%d]: %.2f\n",j,mskEnv->xx[j]);
+ fprintf(logfile, "x[%d]: %.2f\n", j, mskEnv->xx[j]);
#endif
- if(j<n-1) {
- coords[j+1]=-mskEnv->xx[j];
- } else if(j>=n&&j<mskEnv->num_variables-1) {
- hierarchy_boundaries[j-n]=-mskEnv->xx[j];
- }
- }
+ if (j < n - 1) {
+ coords[j + 1] = -mskEnv->xx[j];
+ } else if (j >= n && j < mskEnv->num_variables - 1) {
+ hierarchy_boundaries[j - n] = -mskEnv->xx[j];
+ }
}
- fprintf(logfile,"Return code: %d\n",mskEnv->r);
+ }
+ fprintf(logfile, "Return code: %d\n", mskEnv->r);
}
+
/**********************
lap: the upper RHS of the symmetric graph laplacian matrix which will be transformed
to the hessian of the non-linear part of the optimisation function
has dimensions num_variables, dummy vars do not have entries in lap
cs: array of pointers to separation constraints
***********************/
-MosekEnv* mosek_init_sep(
- float* lap, int num_variables, int num_dummy_vars, Constraint **cs, int num_constraints)
+MosekEnv *mosek_init_sep(float *lap, int num_variables, int num_dummy_vars,
+ Constraint ** cs, int num_constraints)
{
- int i,j;
- MosekEnv *mskEnv = GNEW(MosekEnv);
- // fix var 0
- mskEnv->num_variables=num_variables+num_dummy_vars-1;
-
- fprintf(stderr,"MOSEK!\n");
- logfile = fopen("quad_solve_log","w");
-
- /* nonlinear coefficients matrix of objective function */
- int nonzero_lapsize=num_variables*(num_variables-1)/2;
- mskEnv->qval = N_GNEW(nonzero_lapsize,double);
- mskEnv->qsubi = N_GNEW(nonzero_lapsize,int);
- mskEnv->qsubj = N_GNEW(nonzero_lapsize,int);
-
- /* solution vector */
- mskEnv->xx = N_GNEW(mskEnv->num_variables,double);
-
- // pointer to beginning of nonzero sequence in a column
-
- int count=0;
- for(i=0;i<num_variables-1;i++) {
- for(j=i;j<num_variables-1;j++) {
- mskEnv->qval[count]=-2*lap[count+num_variables];
- //assert(mskEnv->qval[count]!=0);
- mskEnv->qsubi[count]=j;
- mskEnv->qsubj[count]=i;
- count++;
- }
+ int i, j;
+ MosekEnv *mskEnv = GNEW(MosekEnv);
+ int count = 0;
+ int nonzero_lapsize = num_variables * (num_variables - 1) / 2;
+ /* fix var 0 */
+ mskEnv->num_variables = num_variables + num_dummy_vars - 1;
+
+ fprintf(stderr, "MOSEK!\n");
+ logfile = fopen("quad_solve_log", "w");
+
+ /* nonlinear coefficients matrix of objective function */
+ mskEnv->qval = N_GNEW(nonzero_lapsize, double);
+ mskEnv->qsubi = N_GNEW(nonzero_lapsize, int);
+ mskEnv->qsubj = N_GNEW(nonzero_lapsize, int);
+
+ /* solution vector */
+ mskEnv->xx = N_GNEW(mskEnv->num_variables, double);
+
+ /* pointer to beginning of nonzero sequence in a column */
+
+ for (i = 0; i < num_variables - 1; i++) {
+ for (j = i; j < num_variables - 1; j++) {
+ mskEnv->qval[count] = -2 * lap[count + num_variables];
+ /* assert(mskEnv->qval[count]!=0); */
+ mskEnv->qsubi[count] = j;
+ mskEnv->qsubj[count] = i;
+ count++;
}
+ }
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"Q=[");
- count=0;
- for(i=0;i<num_variables-1;i++) {
- if(i!=0) fprintf(logfile,";");
- for(j=0;j<num_variables-1;j++) {
- if(j<i) {
- fprintf(logfile,"0 ");
- } else {
- fprintf(logfile,"%f ",-2*lap[num_variables+count++]);
- }
- }
+ fprintf(logfile, "Q=[");
+ count = 0;
+ for (i = 0; i < num_variables - 1; i++) {
+ if (i != 0)
+ fprintf(logfile, ";");
+ for (j = 0; j < num_variables - 1; j++) {
+ if (j < i) {
+ fprintf(logfile, "0 ");
+ } else {
+ fprintf(logfile, "%f ", -2 * lap[num_variables + count++]);
+ }
}
- fprintf(logfile,"]\nQ=Q-diag(diag(Q))+Q'\n");
+ }
+ fprintf(logfile, "]\nQ=Q-diag(diag(Q))+Q'\n");
#endif
- /* Make the mosek environment. */
- mskEnv->r = MSK_makeenv(&mskEnv->env,NULL,NULL,NULL,NULL);
-
- /* Check whether the return code is ok. */
- if ( mskEnv->r==MSK_RES_OK )
- {
- /* Directs the log stream to the user
- specified procedure 'printstr'. */
- MSK_linkfunctoenvstream(mskEnv->env,MSK_STREAM_LOG,NULL,printstr);
- }
+ /* Make the mosek environment. */
+ mskEnv->r = MSK_makeenv(&mskEnv->env, NULL, NULL, NULL, NULL);
+
+ /* Check whether the return code is ok. */
+ if (mskEnv->r == MSK_RES_OK) {
+ /* Directs the log stream to the user
+ specified procedure 'printstr'. */
+ MSK_linkfunctoenvstream(mskEnv->env, MSK_STREAM_LOG, NULL,
+ printstr);
+ }
+
+ /* Initialize the environment. */
+ mskEnv->r = MSK_initenv(mskEnv->env);
+ if (mskEnv->r == MSK_RES_OK) {
+ /* Make the optimization task. */
+ mskEnv->r =
+ MSK_maketask(mskEnv->env, num_constraints,
+ mskEnv->num_variables, &mskEnv->task);
+
+ if (mskEnv->r == MSK_RES_OK) {
+ mskEnv->r =
+ MSK_linkfunctotaskstream(mskEnv->task, MSK_STREAM_LOG,
+ NULL, printstr);
+ /* Resize the task. */
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r = MSK_resizetask(mskEnv->task, num_constraints, mskEnv->num_variables, 0, /* no cones!! */
+ /* number of non-zero constraint matrix entries:
+ * each constraint applies to 2 vars
+ */
+ 2 * num_constraints,
+ nonzero_lapsize);
- /* Initialize the environment. */
- mskEnv->r = MSK_initenv(mskEnv->env);
- if (mskEnv->r==MSK_RES_OK)
- {
- /* Make the optimization task. */
- mskEnv->r = MSK_maketask(mskEnv->env,num_constraints,mskEnv->num_variables,&mskEnv->task);
-
- if ( mskEnv->r==MSK_RES_OK )
- {
- mskEnv->r = MSK_linkfunctotaskstream(mskEnv->task,MSK_STREAM_LOG,NULL,printstr);
- /* Resize the task. */
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_resizetask(mskEnv->task,num_constraints,mskEnv->num_variables,
- 0, // no cones!!
- // number of non-zero constraint matrix entries:
- // each constraint applies to 2 vars
- 2*num_constraints,
- nonzero_lapsize);
-
- /* Append the constraints. */
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_append(mskEnv->task,1,num_constraints);
-
- /* Append the variables. */
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_append(mskEnv->task,0,mskEnv->num_variables);
- /* Put variable bounds. */
- for(j=0; j<mskEnv->num_variables && mskEnv->r==MSK_RES_OK; j++)
- mskEnv->r = MSK_putbound(mskEnv->task,0,j,MSK_BK_RA,-MSK_INFINITY,MSK_INFINITY);
- for(i=0;i<num_constraints;i++) {
- int u = getLeftVarID(cs[i])-1;
- int v = getRightVarID(cs[i])-1;
- double separation = getSeparation(cs[i]);
- if(u<0) {
- mskEnv->r = MSK_putbound(mskEnv->task,0,v,MSK_BK_RA,-MSK_INFINITY,-separation);
- assert ( mskEnv->r==MSK_RES_OK );
- } else if(v<0) {
- mskEnv->r = MSK_putbound(mskEnv->task,0,u,MSK_BK_RA,separation,MSK_INFINITY);
- assert ( mskEnv->r==MSK_RES_OK );
- } else {
- //fprintf(stderr,"u=%d,v=%d,sep=%f\n",u,v,separation);
- MSKidxt subi[]={u,v};
- double vali[]={1.0,-1.0};
- mskEnv->r=MSK_putavec(mskEnv->task,1,i,2,subi,vali);
- assert ( mskEnv->r==MSK_RES_OK );
- mskEnv->r=MSK_putbound(mskEnv->task,1,i,MSK_BK_LO,separation,MSK_INFINITY);
- assert ( mskEnv->r==MSK_RES_OK );
- }
- }
- if ( mskEnv->r==MSK_RES_OK )
- {
- /*
- * The lower triangular part of the Q
- * matrix in the objective is specified.
- */
- mskEnv->r = MSK_putqobj(mskEnv->task,nonzero_lapsize,mskEnv->qsubi,mskEnv->qsubj,mskEnv->qval);
- assert ( mskEnv->r==MSK_RES_OK );
- }
+ /* Append the constraints. */
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r = MSK_append(mskEnv->task, 1, num_constraints);
+
+ /* Append the variables. */
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r =
+ MSK_append(mskEnv->task, 0, mskEnv->num_variables);
+ /* Put variable bounds. */
+ for (j = 0;
+ j < mskEnv->num_variables && mskEnv->r == MSK_RES_OK; j++)
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 0, j, MSK_BK_RA,
+ -MSK_INFINITY, MSK_INFINITY);
+ for (i = 0; i < num_constraints; i++) {
+ int u = getLeftVarID(cs[i]) - 1;
+ int v = getRightVarID(cs[i]) - 1;
+ double separation = getSeparation(cs[i]);
+ if (u < 0) {
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 0, v, MSK_BK_RA,
+ -MSK_INFINITY, -separation);
+ assert(mskEnv->r == MSK_RES_OK);
+ } else if (v < 0) {
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 0, u, MSK_BK_RA,
+ separation, MSK_INFINITY);
+ assert(mskEnv->r == MSK_RES_OK);
+ } else {
+ /* fprintf(stderr,"u=%d,v=%d,sep=%f\n",u,v,separation); */
+ INIT_sub_val(u,v);
+ mskEnv->r =
+ MSK_putavec(mskEnv->task, 1, i, 2, subi, vali);
+ assert(mskEnv->r == MSK_RES_OK);
+ mskEnv->r =
+ MSK_putbound(mskEnv->task, 1, i, MSK_BK_LO,
+ separation, MSK_INFINITY);
+ assert(mskEnv->r == MSK_RES_OK);
}
+ }
+ if (mskEnv->r == MSK_RES_OK) {
+ /*
+ * The lower triangular part of the Q
+ * matrix in the objective is specified.
+ */
+ mskEnv->r =
+ MSK_putqobj(mskEnv->task, nonzero_lapsize,
+ mskEnv->qsubi, mskEnv->qsubj,
+ mskEnv->qval);
+ assert(mskEnv->r == MSK_RES_OK);
+ }
}
- return mskEnv;
+ }
+ return mskEnv;
}
+
/*
n: size of b and coords, may be smaller than mskEnv->num_variables if we
have dummy vars
b: coefficients of linear part of optimisation function
coords: optimal y* vector, coord[i] is coordinate of node[i]
*/
-void mosek_quad_solve_sep(MosekEnv *mskEnv,int n, float *b,float* coords) {
- int i,j;
- assert(n<=mskEnv->num_variables+1);
- for(i=0;i<n-1 && mskEnv->r==MSK_RES_OK;i++) {
- mskEnv->r = MSK_putcj(mskEnv->task,i,-2*b[i+1]);
- }
- if ( mskEnv->r==MSK_RES_OK )
- mskEnv->r = MSK_optimize(mskEnv->task);
-
- if ( mskEnv->r==MSK_RES_OK )
- {
- MSK_getsolutionslice(mskEnv->task,
- MSK_SOL_ITR,
- MSK_SOL_ITEM_XX,
- 0,
- mskEnv->num_variables,
- mskEnv->xx);
+void mosek_quad_solve_sep(MosekEnv * mskEnv, int n, float *b,
+ float *coords)
+{
+ int i, j;
+ assert(n <= mskEnv->num_variables + 1);
+ for (i = 0; i < n - 1 && mskEnv->r == MSK_RES_OK; i++) {
+ mskEnv->r = MSK_putcj(mskEnv->task, i, -2 * b[i + 1]);
+ }
+ if (mskEnv->r == MSK_RES_OK)
+ mskEnv->r = MSK_optimize(mskEnv->task);
+
+ if (mskEnv->r == MSK_RES_OK) {
+ MSK_getsolutionslice(mskEnv->task,
+ MSK_SOL_ITR,
+ MSK_SOL_ITEM_XX,
+ 0, mskEnv->num_variables, mskEnv->xx);
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"Primal solution\n");
+ fprintf(logfile, "Primal solution\n");
#endif
- coords[0]=0;
- for(j=1; j<=n; j++) {
+ coords[0] = 0;
+ for (j = 1; j <= n; j++) {
#ifdef DUMP_CONSTRAINTS
- fprintf(logfile,"x[%d]: %.2f\n",j,mskEnv->xx[j-1]);
+ fprintf(logfile, "x[%d]: %.2f\n", j, mskEnv->xx[j - 1]);
#endif
- coords[j]=-mskEnv->xx[j-1];
- }
+ coords[j] = -mskEnv->xx[j - 1];
}
- fprintf(logfile,"Return code: %d\n",mskEnv->r);
+ }
+ fprintf(logfile, "Return code: %d\n", mskEnv->r);
}
+
/*
please call to clean up
*/
-void mosek_delete(MosekEnv *mskEnv) {
- MSK_deletetask(&mskEnv->task);
- MSK_deleteenv(&mskEnv->env);
+void mosek_delete(MosekEnv * mskEnv)
+{
+ MSK_deletetask(&mskEnv->task);
+ MSK_deleteenv(&mskEnv->env);
- if(logfile) {
- fclose(logfile);
- logfile=NULL;
- }
- free(mskEnv->qval);
- free(mskEnv->qsubi);
- free(mskEnv->qsubj);
- free(mskEnv->xx);
- free(mskEnv);
+ if (logfile) {
+ fclose(logfile);
+ logfile = NULL;
+ }
+ free(mskEnv->qval);
+ free(mskEnv->qsubi);
+ free(mskEnv->qsubj);
+ free(mskEnv->xx);
+ free(mskEnv);
}
-#endif /* MOSEK */
+#endif /* MOSEK */
+/* vim:set shiftwidth=4 ts=8: */
+
+/**
+ *
+ * Authors:
+ * Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * This version is released under the CPL (Common Public License) with
+ * the Graphviz distribution.
+ * A version is also available under the LGPL as part of the Adaptagrams
+ * project: http://sourceforge.net/projects/adaptagrams.
+ * If you make improvements or bug fixes to this code it would be much
+ * appreciated if you could also contribute those changes back to the
+ * Adaptagrams repository.
+ */
+
#ifdef MOSEK
#ifndef _QSOLVE_H_
#define _QSOLVE_H_
void mosek_delete(MosekEnv*);
-#endif // _QSOLVE_H_
-#endif // MOSEK
+#endif /* _QSOLVE_H_ */
+#endif /* MOSEK */
}
}
-static cluster_data* cluster_map(graph_t *mastergraph, graph_t *g) {
+#ifdef IPSEPCOLA
+static cluster_data* cluster_map(graph_t *mastergraph, graph_t *g)
+{
/* search meta-graph to find clusters */
graph_t *mg, *subg;
node_t *mm, *mn;
+ node_t *n;
edge_t *me;
- // array of arrays of node indices in each cluster
+ /* array of arrays of node indices in each cluster */
int **cs,*cn;
int i,j,nclusters=0;
- bool assigned[agnnodes(g)];
- for(i=0;i<agnnodes(g);i++) {
- assigned[i]=false;
- }
+ bool* assigned = N_NEW(agnnodes(g), bool);
cluster_data *cdata = GNEW(cluster_data);
- cdata->ntoplevel=agnnodes(g);
+
+ cdata->ntoplevel = agnnodes(g);
mm = mastergraph->meta_node;
mg = mm->graph;
for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
}
}
cdata->nvars=0;
- cdata->nclusters=nclusters;
- cs=cdata->clusters=N_GNEW(nclusters,int*);
- cn=cdata->clustersizes=N_GNEW(nclusters,int);
- //fprintf(stderr,"search %d clusters...\n",nclusters);
+ cdata->nclusters = nclusters;
+ cs = cdata->clusters = N_GNEW(nclusters,int*);
+ cn = cdata->clustersizes = N_GNEW(nclusters,int);
+ /* fprintf(stderr,"search %d clusters...\n",nclusters); */
for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
mn = me->head;
subg = agusergraph(mn);
/* clusters are processed by separate calls to ordered_edges */
if (!strncmp(subg->name, "cluster", 7)) {
- *cn=agnnodes(subg);
- cdata->nvars+=*cn;
- int *c=*cs++=N_GNEW(*cn++,int);
- node_t *n;
- //fprintf(stderr,"Cluster with %d nodes...\n",agnnodes(subg));
+ int *c;
+
+ *cn = agnnodes(subg);
+ cdata->nvars += *cn;
+ c = *cs++ = N_GNEW(*cn++,int);
+ /* fprintf(stderr,"Cluster with %d nodes...\n",agnnodes(subg)); */
for (n = agfstnode(subg); n; n = agnxtnode(subg, n)) {
node_t *gn;
- int ind=0;
+ int ind = 0;
for (gn = agfstnode(g); gn; gn = agnxtnode(g, gn)) {
if(gn->id==n->id) break;
ind++;
}
- //fprintf(stderr," node=%s, id=%d, ind=%d\n",n->name,n->id,ind);
+ /* fprintf(stderr," node=%s, id=%d, ind=%d\n",n->name,n->id,ind); */
*c++=ind;
assigned[ind]=true;
cdata->ntoplevel--;
}
}
assert(cdata->ntoplevel==agnnodes(g)-cdata->nvars);
+ free (assigned);
return cdata;
}
+
static void freeClusterData(cluster_data *c) {
if(c->nclusters>0) {
free(c->clusters[0]);
}
free(c);
}
+#endif
+
/* user_spline:
* Attempt to use already existing pos info for spline
* Return 1 if successful, 0 otherwise.
#ifdef DIGCOLA
else if (streq(str, "hier"))
mode = MODE_HIER;
+#ifdef IPSEPCOLA
else if (streq(str, "ipsep"))
mode = MODE_IPSEP;
+#endif
#endif
else
agerr(AGWARN,
{
node_t *np, *hp;
int j, e, f;
+ /* if mode is IPSEP make it an in-edge
+ * at both ends, so that an edge constraint won't be generated!
+ */
+ double x = (mode==MODE_IPSEP?-1.0:1.0);
np = graph[i].np;
ND_mark(np) = TRUE;
j = graph[i].edges[e];
hp = graph[j].np;
if (ND_onstack(hp)) { /* back edge: reverse it */
- // if mode is IPSEP make it an in-edge
- // at both ends, so that an edge constraint won't be generated!
- graph[i].edists[e] = mode==MODE_IPSEP?-1.0:1.0;
+ graph[i].edists[e] = x;
for (f = 1; (f < graph[j].nedges) &&(graph[j].edges[f] != i); f++) ;
assert (f < graph[j].nedges);
graph[j].edists[f] = -1.0;
*ewgts++ = 1.0;
#ifdef DIGCOLA
if (haveDir) {
- char *s=agget(ep,"dir");
+ char *s = agget(ep,"dir");
if(s&&!strncmp(s,"none",4)) {
*edists++ = 0;
} else {
* Solve stress using majorization.
* Old neato attributes to incorporate:
* weight?
+ * model will be MODE_MAJOR, MODE_HIER or MODE_IPSEP
*/
static void
majorization(graph_t *mg, graph_t * g, int nv, int mode, int model, int dim, int steps)
int i;
node_t *v;
vtx_data *gp;
- cluster_data *cs;
int init;
}
gp = makeGraphData(g, nv, &ne, mode, model);
- cs=(cluster_data*)cluster_map(mg,g);
-
if (Verbose) {
fprintf(stderr, "%d nodes %.2f sec\n", nv, elapsed_sec());
}
#ifdef DIGCOLA
- if (mode == MODE_HIER || mode == MODE_IPSEP) {
+ if (mode != MODE_MAJOR) {
double lgap = late_double(g, agfindattr(g, "levelsgap"), 0.0, -MAXDOUBLE);
if (mode == MODE_HIER) {
stress_majorization_with_hierarchy(gp, nv, ne, coords, Ndim,
(init == INIT_SELF), model, MaxIter, lgap);
- } else {
+ }
+#ifdef IPSEPCOLA
+ else {
char* str;
ipsep_options opt;
pointf nsize[nv];
- opt.edge_gap=lgap;
+ cluster_data *cs = (cluster_data*)cluster_map(mg,g);
+ opt.edge_gap = lgap;
#ifdef MOSEK
- opt.mosek=0;
-#endif // MOSEK
- opt.noverlap=opt.diredges=0;
- opt.gap.x=opt.gap.y=0;
- opt.nsize=nsize;
- opt.clusters=cs;
+ opt.mosek = 0;
+#endif /* MOSEK */
+ opt.nsize = nsize;
+ opt.clusters = cs;
str = agget(g, "diredgeconstraints");
- if(str && !strncmp(str,"true",4)) {
+ if (mapbool(str)) {
opt.diredges = 1;
if(Verbose)
fprintf(stderr,"Generating Edge Constraints...\n");
- } else if(str && !strncmp(str,"hier",4)) {
+ } else if (str && !strncasecmp(str,"hier",4)) {
opt.diredges = 2;
if(Verbose)
fprintf(stderr,"Generating DiG-CoLa Edge Constraints...\n");
}
+ else opt.diredges = 0;
str = agget(g, "overlapconstraints");
- if(str && !strncmp(str,"true",4)) {
+ if (mapbool(str)) {
opt.noverlap = 1;
if(Verbose)
fprintf(stderr,"Generating Non-overlap Constraints...\n");
- } else if(str && !strncmp(str,"post",4)) {
+ } else if (str && !strncasecmp(str,"post",4)) {
opt.noverlap = 2;
if(Verbose)
fprintf(stderr,"Removing overlaps as postprocess...\n");
}
+ else opt.noverlap = 0;
#ifdef MOSEK
str = agget(g, "mosek");
if(str && !strncmp(str,"true",4)) {
if(Verbose)
fprintf(stderr,"Using Mosek for constraint optimization...\n");
}
-#endif // MOSEK
- if ((str = agget(g, "sep"))) {
- char* c;
- if ((c=strchr(str,','))) {
- i=c-str;
- if(Verbose)
- fprintf(stderr,"gap=%s,%d\n",str,i);
- char s[strlen(str)];
- strncpy(s,str,i);
- s[i]=0;
- opt.gap.x=atof(s);
- strcpy(s,c+1);
- opt.gap.y=atof(s);
+#endif /* MOSEK */
+ if ((str = agget(g, "sep")) &&
+ (i = sscanf(str, "%lf,%lf", &opt.gap.x, &opt.gap.y))) {
+ if (i == 1) opt.gap.y = opt.gap.x;
if(Verbose)
fprintf(stderr,"gap=%f,%f\n",opt.gap.x,opt.gap.y);
- } else {
- opt.gap.x =opt.gap.y = atof(str);
- }
}
+ else opt.gap.x = opt.gap.y = 0;
for (i=0, v = agfstnode(g); v; v = agnxtnode(g, v),i++) {
- nsize[i].x=ND_width(v);
- nsize[i].y=ND_height(v);
+ nsize[i].x = ND_width(v);
+ nsize[i].y = ND_height(v);
}
stress_majorization_cola(gp, nv, ne, coords, Ndim, model, MaxIter, &opt);
+ freeClusterData(cs);
}
+#endif
}
else
#endif
}
}
freeGraphData(gp);
- freeClusterData(cs);
free(coords[0]);
free(coords);
}
free_scan_graph(gc);
agdelete(g, gc);
}
+#ifdef IPSEPCOLA
{
graph_t *mg, *subg;
node_t *mm, *mn;
}
}
}
+#endif
} else {
neatoLayout(g, g, layoutMode, model);
adjustNodes(g);
(verts[0].y == verts[3].y) && (verts[1].y == verts[2].y));
}
-#define DFLT_SAMPLE 20
-
static Point makeScaledPoint(int x, int y)
{
Point rv;
+/* vim:set shiftwidth=4 ts=8: */
/**********************************************************
* This software is part of the graphviz package *
* http://www.graphviz.org/ *
#define quad_prog_tol 1e-2
-float**
-unpackMatrix(float * packedMat, int n)
+float **unpackMatrix(float *packedMat, int n)
{
float **mat;
- int i,j,k;
-
- mat = N_GNEW(n, float*);
- mat[0] = N_GNEW(n*n, float);
- set_vector_valf(n*n, 0, mat[0]);
- for (i=1; i<n; i++) {
- mat[i] = mat[0]+i*n;
- }
-
- for (i=0,k=0; i<n; i++) {
- for (j=i; j<n; j++,k++) {
- mat[j][i]=mat[i][j]=packedMat[k];
- }
+ int i, j, k;
+
+ mat = N_GNEW(n, float *);
+ mat[0] = N_GNEW(n * n, float);
+ set_vector_valf(n * n, 0, mat[0]);
+ for (i = 1; i < n; i++) {
+ mat[i] = mat[0] + i * n;
+ }
+
+ for (i = 0, k = 0; i < n; i++) {
+ for (j = i; j < n; j++, k++) {
+ mat[j][i] = mat[i][j] = packedMat[k];
}
+ }
return mat;
}
-static void
-ensureMonotonicOrdering(float * place, int n, int * ordering)
+static void ensureMonotonicOrdering(float *place, int n, int *ordering)
{
- /* ensure that 'ordering' is monotonically increasing by 'place', */
- /* this also implies that levels are separated in the initial layout */
- int i, node;
- float lower_bound = place[ordering[0]];
- for (i=1; i<n; i++) {
- node = ordering[i];
- if (place[node]<lower_bound) {
- place[node]=lower_bound;
- }
- lower_bound = place[node];
+ /* ensure that 'ordering' is monotonically increasing by 'place', */
+ /* this also implies that levels are separated in the initial layout */
+ int i, node;
+ float lower_bound = place[ordering[0]];
+ for (i = 1; i < n; i++) {
+ node = ordering[i];
+ if (place[node] < lower_bound) {
+ place[node] = lower_bound;
}
+ lower_bound = place[node];
+ }
}
-
+
static void
-ensureMonotonicOrderingWithGaps(float* place, int n, int* ordering,
- int* levels, int num_levels, float levels_gap)
+ensureMonotonicOrderingWithGaps(float *place, int n, int *ordering,
+ int *levels, int num_levels,
+ float levels_gap)
{
- /* ensure that levels are separated in the initial layout and that
+ /* ensure that levels are separated in the initial layout and that
* places are monotonic increasing within layer
*/
-
- int i;
- int node, level, max_in_level;
- float lower_bound = (float)-1e9;
-
- level = -1; max_in_level=0;
- for (i=0; i<n; i++) {
- if (i>=max_in_level) {
- /* we are entering a new level */
- level++;
- if (level==num_levels) {
- /* last_level */
- max_in_level=n;
- }
- else {
- max_in_level = levels[level];
- }
- lower_bound = i>0 ? place[ordering[i-1]]+levels_gap : (float)-1e9;
- quicksort_placef(place, ordering, i, max_in_level-1);
- }
- node = ordering[i];
- if (place[node]<lower_bound) {
- place[node]=lower_bound;
- }
+ int i;
+ int node, level, max_in_level;
+ float lower_bound = (float) -1e9;
+
+ level = -1;
+ max_in_level = 0;
+ for (i = 0; i < n; i++) {
+ if (i >= max_in_level) {
+ /* we are entering a new level */
+ level++;
+ if (level == num_levels) {
+ /* last_level */
+ max_in_level = n;
+ } else {
+ max_in_level = levels[level];
+ }
+ lower_bound =
+ i > 0 ? place[ordering[i - 1]] + levels_gap : (float) -1e9;
+ quicksort_placef(place, ordering, i, max_in_level - 1);
+ }
+
+ node = ordering[i];
+ if (place[node] < lower_bound) {
+ place[node] = lower_bound;
}
+ }
}
static void
-computeHierarchyBoundaries(float* place, int n, int* ordering, int* levels,
- int num_levels, float* hierarchy_boundaries)
+computeHierarchyBoundaries(float *place, int n, int *ordering, int *levels,
+ int num_levels, float *hierarchy_boundaries)
{
- int i;
- for (i=0; i<num_levels; i++) {
- hierarchy_boundaries[i] = place[ordering[levels[i]-1]];
- }
+ int i;
+ for (i = 0; i < num_levels; i++) {
+ hierarchy_boundaries[i] = place[ordering[levels[i] - 1]];
+ }
}
-int
-constrained_majorization_new(CMajEnv *e, float * b, float **coords,
- int cur_axis, int dims, int max_iterations,
- float* hierarchy_boundaries, float levels_gap)
+int
+constrained_majorization_new(CMajEnv * e, float *b, float **coords,
+ int cur_axis, int dims, int max_iterations,
+ float *hierarchy_boundaries, float levels_gap)
{
- int n=e->n;
- float *place=coords[cur_axis];
- float **lap=e->A;
- int * ordering=e->ordering;
- int * levels=e->levels;
- int num_levels=e->num_levels;
- int i,j;
- float new_place_i;
- bool converged=false;
- float upper_bound, lower_bound;
- int node;
- int left, right;
- float cur_place;
- float des_place_block;
- float block_deg;
- float toBlockConnectivity;
- float * lap_node;
- int block_len;
- int first_next_level;
- int level=-1, max_in_level=0;
- float* desired_place;
- float* prefix_desired_place;
- float* suffix_desired_place;
- int* block;
- int* lev;
+ int n = e->n;
+ float *place = coords[cur_axis];
+ float **lap = e->A;
+ int *ordering = e->ordering;
+ int *levels = e->levels;
+ int num_levels = e->num_levels;
+ int i, j;
+ float new_place_i;
+ bool converged = false;
+ float upper_bound, lower_bound;
+ int node;
+ int left, right;
+ float cur_place;
+ float des_place_block;
+ float block_deg;
+ float toBlockConnectivity;
+ float *lap_node;
+ int block_len;
+ int first_next_level;
+ int level = -1, max_in_level = 0;
+ float *desired_place;
+ float *prefix_desired_place;
+ float *suffix_desired_place;
+ int *block;
+ int *lev;
int counter;
- if (max_iterations<=0) {
- return 0;
- }
- if (levels_gap!=0) {
- return constrained_majorization_new_with_gaps(e, b, coords,cur_axis,dims, max_iterations, hierarchy_boundaries, levels_gap);
- }
-
- /* ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels); */
- ensureMonotonicOrdering(place, n, ordering);
- /* it is important that in 'ordering' nodes are always sorted by layers,
+ if (max_iterations <= 0) {
+ return 0;
+ }
+ if (levels_gap != 0) {
+ return constrained_majorization_new_with_gaps(e, b, coords,
+ cur_axis, dims,
+ max_iterations,
+ hierarchy_boundaries,
+ levels_gap);
+ }
+
+ /* ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels); */
+ ensureMonotonicOrdering(place, n, ordering);
+ /* it is important that in 'ordering' nodes are always sorted by layers,
* and within a layer by 'place'
*/
-
+
/* the desired place of each individual node in the current block */
- desired_place = e->fArray1;
+ desired_place = e->fArray1;
/* the desired place of each prefix of current block */
- prefix_desired_place = e->fArray2;
+ prefix_desired_place = e->fArray2;
/* the desired place of each suffix of current block */
- suffix_desired_place = e->fArray3;
+ suffix_desired_place = e->fArray3;
/* current block (nodes connected by active constraints) */
- block = e->iArray1;
-
- lev = e->iArray2; /* level of each node */
- for (i=0; i<n; i++) {
- if (i>=max_in_level) {
- /* we are entering a new level */
- level++;
- if (level==num_levels) {
- /* last_level */
- max_in_level=n;
- }
- else {
- max_in_level = levels[level];
- }
+ block = e->iArray1;
+
+ lev = e->iArray2; /* level of each node */
+ for (i = 0; i < n; i++) {
+ if (i >= max_in_level) {
+ /* we are entering a new level */
+ level++;
+ if (level == num_levels) {
+ /* last_level */
+ max_in_level = n;
+ } else {
+ max_in_level = levels[level];
+ }
+ }
+ node = ordering[i];
+ lev[node] = level;
+ }
+
+ for (counter = 0; counter < max_iterations && !converged; counter++) {
+ converged = true;
+ lower_bound = -1e9; /* no lower bound for first level */
+ for (left = 0; left < n; left = right) {
+ int best_i;
+ double max_movement;
+ double movement;
+ float prefix_des_place, suffix_des_place;
+ /* compute a block 'ordering[left]...ordering[right-1]' of
+ * nodes with the same coordinate:
+ */
+ cur_place = place[ordering[left]];
+ for (right = left + 1; right < n; right++) {
+ if (place[ordering[right]] != cur_place) {
+ break;
}
+ }
+
+ /* compute desired place of nodes in block: */
+ for (i = left; i < right; i++) {
node = ordering[i];
- lev[node]=level;
- }
+ new_place_i = -b[node];
+ lap_node = lap[node];
+ for (j = 0; j < n; j++) {
+ if (j == node) {
+ continue;
+ }
+ new_place_i += lap_node[j] * place[j];
+ }
+ desired_place[node] = new_place_i / (-lap_node[node]);
+ }
+
+ /* reorder block by levels, and within levels by "relaxed" desired position */
+ block_len = 0;
+ first_next_level = 0;
+ for (i = left; i < right; i = first_next_level) {
+ level = lev[ordering[i]];
+ if (level == num_levels) {
+ /* last_level */
+ first_next_level = right;
+ } else {
+ first_next_level = MIN(right, levels[level]);
+ }
- for (counter=0; counter<max_iterations&&!converged; counter++) {
- converged=true;
- lower_bound = -1e9; /* no lower bound for first level */
- for (left=0; left<n; left=right) {
- int best_i;
- double max_movement;
- double movement;
- float prefix_des_place, suffix_des_place;
- /* compute a block 'ordering[left]...ordering[right-1]' of
- * nodes with the same coordinate:
- */
- cur_place = place[ordering[left]];
- for (right=left+1; right<n; right++) {
- if (place[ordering[right]]!=cur_place) {
- break;
- }
+ /* First, collect all nodes with desired places smaller than current place */
+ for (j = i; j < first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] < cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Second, collect all nodes with desired places equal to current place */
+ for (j = i; j < first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] == cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Third, collect all nodes with desired places greater than current place */
+ for (j = i; j < first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] > cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ }
+
+ /* loop through block and compute desired places of its prefixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i = 0; i < block_len; i++) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j = 0; j < i; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity *= 2;
+ /* update block stats */
+ des_place_block =
+ (block_deg * des_place_block +
+ (-lap_node[node]) * desired_place[node] +
+ toBlockConnectivity * cur_place) / (block_deg -
+ lap_node[node] +
+ toBlockConnectivity);
+ prefix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity - lap_node[node];
+ }
+
+ /* loop through block and compute desired places of its suffixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i = block_len - 1; i >= 0; i--) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j = i + 1; j < block_len; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity *= 2;
+ /* update block stats */
+ des_place_block =
+ (block_deg * des_place_block +
+ (-lap_node[node]) * desired_place[node] +
+ toBlockConnectivity * cur_place) / (block_deg -
+ lap_node[node] +
+ toBlockConnectivity);
+ suffix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity - lap_node[node];
+ }
+
+
+ /* now, find best place to split block */
+ best_i = -1;
+ max_movement = 0;
+ for (i = 0; i < block_len; i++) {
+ suffix_des_place = suffix_desired_place[i];
+ prefix_des_place =
+ i > 0 ? prefix_desired_place[i - 1] : suffix_des_place;
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place < prefix_des_place) {
+ if (suffix_des_place < cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
}
-
- /* compute desired place of nodes in block: */
- for (i=left; i<right; i++) {
- node = ordering[i];
- new_place_i = -b[node];
- lap_node = lap[node];
- for (j=0; j<n; j++) {
- if (j==node) {
- continue;
- }
- new_place_i += lap_node[j]*place[j];
- }
- desired_place[node] = new_place_i/(-lap_node[node]);
+ suffix_des_place = prefix_des_place;
+ } else if (prefix_des_place > cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
+ movement =
+ (block_len - i) * fabs(suffix_des_place - cur_place) +
+ i * fabs(prefix_des_place - cur_place);
+ if (movement > max_movement) {
+ max_movement = movement;
+ best_i = i;
+ }
+ }
+ /* Actually move prefix and suffix */
+ if (best_i >= 0) {
+ suffix_des_place = suffix_desired_place[best_i];
+ prefix_des_place =
+ best_i >
+ 0 ? prefix_desired_place[best_i -
+ 1] : suffix_des_place;
+
+ /* compute right border of feasible move */
+ if (right >= n) {
+ /* no nodes after current block */
+ upper_bound = 1e9;
+ } else {
+ upper_bound = place[ordering[right]];
+ }
+ suffix_des_place = MIN(suffix_des_place, upper_bound);
+ prefix_des_place = MAX(prefix_des_place, lower_bound);
+
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place < prefix_des_place) {
+ if (suffix_des_place < cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
}
+ suffix_des_place = prefix_des_place;
+ } else if (prefix_des_place > cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
- /* reorder block by levels, and within levels by "relaxed" desired position */
- block_len = 0;
- first_next_level=0;
- for (i=left; i<right; i=first_next_level) {
- level = lev[ordering[i]];
- if (level==num_levels) {
- /* last_level */
- first_next_level=right;
- }
- else {
- first_next_level = MIN(right,levels[level]);
- }
-
- /* First, collect all nodes with desired places smaller than current place */
- for (j=i; j<first_next_level; j++) {
- node = ordering[j];
- if (desired_place[node] < cur_place) {
- block[block_len++] = node;
- }
- }
- /* Second, collect all nodes with desired places equal to current place */
- for (j=i; j<first_next_level; j++) {
- node = ordering[j];
- if (desired_place[node] == cur_place) {
- block[block_len++] = node;
- }
- }
- /* Third, collect all nodes with desired places greater than current place */
- for (j=i; j<first_next_level; j++) {
- node = ordering[j];
- if (desired_place[node] > cur_place) {
- block[block_len++] = node;
- }
- }
- }
-
- /* loop through block and compute desired places of its prefixes */
- des_place_block = 0;
- block_deg = 0;
- for (i=0; i<block_len; i++) {
- node = block[i];
- toBlockConnectivity = 0;
- lap_node = lap[node];
- for (j=0; j<i; j++) {
- toBlockConnectivity -= lap_node[block[j]];
- }
- toBlockConnectivity*=2;
- /* update block stats */
- des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
- prefix_desired_place[i] = des_place_block;
- block_deg += toBlockConnectivity-lap_node[node];
- }
+ /* move prefix: */
+ for (i = 0; i < best_i; i++) {
+ place[block[i]] = prefix_des_place;
+ }
+ /* move suffix: */
+ for (i = best_i; i < block_len; i++) {
+ place[block[i]] = suffix_des_place;
+ }
- /* loop through block and compute desired places of its suffixes */
- des_place_block = 0;
- block_deg = 0;
- for (i=block_len-1; i>=0; i--) {
- node = block[i];
- toBlockConnectivity = 0;
- lap_node = lap[node];
- for (j=i+1; j<block_len; j++) {
- toBlockConnectivity -= lap_node[block[j]];
- }
- toBlockConnectivity*=2;
- /* update block stats */
- des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
- suffix_desired_place[i] = des_place_block;
- block_deg += toBlockConnectivity-lap_node[node];
- }
+ lower_bound = suffix_des_place; /* lower bound for next block */
-
- /* now, find best place to split block */
- best_i = -1;
- max_movement = 0;
- for (i=0; i<block_len; i++) {
- suffix_des_place = suffix_desired_place[i];
- prefix_des_place = i>0 ? prefix_desired_place[i-1] : suffix_des_place;
- /* limit moves to ensure that the prefix is placed before the suffix */
- if (suffix_des_place<prefix_des_place) {
- if (suffix_des_place<cur_place) {
- if (prefix_des_place > cur_place) {
- prefix_des_place = cur_place;
- }
- suffix_des_place = prefix_des_place;
- }
- else if (prefix_des_place>cur_place) {
- prefix_des_place = suffix_des_place;
- }
- }
- movement = (block_len-i)*fabs(suffix_des_place-cur_place) + i*fabs(prefix_des_place-cur_place);
- if (movement>max_movement) {
- max_movement = movement;
- best_i = i;
- }
- }
- /* Actually move prefix and suffix */
- if (best_i>=0) {
- suffix_des_place = suffix_desired_place[best_i];
- prefix_des_place = best_i>0 ? prefix_desired_place[best_i-1] : suffix_des_place;
-
- /* compute right border of feasible move */
- if (right>=n) {
- /* no nodes after current block */
- upper_bound = 1e9;
- }
- else {
- upper_bound = place[ordering[right]];
- }
- suffix_des_place = MIN(suffix_des_place, upper_bound);
- prefix_des_place = MAX(prefix_des_place, lower_bound);
-
- /* limit moves to ensure that the prefix is placed before the suffix */
- if (suffix_des_place<prefix_des_place) {
- if (suffix_des_place<cur_place) {
- if (prefix_des_place > cur_place) {
- prefix_des_place = cur_place;
- }
- suffix_des_place = prefix_des_place;
- }
- else if (prefix_des_place>cur_place) {
- prefix_des_place = suffix_des_place;
- }
- }
-
- /* move prefix: */
- for (i=0; i<best_i; i++) {
- place[block[i]] = prefix_des_place;
- }
- /* move suffix: */
- for (i=best_i; i<block_len; i++) {
- place[block[i]] = suffix_des_place;
- }
-
- lower_bound = suffix_des_place; // lower bound for next block
-
- /* reorder 'ordering' to reflect change of places
- * Note that it is enough to reorder the level where
- * the split was done
- */
+ /* reorder 'ordering' to reflect change of places
+ * Note that it is enough to reorder the level where
+ * the split was done
+ */
#if 0
- int max_in_level, min_in_level;
-
- level = lev[best_i];
- if (level==num_levels) {
- // last_level
- max_in_level = MIN(right,n);
- }
- else {
- max_in_level = MIN(right,levels[level]);
- }
- if (level==0) {
- // first level
- min_in_level = MAX(left,0);
- }
- else {
- min_in_level = MAX(left,levels[level-1]);
- }
+ int max_in_level, min_in_level;
+
+ level = lev[best_i];
+ if (level == num_levels) {
+ /* last_level */
+ max_in_level = MIN(right, n);
+ } else {
+ max_in_level = MIN(right, levels[level]);
+ }
+ if (level == 0) {
+ /* first level */
+ min_in_level = MAX(left, 0);
+ } else {
+ min_in_level = MAX(left, levels[level - 1]);
+ }
#endif
- for (i=left; i<right; i++) {
- ordering[i] = block[i-left];
- }
- converged = converged && fabs(prefix_des_place-cur_place)<quad_prog_tol
- && fabs(suffix_des_place-cur_place)<quad_prog_tol;
-
- }
- else {
- /* no movement */
- lower_bound = cur_place; // lower bound for next block
- }
-
+ for (i = left; i < right; i++) {
+ ordering[i] = block[i - left];
}
+ converged = converged
+ && fabs(prefix_des_place - cur_place) < quad_prog_tol
+ && fabs(suffix_des_place - cur_place) < quad_prog_tol;
+
+ } else {
+ /* no movement */
+ lower_bound = cur_place; /* lower bound for next block */
+ }
+
}
+ }
- computeHierarchyBoundaries(place, n, ordering, levels, num_levels, hierarchy_boundaries);
+ computeHierarchyBoundaries(place, n, ordering, levels, num_levels,
+ hierarchy_boundaries);
- return counter;
+ return counter;
}
+#ifdef IPSEPCOLA
static float *place;
-static int compare_incr(const void *a, const void *b) {
- if(place[*(int*)a]>place[*(int*)b]) {
- return 1;
- } else if(place[*(int*)a]<place[*(int*)b]) {
- return -1;
- }
- return 0;
+static int compare_incr(const void *a, const void *b)
+{
+ if (place[*(int *) a] > place[*(int *) b]) {
+ return 1;
+ } else if (place[*(int *) a] < place[*(int *) b]) {
+ return -1;
+ }
+ return 0;
}
+
/*
While not converged: move everything towards the optimum, then satisfy constraints with as little displacement as possible.
Returns number of iterations before convergence.
*/
-int constrained_majorization_gradient_projection(CMajEnv *e,
- float * b, float ** coords, int ndims, int cur_axis, int max_iterations,
- float * hierarchy_boundaries,float levels_gap) {
-
- int i,j,counter;
- if(max_iterations==0) return 0;
-
- int * ordering=e->ordering;
- int * levels=e->levels;
- int num_levels=e->num_levels;
- bool converged=false;
- place = coords[cur_axis];
- float *g = e->fArray1;
- float *old_place = e->fArray2;
- float *d = e->fArray4;
+int constrained_majorization_gradient_projection(CMajEnv * e,
+ float *b, float **coords,
+ int ndims, int cur_axis,
+ int max_iterations,
+ float
+ *hierarchy_boundaries,
+ float levels_gap)
+{
+
+ int i, j, counter;
+ int *ordering = e->ordering;
+ int *levels = e->levels;
+ int num_levels = e->num_levels;
+ bool converged = false;
+ float *g = e->fArray1;
+ float *old_place = e->fArray2;
+ float *d = e->fArray4;
+ float test = 0, tmptest = 0;
+ float beta;
+
+ if (max_iterations == 0)
+ return 0;
+ place = coords[cur_axis];
#ifdef CONMAJ_LOGGING
- double prev_stress=0;
- static int call_no=0;
- for (i=0; i<e->n; i++) {
- prev_stress += 2*b[i]*place[i];
- for (j=0; j<e->n; j++) {
- prev_stress -= e->A[i][j]*place[j]*place[i];
- }
+ double prev_stress = 0;
+ static int call_no = 0;
+ for (i = 0; i < e->n; i++) {
+ prev_stress += 2 * b[i] * place[i];
+ for (j = 0; j < e->n; j++) {
+ prev_stress -= e->A[i][j] * place[j] * place[i];
}
- FILE *logfile = fopen("constrained_majorization_log","a");
+ }
+ FILE *logfile = fopen("constrained_majorization_log", "a");
- fprintf(logfile,"grad proj %d: stress=%f\n",call_no,prev_stress);
+ fprintf(logfile, "grad proj %d: stress=%f\n", call_no, prev_stress);
#endif
- for (counter=0; counter<max_iterations&&!converged; counter++) {
- converged=true;
- // find steepest descent direction
- for (i=0; i<e->n; i++) {
- old_place[i]=place[i];
- g[i] = 2*b[i];
- for (j=0; j<e->n; j++) {
- g[i] -= 2*e->A[i][j]*place[j];
- }
- }
- float numerator = 0, denominator = 0, r;
- for (i=0; i<e->n; i++) {
- numerator += g[i]*g[i];
- r=0;
- for (j=0; j<e->n; j++) {
- r += 2*e->A[i][j]*g[j];
- }
- denominator -= r*g[i];
- }
- float alpha = numerator/denominator;
- for (i=0; i<e->n; i++) {
- if(alpha > 0 && alpha < 1000) {
- place[i]-=alpha*g[i];
- }
- }
- if(num_levels)
- qsort((int *)ordering, (size_t)levels[0], sizeof(int), compare_incr );
- //project to constraint boundary
- for (i=0; i<num_levels; i++) {
- int endOfLevel=i==num_levels-1?e->n:levels[i+1];
- // ensure monotic increase in position within levels
- qsort((int *)ordering+levels[i], (size_t)endOfLevel-levels[i], sizeof(int), compare_incr );
- // If there are overlapping levels find offending nodes and place at average position
- int ui=levels[i],li=ui-1;
- int l = ordering[li--],u = ordering[ui++];
- if(place[l]+levels_gap>place[u]) {
- float sum=place[l]+place[u]-levels_gap*(e->lev[l]+e->lev[u]),w=2;
- float avgPos=sum/w;
- float pos;
- bool finished;
- do {
- finished=true;
- if(ui<endOfLevel) {
- u=ordering[ui];
- pos=place[u]-levels_gap*e->lev[u];
- if(pos<avgPos) {
- ui++; w++;
- sum+=pos;
- avgPos=sum/w;
- finished=false;
- }
- }
-
- if(li>=0) {
- l=ordering[li];
- pos=place[l]-levels_gap*e->lev[l];
- if(pos>avgPos) {
- li--; w++;
- sum+=pos;
- avgPos=sum/w;
- finished=false;
- }
- }
- } while (!finished);
- for(j=li+1;j<ui;j++) {
- place[ordering[j]]=avgPos+levels_gap*e->lev[ordering[j]];
- }
+ for (counter = 0; counter < max_iterations && !converged; counter++) {
+ float alpha;
+ float numerator = 0, denominator = 0, r;
+ converged = true;
+ /* find steepest descent direction */
+ for (i = 0; i < e->n; i++) {
+ old_place[i] = place[i];
+ g[i] = 2 * b[i];
+ for (j = 0; j < e->n; j++) {
+ g[i] -= 2 * e->A[i][j] * place[j];
+ }
+ }
+ for (i = 0; i < e->n; i++) {
+ numerator += g[i] * g[i];
+ r = 0;
+ for (j = 0; j < e->n; j++) {
+ r += 2 * e->A[i][j] * g[j];
+ }
+ denominator -= r * g[i];
+ }
+ alpha = numerator / denominator;
+ for (i = 0; i < e->n; i++) {
+ if (alpha > 0 && alpha < 1000) {
+ place[i] -= alpha * g[i];
+ }
+ }
+ if (num_levels)
+ qsort((int *) ordering, (size_t) levels[0], sizeof(int),
+ compare_incr);
+ /* project to constraint boundary */
+ for (i = 0; i < num_levels; i++) {
+ int endOfLevel = i == num_levels - 1 ? e->n : levels[i + 1];
+ int ui, li, u, l;
+
+ /* ensure monotic increase in position within levels */
+ qsort((int *) ordering + levels[i],
+ (size_t) endOfLevel - levels[i], sizeof(int),
+ compare_incr);
+ /* If there are overlapping levels find offending nodes and place at average position */
+ ui = levels[i]; li = ui - 1;
+ l = ordering[li--]; u = ordering[ui++];
+ if (place[l] + levels_gap > place[u]) {
+ float sum =
+ place[l] + place[u] - levels_gap * (e->lev[l] +
+ e->lev[u]), w = 2;
+ float avgPos = sum / w;
+ float pos;
+ bool finished;
+ do {
+ finished = true;
+ if (ui < endOfLevel) {
+ u = ordering[ui];
+ pos = place[u] - levels_gap * e->lev[u];
+ if (pos < avgPos) {
+ ui++;
+ w++;
+ sum += pos;
+ avgPos = sum / w;
+ finished = false;
}
- }
- // set place to the intersection of old_place-g and boundary and compute d, the vector from intersection pnt to projection pnt
- for (i=0; i<e->n; i++) {
- d[i]=place[i]-old_place[i];
- }
- // now compute beta
- numerator = 0, denominator = 0;
- for (i=0; i<e->n; i++) {
- numerator += g[i]*d[i];
- r=0;
- for (j=0; j<e->n; j++) {
- r += 2*e->A[i][j]*d[j];
+ }
+
+ if (li >= 0) {
+ l = ordering[li];
+ pos = place[l] - levels_gap * e->lev[l];
+ if (pos > avgPos) {
+ li--;
+ w++;
+ sum += pos;
+ avgPos = sum / w;
+ finished = false;
}
- denominator += r*d[i];
+ }
+ } while (!finished);
+ for (j = li + 1; j < ui; j++) {
+ place[ordering[j]] =
+ avgPos + levels_gap * e->lev[ordering[j]];
}
- float beta = numerator/denominator;
-
- float test=0,tmptest=0;
- for (i=0; i<e->n; i++) {
- if(beta>0&&beta<1.0) {
- place[i]=old_place[i]+beta*d[i];
- }
- tmptest = fabs(place[i]-old_place[i]);
- if(test<tmptest) test=tmptest;
- }
- computeHierarchyBoundaries(place, e->n, ordering, levels, num_levels, hierarchy_boundaries);
- /*
- if(num_levels)
- qsort((int *)ordering, (size_t)levels[0], sizeof(int), compare_incr );
- for (i=0; i<num_levels; i++) {
- int endOfLevel=i==num_levels-1?e->n:levels[i+1];
- // ensure monotic increase in position within levels
- qsort((int *)ordering+levels[i], (size_t)endOfLevel-levels[i], sizeof(int), compare_incr );
- // If there are overlapping levels find offending nodes and place at average position
- int l = ordering[levels[i]-1],u = ordering[levels[i]];
- //assert(place[l]+levels_gap<=place[u]+0.00001);
- }
- */
+ }
+ }
+ /* set place to the intersection of old_place-g and boundary and compute d, the vector from intersection pnt to projection pnt */
+ for (i = 0; i < e->n; i++) {
+ d[i] = place[i] - old_place[i];
+ }
+ /* now compute beta */
+ numerator = 0, denominator = 0;
+ for (i = 0; i < e->n; i++) {
+ numerator += g[i] * d[i];
+ r = 0;
+ for (j = 0; j < e->n; j++) {
+ r += 2 * e->A[i][j] * d[j];
+ }
+ denominator += r * d[i];
+ }
+ beta = numerator / denominator;
+
+ for (i = 0; i < e->n; i++) {
+ if (beta > 0 && beta < 1.0) {
+ place[i] = old_place[i] + beta * d[i];
+ }
+ tmptest = fabs(place[i] - old_place[i]);
+ if (test < tmptest)
+ test = tmptest;
+ }
+ computeHierarchyBoundaries(place, e->n, ordering, levels,
+ num_levels, hierarchy_boundaries);
+#if 0
+ if (num_levels)
+ qsort((int *) ordering, (size_t) levels[0], sizeof(int),
+ compare_incr);
+ for (i = 0; i < num_levels; i++) {
+ int endOfLevel = i == num_levels - 1 ? e->n : levels[i + 1];
+ /* ensure monotic increase in position within levels */
+ qsort((int *) ordering + levels[i],
+ (size_t) endOfLevel - levels[i], sizeof(int),
+ compare_incr);
+ /* If there are overlapping levels find offending nodes and place at average position */
+ int l = ordering[levels[i] - 1], u = ordering[levels[i]];
+ /* assert(place[l]+levels_gap<=place[u]+0.00001); */
+ }
+#endif
#ifdef CONMAJ_LOGGING
- double stress=0;
- for (i=0; i<e->n; i++) {
- stress += 2*b[i]*place[i];
- for (j=0; j<e->n; j++) {
- stress -= e->A[i][j]*place[j]*place[i];
- }
- }
- fprintf(logfile,"%d: stress=%f, test=%f, %s\n",call_no,stress,test,(stress>=prev_stress)?"No Improvement":"");
- prev_stress=stress;
+ double stress = 0;
+ for (i = 0; i < e->n; i++) {
+ stress += 2 * b[i] * place[i];
+ for (j = 0; j < e->n; j++) {
+ stress -= e->A[i][j] * place[j] * place[i];
+ }
+ }
+ fprintf(logfile, "%d: stress=%f, test=%f, %s\n", call_no, stress,
+ test, (stress >= prev_stress) ? "No Improvement" : "");
+ prev_stress = stress;
#endif
- if(test>quad_prog_tol) {
- converged=false;
- }
+ if (test > quad_prog_tol) {
+ converged = false;
}
+ }
#ifdef CONMAJ_LOGGING
- call_no++;
- fclose(logfile);
+ call_no++;
+ fclose(logfile);
#endif
- return counter;
+ return counter;
}
+#endif
-int
-constrained_majorization_new_with_gaps(CMajEnv *e, float* b, float** coords,
- int ndims, int cur_axis, int max_iterations, float* hierarchy_boundaries,
- float levels_gap)
+int
+constrained_majorization_new_with_gaps(CMajEnv * e, float *b,
+ float **coords, int ndims,
+ int cur_axis, int max_iterations,
+ float *hierarchy_boundaries,
+ float levels_gap)
{
- float *place=coords[cur_axis];
- int i,j;
- int n=e->n;
- float** lap = e->A;
- int* ordering = e->ordering;
- int* levels = e->levels;
- int num_levels = e->num_levels;
- float new_place_i;
- bool converged=false;
- float upper_bound, lower_bound;
- int node;
- int left, right;
- float cur_place;
- float des_place_block;
- float block_deg;
- float toBlockConnectivity;
- float * lap_node;
- float * desired_place;
- float * prefix_desired_place;
- float * suffix_desired_place;
- int* block;
- int block_len;
- int first_next_level;
- int * lev;
- int level=-1, max_in_level=0;
+ float *place = coords[cur_axis];
+ int i, j;
+ int n = e->n;
+ float **lap = e->A;
+ int *ordering = e->ordering;
+ int *levels = e->levels;
+ int num_levels = e->num_levels;
+ float new_place_i;
+ bool converged = false;
+ float upper_bound, lower_bound;
+ int node;
+ int left, right;
+ float cur_place;
+ float des_place_block;
+ float block_deg;
+ float toBlockConnectivity;
+ float *lap_node;
+ float *desired_place;
+ float *prefix_desired_place;
+ float *suffix_desired_place;
+ int *block;
+ int block_len;
+ int first_next_level;
+ int *lev;
+ int level = -1, max_in_level = 0;
int counter;
- float* gap;
- float total_gap, target_place;
+ float *gap;
+ float total_gap, target_place;
- if (max_iterations<=0) {
- return 0;
- }
+ if (max_iterations <= 0) {
+ return 0;
+ }
- ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels, levels_gap);
- /* it is important that in 'ordering' nodes are always sorted by layers,
+ ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels,
+ levels_gap);
+ /* it is important that in 'ordering' nodes are always sorted by layers,
* and within a layer by 'place'
*/
-
+
/* the desired place of each individual node in the current block */
- desired_place = e->fArray1;
+ desired_place = e->fArray1;
/* the desired place of each prefix of current block */
- prefix_desired_place = e->fArray2;
+ prefix_desired_place = e->fArray2;
/* the desired place of each suffix of current block */
- suffix_desired_place = e->fArray3;
+ suffix_desired_place = e->fArray3;
/* current block (nodes connected by active constraints) */
- block = e->iArray1;
-
- lev = e->iArray2; /* level of each node */
- for (i=0; i<n; i++) {
- if (i>=max_in_level) {
- /* we are entering a new level */
- level++;
- if (level==num_levels) {
- /* last_level */
- max_in_level=n;
- }
- else {
- max_in_level = levels[level];
- }
- }
- node = ordering[i];
- lev[node]=level;
+ block = e->iArray1;
+
+ lev = e->iArray2; /* level of each node */
+ for (i = 0; i < n; i++) {
+ if (i >= max_in_level) {
+ /* we are entering a new level */
+ level++;
+ if (level == num_levels) {
+ /* last_level */
+ max_in_level = n;
+ } else {
+ max_in_level = levels[level];
+ }
}
+ node = ordering[i];
+ lev[node] = level;
+ }
/* displacement of block's nodes from block's reference point */
- gap = e->fArray4;
-
- for (counter=0; counter<max_iterations&&!converged; counter++) {
- converged=true;
- lower_bound = -1e9; // no lower bound for first level
- for (left=0; left<n; left=right) {
- int best_i;
- double max_movement;
- double movement;
- float prefix_des_place, suffix_des_place;
- /* compute a block 'ordering[left]...ordering[right-1]' of
- * nodes connected with active constraints
- */
- cur_place = place[ordering[left]];
- total_gap = 0;
- target_place = cur_place;
- gap[ordering[left]] = 0;
- for (right=left+1; right<n; right++) {
- if (lev[right]>lev[right-1]) {
- /* we are entering a new level */
- target_place += levels_gap; /* note that 'levels_gap' may be negative */
- total_gap += levels_gap;
- }
- node = ordering[right];
+ gap = e->fArray4;
+
+ for (counter = 0; counter < max_iterations && !converged; counter++) {
+ converged = true;
+ lower_bound = -1e9; /* no lower bound for first level */
+ for (left = 0; left < n; left = right) {
+ int best_i;
+ double max_movement;
+ double movement;
+ float prefix_des_place, suffix_des_place;
+ /* compute a block 'ordering[left]...ordering[right-1]' of
+ * nodes connected with active constraints
+ */
+ cur_place = place[ordering[left]];
+ total_gap = 0;
+ target_place = cur_place;
+ gap[ordering[left]] = 0;
+ for (right = left + 1; right < n; right++) {
+ if (lev[right] > lev[right - 1]) {
+ /* we are entering a new level */
+ target_place += levels_gap; /* note that 'levels_gap' may be negative */
+ total_gap += levels_gap;
+ }
+ node = ordering[right];
#if 0
- if (place[node]!=target_place)
+ if (place[node] != target_place)
#endif
- /* not sure if this is better than 'place[node]!=target_place' */
- if (fabs(place[node]-target_place)>1e-9)
- {
- break;
- }
- gap[node] = place[node]-cur_place;
- }
-
- /* compute desired place of block's reference point according
- * to each node in the block:
- */
- for (i=left; i<right; i++) {
- node = ordering[i];
- new_place_i = -b[node];
- lap_node = lap[node];
- for (j=0; j<n; j++) {
- if (j==node) {
- continue;
- }
- new_place_i += lap_node[j]*place[j];
- }
- desired_place[node] = new_place_i/(-lap_node[node]) - gap[node];
- }
+ /* not sure if this is better than 'place[node]!=target_place' */
+ if (fabs(place[node] - target_place) > 1e-9) {
+ break;
+ }
+ gap[node] = place[node] - cur_place;
+ }
+
+ /* compute desired place of block's reference point according
+ * to each node in the block:
+ */
+ for (i = left; i < right; i++) {
+ node = ordering[i];
+ new_place_i = -b[node];
+ lap_node = lap[node];
+ for (j = 0; j < n; j++) {
+ if (j == node) {
+ continue;
+ }
+ new_place_i += lap_node[j] * place[j];
+ }
+ desired_place[node] =
+ new_place_i / (-lap_node[node]) - gap[node];
+ }
+
+ /* reorder block by levels, and within levels
+ * by "relaxed" desired position
+ */
+ block_len = 0;
+ first_next_level = 0;
+ for (i = left; i < right; i = first_next_level) {
+ level = lev[ordering[i]];
+ if (level == num_levels) {
+ /* last_level */
+ first_next_level = right;
+ } else {
+ first_next_level = MIN(right, levels[level]);
+ }
- /* reorder block by levels, and within levels
- * by "relaxed" desired position
- */
- block_len = 0;
- first_next_level=0;
- for (i=left; i<right; i=first_next_level) {
- level = lev[ordering[i]];
- if (level==num_levels) {
- /* last_level */
- first_next_level=right;
- }
- else {
- first_next_level = MIN(right,levels[level]);
- }
-
- /* First, collect all nodes with desired places smaller
- * than current place
- */
- for (j=i; j<first_next_level; j++) {
- node = ordering[j];
- if (desired_place[node] < cur_place) {
- block[block_len++] = node;
- }
- }
- /* Second, collect all nodes with desired places equal
- * to current place
- */
- for (j=i; j<first_next_level; j++) {
- node = ordering[j];
- if (desired_place[node] == cur_place) {
- block[block_len++] = node;
- }
- }
- /* Third, collect all nodes with desired places greater
- * than current place
- */
- for (j=i; j<first_next_level; j++) {
- node = ordering[j];
- if (desired_place[node] > cur_place) {
- block[block_len++] = node;
- }
- }
+ /* First, collect all nodes with desired places smaller
+ * than current place
+ */
+ for (j = i; j < first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] < cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Second, collect all nodes with desired places equal
+ * to current place
+ */
+ for (j = i; j < first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] == cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Third, collect all nodes with desired places greater
+ * than current place
+ */
+ for (j = i; j < first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] > cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ }
+
+ /* loop through block and compute desired places of its prefixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i = 0; i < block_len; i++) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j = 0; j < i; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity *= 2;
+ /* update block stats */
+ des_place_block =
+ (block_deg * des_place_block +
+ (-lap_node[node]) * desired_place[node] +
+ toBlockConnectivity * cur_place) / (block_deg -
+ lap_node[node] +
+ toBlockConnectivity);
+ prefix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity - lap_node[node];
+ }
+
+ if (block_len == n) {
+ /* fix is needed since denominator was 0 in this case */
+ prefix_desired_place[n - 1] = cur_place; /* a "neutral" value */
+ }
+
+ /* loop through block and compute desired places of its suffixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i = block_len - 1; i >= 0; i--) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j = i + 1; j < block_len; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity *= 2;
+ /* update block stats */
+ des_place_block =
+ (block_deg * des_place_block +
+ (-lap_node[node]) * desired_place[node] +
+ toBlockConnectivity * cur_place) / (block_deg -
+ lap_node[node] +
+ toBlockConnectivity);
+ suffix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity - lap_node[node];
+ }
+
+ if (block_len == n) {
+ /* fix is needed since denominator was 0 in this case */
+ suffix_desired_place[0] = cur_place; /* a "neutral" value? */
+ }
+
+
+ /* now, find best place to split block */
+ best_i = -1;
+ max_movement = 0;
+ for (i = 0; i < block_len; i++) {
+ suffix_des_place = suffix_desired_place[i];
+ prefix_des_place =
+ i > 0 ? prefix_desired_place[i - 1] : suffix_des_place;
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place < prefix_des_place) {
+ if (suffix_des_place < cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
}
-
- /* loop through block and compute desired places of its prefixes */
- des_place_block = 0;
- block_deg = 0;
- for (i=0; i<block_len; i++) {
- node = block[i];
- toBlockConnectivity = 0;
- lap_node = lap[node];
- for (j=0; j<i; j++) {
- toBlockConnectivity -= lap_node[block[j]];
- }
- toBlockConnectivity*=2;
- /* update block stats */
- des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
- prefix_desired_place[i] = des_place_block;
- block_deg += toBlockConnectivity-lap_node[node];
+ suffix_des_place = prefix_des_place;
+ } else if (prefix_des_place > cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
+ movement =
+ (block_len - i) * fabs(suffix_des_place - cur_place) +
+ i * fabs(prefix_des_place - cur_place);
+ if (movement > max_movement) {
+ max_movement = movement;
+ best_i = i;
+ }
+ }
+ /* Actually move prefix and suffix */
+ if (best_i >= 0) {
+ suffix_des_place = suffix_desired_place[best_i];
+ prefix_des_place =
+ best_i >
+ 0 ? prefix_desired_place[best_i -
+ 1] : suffix_des_place;
+
+ /* compute right border of feasible move */
+ if (right >= n) {
+ /* no nodes after current block */
+ upper_bound = 1e9;
+ } else {
+ /* notice that we have to deduct 'gap[block[block_len-1]]'
+ * since all computations should be relative to
+ * the block's reference point
+ */
+ if (lev[ordering[right]] > lev[ordering[right - 1]]) {
+ upper_bound =
+ place[ordering[right]] - levels_gap -
+ gap[block[block_len - 1]];
+ } else {
+ upper_bound =
+ place[ordering[right]] -
+ gap[block[block_len - 1]];
+ }
+ }
+ suffix_des_place = MIN(suffix_des_place, upper_bound);
+ prefix_des_place = MAX(prefix_des_place, lower_bound);
+
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place < prefix_des_place) {
+ if (suffix_des_place < cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
}
+ suffix_des_place = prefix_des_place;
+ } else if (prefix_des_place > cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
- if (block_len==n) {
- /* fix is needed since denominator was 0 in this case */
- prefix_desired_place[n-1] = cur_place; /* a "neutral" value */
- }
+ /* move prefix: */
+ for (i = 0; i < best_i; i++) {
+ place[block[i]] = prefix_des_place + gap[block[i]];
+ }
+ /* move suffix: */
+ for (i = best_i; i < block_len; i++) {
+ place[block[i]] = suffix_des_place + gap[block[i]];
+ }
- /* loop through block and compute desired places of its suffixes */
- des_place_block = 0;
- block_deg = 0;
- for (i=block_len-1; i>=0; i--) {
- node = block[i];
- toBlockConnectivity = 0;
- lap_node = lap[node];
- for (j=i+1; j<block_len; j++) {
- toBlockConnectivity -= lap_node[block[j]];
- }
- toBlockConnectivity*=2;
- /* update block stats */
- des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
- suffix_desired_place[i] = des_place_block;
- block_deg += toBlockConnectivity-lap_node[node];
- }
-
- if (block_len==n) {
- /* fix is needed since denominator was 0 in this case */
- suffix_desired_place[0] = cur_place; /* a "neutral" value? */
- }
-
- /* now, find best place to split block */
- best_i = -1;
- max_movement = 0;
- for (i=0; i<block_len; i++) {
- suffix_des_place = suffix_desired_place[i];
- prefix_des_place = i>0 ? prefix_desired_place[i-1] : suffix_des_place;
- /* limit moves to ensure that the prefix is placed before the suffix */
- if (suffix_des_place<prefix_des_place) {
- if (suffix_des_place<cur_place) {
- if (prefix_des_place > cur_place) {
- prefix_des_place = cur_place;
- }
- suffix_des_place = prefix_des_place;
- }
- else if (prefix_des_place>cur_place) {
- prefix_des_place = suffix_des_place;
- }
- }
- movement = (block_len-i)*fabs(suffix_des_place-cur_place) + i*fabs(prefix_des_place-cur_place);
- if (movement>max_movement) {
- max_movement = movement;
- best_i = i;
- }
- }
- /* Actually move prefix and suffix */
- if (best_i>=0) {
- suffix_des_place = suffix_desired_place[best_i];
- prefix_des_place = best_i>0 ? prefix_desired_place[best_i-1] : suffix_des_place;
-
- /* compute right border of feasible move */
- if (right>=n) {
- /* no nodes after current block */
- upper_bound = 1e9;
- }
- else {
- /* notice that we have to deduct 'gap[block[block_len-1]]'
- * since all computations should be relative to
- * the block's reference point
- */
- if (lev[ordering[right]]>lev[ordering[right-1]]) {
- upper_bound = place[ordering[right]]-levels_gap-gap[block[block_len-1]];
- }
- else {
- upper_bound = place[ordering[right]]-gap[block[block_len-1]];
- }
- }
- suffix_des_place = MIN(suffix_des_place, upper_bound);
- prefix_des_place = MAX(prefix_des_place, lower_bound);
-
- /* limit moves to ensure that the prefix is placed before the suffix */
- if (suffix_des_place<prefix_des_place) {
- if (suffix_des_place<cur_place) {
- if (prefix_des_place > cur_place) {
- prefix_des_place = cur_place;
- }
- suffix_des_place = prefix_des_place;
- }
- else if (prefix_des_place>cur_place) {
- prefix_des_place = suffix_des_place;
- }
- }
-
- /* move prefix: */
- for (i=0; i<best_i; i++) {
- place[block[i]] = prefix_des_place+gap[block[i]];
- }
- /* move suffix: */
- for (i=best_i; i<block_len; i++) {
- place[block[i]] = suffix_des_place+gap[block[i]];
- }
-
-
- /* compute lower bound for next block */
- if (right<n && lev[ordering[right]]>lev[ordering[right-1]]) {
- lower_bound = place[block[block_len-1]]+levels_gap;
- }
- else {
- lower_bound = place[block[block_len-1]];
- }
-
-
- /* reorder 'ordering' to reflect change of places
- * Note that it is enough to reorder the level where
- * the split was done
- */
+ /* compute lower bound for next block */
+ if (right < n
+ && lev[ordering[right]] > lev[ordering[right - 1]]) {
+ lower_bound = place[block[block_len - 1]] + levels_gap;
+ } else {
+ lower_bound = place[block[block_len - 1]];
+ }
+
+
+ /* reorder 'ordering' to reflect change of places
+ * Note that it is enough to reorder the level where
+ * the split was done
+ */
#if 0
- int max_in_level, min_in_level;
-
- level = lev[best_i];
- if (level==num_levels) {
- /* last_level */
- max_in_level = MIN(right,n);
- }
- else {
- max_in_level = MIN(right,levels[level]);
- }
- if (level==0) {
- /* first level */
- min_in_level = MAX(left,0);
- }
- else {
- min_in_level = MAX(left,levels[level-1]);
- }
+ int max_in_level, min_in_level;
+
+ level = lev[best_i];
+ if (level == num_levels) {
+ /* last_level */
+ max_in_level = MIN(right, n);
+ } else {
+ max_in_level = MIN(right, levels[level]);
+ }
+ if (level == 0) {
+ /* first level */
+ min_in_level = MAX(left, 0);
+ } else {
+ min_in_level = MAX(left, levels[level - 1]);
+ }
#endif
- for (i=left; i<right; i++) {
- ordering[i] = block[i-left];
- }
- converged = converged && fabs(prefix_des_place-cur_place)<quad_prog_tol
- && fabs(suffix_des_place-cur_place)<quad_prog_tol;
-
-
- }
- else {
- /* no movement */
- /* compute lower bound for next block */
- if (right<n && lev[ordering[right]]>lev[ordering[right-1]]) {
- lower_bound = place[block[block_len-1]]+levels_gap;
- }
- else {
- lower_bound = place[block[block_len-1]];
- }
- }
- }
- orthog1f(n, place); /* for numerical stability, keep ||place|| small */
- computeHierarchyBoundaries(place, n, ordering, levels, num_levels, hierarchy_boundaries);
+ for (i = left; i < right; i++) {
+ ordering[i] = block[i - left];
+ }
+ converged = converged
+ && fabs(prefix_des_place - cur_place) < quad_prog_tol
+ && fabs(suffix_des_place - cur_place) < quad_prog_tol;
+
+
+ } else {
+ /* no movement */
+ /* compute lower bound for next block */
+ if (right < n
+ && lev[ordering[right]] > lev[ordering[right - 1]]) {
+ lower_bound = place[block[block_len - 1]] + levels_gap;
+ } else {
+ lower_bound = place[block[block_len - 1]];
+ }
+ }
}
-
- return counter;
+ orthog1f(n, place); /* for numerical stability, keep ||place|| small */
+ computeHierarchyBoundaries(place, n, ordering, levels, num_levels,
+ hierarchy_boundaries);
+ }
+
+ return counter;
}
-void
-deleteCMajEnv(CMajEnv *e)
+void deleteCMajEnv(CMajEnv * e)
{
- free (e->A[0]);
- free (e->A);
- free (e->lev);
- free (e->fArray1);
- free (e->fArray2);
- free (e->fArray3);
- free (e->fArray4);
- free (e->iArray1);
- free (e->iArray2);
- free (e->iArray3);
- free (e->iArray4);
- free (e);
+ free(e->A[0]);
+ free(e->A);
+ free(e->lev);
+ free(e->fArray1);
+ free(e->fArray2);
+ free(e->fArray3);
+ free(e->fArray4);
+ free(e->iArray1);
+ free(e->iArray2);
+ free(e->iArray3);
+ free(e->iArray4);
+ free(e);
}
-CMajEnv*
-initConstrainedMajorization(float* packedMat, int n, int* ordering,
- int* levels, int num_levels)
+CMajEnv *initConstrainedMajorization(float *packedMat, int n,
+ int *ordering, int *levels,
+ int num_levels)
{
- int i, level=-1, start_of_level_above=0;
- CMajEnv *e = GNEW(CMajEnv);
- e->A=NULL;
- e->n=n;
- e->ordering=ordering;
- e->levels=levels;
- e->num_levels=num_levels;
- e->A = unpackMatrix(packedMat,n);
- e->lev = N_GNEW(n, int);
- for (i=0; i<e->n; i++) {
- if (i>=start_of_level_above) {
- level++;
- start_of_level_above=(level==num_levels)?e->n:levels[level];
- }
- e->lev[ordering[i]]=level;
- }
- e->fArray1 = N_GNEW(n,float);
- e->fArray2 = N_GNEW(n,float);
- e->fArray3 = N_GNEW(n,float);
- e->fArray4 = N_GNEW(n,float);
- e->iArray1 = N_GNEW(n,int);
- e->iArray2 = N_GNEW(n,int);
- e->iArray3 = N_GNEW(n,int);
- e->iArray4 = N_GNEW(n,int);
- return e;
+ int i, level = -1, start_of_level_above = 0;
+ CMajEnv *e = GNEW(CMajEnv);
+ e->A = NULL;
+ e->n = n;
+ e->ordering = ordering;
+ e->levels = levels;
+ e->num_levels = num_levels;
+ e->A = unpackMatrix(packedMat, n);
+ e->lev = N_GNEW(n, int);
+ for (i = 0; i < e->n; i++) {
+ if (i >= start_of_level_above) {
+ level++;
+ start_of_level_above =
+ (level == num_levels) ? e->n : levels[level];
+ }
+ e->lev[ordering[i]] = level;
+ }
+ e->fArray1 = N_GNEW(n, float);
+ e->fArray2 = N_GNEW(n, float);
+ e->fArray3 = N_GNEW(n, float);
+ e->fArray4 = N_GNEW(n, float);
+ e->iArray1 = N_GNEW(n, int);
+ e->iArray2 = N_GNEW(n, int);
+ e->iArray3 = N_GNEW(n, int);
+ e->iArray4 = N_GNEW(n, int);
+ return e;
}
-#endif /* DIGCOLA */
-
+#endif /* DIGCOLA */
extern int constrained_majorization_new_with_gaps(CMajEnv*, float*, float**,
int, int, int, float*, float);
+#ifdef IPSEPCOLA
extern int constrained_majorization_gradient_projection(CMajEnv *e,
float * b, float ** coords, int ndims, int cur_axis, int max_iterations,
float * hierarchy_boundaries,float levels_gap);
-
-extern int constrained_majorization_new_with_gaps(CMajEnv*, float*, float**,
- int, int, int, float*, float);
+#endif
extern void deleteCMajEnv(CMajEnv *e);
extern float** unpackMatrix(float * packedMat, int n);
+/* vim:set shiftwidth=4 ts=8: */
+
+/**
+ *
+ * Authors:
+ * Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * This version is released under the CPL (Common Public License) with
+ * the Graphviz distribution.
+ * A version is also available under the LGPL as part of the Adaptagrams
+ * project: http://sourceforge.net/projects/adaptagrams.
+ * If you make improvements or bug fixes to this code it would be much
+ * appreciated if you could also contribute those changes back to the
+ * Adaptagrams repository.
+ */
+
/**********************************************************
*
* Solve a quadratic function f(X) = X' e->A X + b X
**********************************************************/
#include "digcola.h"
-#ifdef DIGCOLA
+#ifdef IPSEPCOLA
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "quad_prog_vpsc.h"
#include "quad_prog_solver.h"
-//#define CONMAJ_LOGGING 1
+/* #define CONMAJ_LOGGING 1 */
#define quad_prog_tol 1e-4
/*
* Use gradient-projection to solve Variable Placement with Separation Constraints problem.
*/
-int
-constrained_majorization_vpsc(CMajEnvVPSC *e, float * b, float *place,
- int max_iterations)
+int
+constrained_majorization_vpsc(CMajEnvVPSC * e, float *b, float *place,
+ int max_iterations)
{
- int i,j,counter;
- int n=e->nv+e->nldv; // for laplacian computation need number of real vars and those
- // dummy vars included in lap
- if(max_iterations==0) return 0;
-
- bool converged=false;
+ int i, j, counter;
+ float *g, *old_place, *d;
+ /* for laplacian computation need number of real vars and those
+ * dummy vars included in lap
+ */
+ int n = e->nv + e->nldv;
+ bool converged = false;
#ifdef CONMAJ_LOGGING
- static int call_no=0;
-#endif // CONMAJ_LOGGING
- float *g = e->fArray1;
- float *old_place = e->fArray2;
- float *d = e->fArray3;
- //fprintf(stderr,"Entered: constrained_majorization_vpsc, #constraints=%d\n",e->m);
- if(e->m>0) {
- for (i=0;i<n;i++) {
- setVariableDesiredPos(e->vs[i],place[i]);
- }
- //fprintf(stderr," calling satisfyVPSC...\n");
- satisfyVPSC(e->vpsc);
- for (i=0;i<n;i++) {
- place[i]=getVariablePos(e->vs[i]);
- //fprintf(stderr,"vs[%d]=%f\n",i,place[i]);
- }
- //fprintf(stderr," done.\n");
+ static int call_no = 0;
+#endif /* CONMAJ_LOGGING */
+
+ if (max_iterations == 0)
+ return 0;
+ g = e->fArray1;
+ old_place = e->fArray2;
+ d = e->fArray3;
+ /* fprintf(stderr,"Entered: constrained_majorization_vpsc, #constraints=%d\n",e->m); */
+ if (e->m > 0) {
+ for (i = 0; i < n; i++) {
+ setVariableDesiredPos(e->vs[i], place[i]);
+ }
+ /* fprintf(stderr," calling satisfyVPSC...\n"); */
+ satisfyVPSC(e->vpsc);
+ for (i = 0; i < n; i++) {
+ place[i] = getVariablePos(e->vs[i]);
+ /* fprintf(stderr,"vs[%d]=%f\n",i,place[i]); */
+ }
+ /* fprintf(stderr," done.\n"); */
}
#ifdef CONMAJ_LOGGING
- float prev_stress=0;
- for (i=0; i<n; i++) {
- prev_stress += 2*b[i]*place[i];
- for (j=0; j<n; j++) {
- prev_stress -= e->A[i][j]*place[j]*place[i];
- }
+ float prev_stress = 0;
+ for (i = 0; i < n; i++) {
+ prev_stress += 2 * b[i] * place[i];
+ for (j = 0; j < n; j++) {
+ prev_stress -= e->A[i][j] * place[j] * place[i];
}
- FILE *logfile = fopen("constrained_majorization_log","a");
+ }
+ FILE *logfile = fopen("constrained_majorization_log", "a");
- //fprintf(logfile,"grad proj %d: stress=%f\n",call_no,prev_stress);
+ /* fprintf(logfile,"grad proj %d: stress=%f\n",call_no,prev_stress); */
#endif
-
- for (counter=0; counter<max_iterations&&!converged; counter++) {
- //fprintf(stderr,".");
- converged=true;
- // find steepest descent direction
- for (i=0; i<n; i++) {
- old_place[i]=place[i];
- g[i] = 2*b[i];
- for (j=0; j<n; j++) {
- g[i] -= 2*e->A[i][j]*place[j];
- }
- }
- float numerator = 0, denominator = 0, r;
- for (i=0; i<n; i++) {
- numerator += g[i]*g[i];
- r=0;
- for (j=0; j<n; j++) {
- r += 2*e->A[i][j]*g[j];
- }
- denominator -= r*g[i];
- }
- float alpha = numerator/denominator;
- for (i=0; i<n; i++) {
- place[i]-=alpha*g[i];
- }
- if(e->m>0) {
- //project to constraint boundary
- for (i=0;i<n;i++) {
- setVariableDesiredPos(e->vs[i],place[i]);
- }
- satisfyVPSC(e->vpsc);
- for (i=0;i<n;i++) {
- place[i]=getVariablePos(e->vs[i]);
- }
- }
- // set place to the intersection of old_place-g and boundary and compute d, the vector from intersection pnt to projection pnt
- for (i=0; i<n; i++) {
- d[i]=place[i]-old_place[i];
- }
- // now compute beta
- numerator = 0, denominator = 0;
- for (i=0; i<n; i++) {
- numerator += g[i]*d[i];
- r=0;
- for (j=0; j<n; j++) {
- r += 2*e->A[i][j]*d[j];
- }
- denominator += r*d[i];
- }
- float beta = numerator/denominator;
- float test=0;
- for (i=0; i<n; i++) {
- // beta > 1.0 takes us back outside the feasible region
- // beta < 0 clearly not useful and may happen due to numerical imp.
- if(beta>0&&beta<1.0) {
- place[i]=old_place[i]+beta*d[i];
- }
- test+= fabs(place[i]-old_place[i]);
- }
+ for (counter = 0; counter < max_iterations && !converged; counter++) {
+ float test = 0;
+ float alpha, beta;
+ float numerator = 0, denominator = 0, r;
+ /* fprintf(stderr,"."); */
+ converged = true;
+ /* find steepest descent direction */
+ for (i = 0; i < n; i++) {
+ old_place[i] = place[i];
+ g[i] = 2 * b[i];
+ for (j = 0; j < n; j++) {
+ g[i] -= 2 * e->A[i][j] * place[j];
+ }
+ }
+ for (i = 0; i < n; i++) {
+ numerator += g[i] * g[i];
+ r = 0;
+ for (j = 0; j < n; j++) {
+ r += 2 * e->A[i][j] * g[j];
+ }
+ denominator -= r * g[i];
+ }
+ alpha = numerator / denominator;
+ for (i = 0; i < n; i++) {
+ place[i] -= alpha * g[i];
+ }
+ if (e->m > 0) {
+ /* project to constraint boundary */
+ for (i = 0; i < n; i++) {
+ setVariableDesiredPos(e->vs[i], place[i]);
+ }
+ satisfyVPSC(e->vpsc);
+ for (i = 0; i < n; i++) {
+ place[i] = getVariablePos(e->vs[i]);
+ }
+ }
+ /* set place to the intersection of old_place-g and boundary and
+ * compute d, the vector from intersection pnt to projection pnt
+ */
+ for (i = 0; i < n; i++) {
+ d[i] = place[i] - old_place[i];
+ }
+ /* now compute beta */
+ numerator = 0, denominator = 0;
+ for (i = 0; i < n; i++) {
+ numerator += g[i] * d[i];
+ r = 0;
+ for (j = 0; j < n; j++) {
+ r += 2 * e->A[i][j] * d[j];
+ }
+ denominator += r * d[i];
+ }
+ beta = numerator / denominator;
+
+ for (i = 0; i < n; i++) {
+ /* beta > 1.0 takes us back outside the feasible region
+ * beta < 0 clearly not useful and may happen due to numerical imp.
+ */
+ if (beta > 0 && beta < 1.0) {
+ place[i] = old_place[i] + beta * d[i];
+ }
+ test += fabs(place[i] - old_place[i]);
+ }
#ifdef CONMAJ_LOGGING
- float stress=0;
- for (i=0; i<n; i++) {
- stress += 2*b[i]*place[i];
- for (j=0; j<n; j++) {
- stress -= e->A[i][j]*place[j]*place[i];
- }
- }
- fprintf(logfile,"%d: stress=%f, test=%f, %s\n",call_no,stress,test,(stress>=prev_stress)?"No Improvement":"");
- prev_stress=stress;
+ float stress = 0;
+ for (i = 0; i < n; i++) {
+ stress += 2 * b[i] * place[i];
+ for (j = 0; j < n; j++) {
+ stress -= e->A[i][j] * place[j] * place[i];
+ }
+ }
+ fprintf(logfile, "%d: stress=%f, test=%f, %s\n", call_no, stress,
+ test, (stress >= prev_stress) ? "No Improvement" : "");
+ prev_stress = stress;
#endif
- if(test>quad_prog_tol) {
- converged=false;
- }
+ if (test > quad_prog_tol) {
+ converged = false;
}
+ }
#ifdef CONMAJ_LOGGING
- call_no++;
- fclose(logfile);
+ call_no++;
+ fclose(logfile);
#endif
- return counter;
+ return counter;
}
/*
* 1=one separation constraint for each edge (in acyclic subgraph)
* 2=DiG-CoLa level constraints
*/
-CMajEnvVPSC*
-initCMajVPSC(int n, float* packedMat, vtx_data* graph,
- ipsep_options *opt, int diredges)
+CMajEnvVPSC *initCMajVPSC(int n, float *packedMat, vtx_data * graph,
+ ipsep_options * opt, int diredges)
{
- int i,j;
- // nv is the number of real nodes
+ int i, j;
+ /* nv is the number of real nodes */
int nConCs;
- //fprintf(stderr,"Entered initCMajVPSC\n");
- CMajEnvVPSC *e = GNEW(CMajEnvVPSC);
- e->A=NULL;
- e->packedMat=packedMat;
+ /* fprintf(stderr,"Entered initCMajVPSC\n"); */
+ CMajEnvVPSC *e = GNEW(CMajEnvVPSC);
+ e->A = NULL;
+ e->packedMat = packedMat;
/* if we have clusters then we'll need two constraints for each var in
* a cluster */
- e->nldv = 2*opt->clusters->nclusters;
+ e->nldv = 2 * opt->clusters->nclusters;
e->nv = n - e->nldv;
e->ndv = 0;
- e->gcs=NULL;
- e->vs=N_GNEW(n, Variable*);
- for(i=0;i<n;i++) {
- e->vs[i]=newVariable(i,1.0,1.0);
+ e->gcs = NULL;
+ e->vs = N_GNEW(n, Variable *);
+ for (i = 0; i < n; i++) {
+ e->vs[i] = newVariable(i, 1.0, 1.0);
}
- e->gm=0;
- if(diredges==1) {
- if(Verbose)
- fprintf(stderr," generate edge constraints...\n");
- for(i=0;i<e->nv;i++) {
- for(j=1;j<graph[i].nedges;j++) {
- //fprintf(stderr,"edist=%f\n",graph[i].edists[j]);
- if(graph[i].edists[j]>0.01) {
- e->gm++;
- }
- }
- }
- e->gcs=newConstraints(e->gm);
- e->gm=0;
- for(i=0;i<e->nv;i++) {
- for(j=1;j<graph[i].nedges;j++) {
- int u=i,v=graph[i].edges[j];
- if(graph[i].edists[j]>0) {
- e->gcs[e->gm++]=newConstraint(e->vs[u],e->vs[v],opt->edge_gap);
- }
- }
- }
- } else if(diredges==2) {
- int *ordering=NULL, *ls=NULL, cvar;
- DigColaLevel* levels;
- Variable** vs = e->vs;
- // e->ndv is the number of dummy variables required, one for each boundary
- compute_hierarchy(graph, n, 1e-2, 1e-1, NULL, &ordering, &ls, &e->ndv);
- levels=assign_digcola_levels(ordering,n,ls,e->ndv);
- if(Verbose)
- fprintf(stderr,"Found %d DiG-CoLa boundaries\n",e->ndv);
- e->gm=get_num_digcola_constraints(levels,e->ndv+1)+e->ndv-1;
- e->gcs=newConstraints(e->gm);
- e->gm=0;
- e->vs = N_GNEW(n+e->ndv,Variable*);
- for(i=0;i<n;i++) {
- e->vs[i]=vs[i];
- }
- free(vs);
- // create dummy vars
- for(i=0;i<e->ndv;i++) {
- // dummy vars should have 0 weight
- cvar=n+i;
- e->vs[cvar]=newVariable(cvar,1.0,0.000001);
- }
- double halfgap=opt->edge_gap;
- for(i=0;i<e->ndv;i++) {
- cvar=n+i;
- // outgoing constraints for each var in level below boundary
- for(j=0;j<levels[i].num_nodes;j++) {
- e->gcs[e->gm++]=newConstraint(e->vs[levels[i].nodes[j]],e->vs[cvar],halfgap);
- }
- // incoming constraints for each var in level above boundary
- for(j=0;j<levels[i+1].num_nodes;j++) {
- e->gcs[e->gm++]=newConstraint(e->vs[cvar],e->vs[levels[i+1].nodes[j]],halfgap);
- }
- }
- // constraints between adjacent boundary dummy vars
- for(i=0;i<e->ndv-1;i++) {
- e->gcs[e->gm++]=newConstraint(e->vs[n+i],e->vs[n+i+1],0);
- }
+ e->gm = 0;
+ if (diredges == 1) {
+ if (Verbose)
+ fprintf(stderr, " generate edge constraints...\n");
+ for (i = 0; i < e->nv; i++) {
+ for (j = 1; j < graph[i].nedges; j++) {
+ /* fprintf(stderr,"edist=%f\n",graph[i].edists[j]); */
+ if (graph[i].edists[j] > 0.01) {
+ e->gm++;
+ }
+ }
+ }
+ e->gcs = newConstraints(e->gm);
+ e->gm = 0;
+ for (i = 0; i < e->nv; i++) {
+ for (j = 1; j < graph[i].nedges; j++) {
+ int u = i, v = graph[i].edges[j];
+ if (graph[i].edists[j] > 0) {
+ e->gcs[e->gm++] =
+ newConstraint(e->vs[u], e->vs[v], opt->edge_gap);
+ }
+ }
+ }
+ } else if (diredges == 2) {
+ int *ordering = NULL, *ls = NULL, cvar;
+ double halfgap;
+ DigColaLevel *levels;
+ Variable **vs = e->vs;
+ /* e->ndv is the number of dummy variables required, one for each boundary */
+ compute_hierarchy(graph, n, 1e-2, 1e-1, NULL, &ordering, &ls,
+ &e->ndv);
+ levels = assign_digcola_levels(ordering, n, ls, e->ndv);
+ if (Verbose)
+ fprintf(stderr, "Found %d DiG-CoLa boundaries\n", e->ndv);
+ e->gm =
+ get_num_digcola_constraints(levels, e->ndv + 1) + e->ndv - 1;
+ e->gcs = newConstraints(e->gm);
+ e->gm = 0;
+ e->vs = N_GNEW(n + e->ndv, Variable *);
+ for (i = 0; i < n; i++) {
+ e->vs[i] = vs[i];
+ }
+ free(vs);
+ /* create dummy vars */
+ for (i = 0; i < e->ndv; i++) {
+ /* dummy vars should have 0 weight */
+ cvar = n + i;
+ e->vs[cvar] = newVariable(cvar, 1.0, 0.000001);
+ }
+ halfgap = opt->edge_gap;
+ for (i = 0; i < e->ndv; i++) {
+ cvar = n + i;
+ /* outgoing constraints for each var in level below boundary */
+ for (j = 0; j < levels[i].num_nodes; j++) {
+ e->gcs[e->gm++] =
+ newConstraint(e->vs[levels[i].nodes[j]], e->vs[cvar],
+ halfgap);
+ }
+ /* incoming constraints for each var in level above boundary */
+ for (j = 0; j < levels[i + 1].num_nodes; j++) {
+ e->gcs[e->gm++] =
+ newConstraint(e->vs[cvar],
+ e->vs[levels[i + 1].nodes[j]], halfgap);
+ }
+ }
+ /* constraints between adjacent boundary dummy vars */
+ for (i = 0; i < e->ndv - 1; i++) {
+ e->gcs[e->gm++] =
+ newConstraint(e->vs[n + i], e->vs[n + i + 1], 0);
+ }
}
- //fprintf(stderr," generate edge constraints... done: n=%d,m=%d\n",e->n,e->gm);
- if(opt->clusters->nclusters>0) {
- //fprintf(stderr," generate cluster containment constraints...\n");
- Constraint** ecs=e->gcs;
- nConCs=2*opt->clusters->nvars;
- e->gcs=newConstraints(e->gm+nConCs);
- for(i=0;i<e->gm;i++) {
- e->gcs[i]=ecs[i];
- }
- if(ecs!=NULL) deleteConstraints(0,ecs);
- for(i=0;i<opt->clusters->nclusters;i++) {
- for(j=0;j<opt->clusters->clustersizes[i];j++) {
- Variable* v=e->vs[opt->clusters->clusters[i][j]];
- Variable* cl=e->vs[e->nv+2*i];
- Variable* cr=e->vs[e->nv+2*i+1];
- e->gcs[e->gm++]=newConstraint(cl,v,0);
- e->gcs[e->gm++]=newConstraint(v,cr,0);
- }
- }
- //fprintf(stderr," containment constraints... done: \n");
+ /* fprintf(stderr," generate edge constraints... done: n=%d,m=%d\n",e->n,e->gm); */
+ if (opt->clusters->nclusters > 0) {
+ /* fprintf(stderr," generate cluster containment constraints...\n"); */
+ Constraint **ecs = e->gcs;
+ nConCs = 2 * opt->clusters->nvars;
+ e->gcs = newConstraints(e->gm + nConCs);
+ for (i = 0; i < e->gm; i++) {
+ e->gcs[i] = ecs[i];
+ }
+ if (ecs != NULL)
+ deleteConstraints(0, ecs);
+ for (i = 0; i < opt->clusters->nclusters; i++) {
+ for (j = 0; j < opt->clusters->clustersizes[i]; j++) {
+ Variable *v = e->vs[opt->clusters->clusters[i][j]];
+ Variable *cl = e->vs[e->nv + 2 * i];
+ Variable *cr = e->vs[e->nv + 2 * i + 1];
+ e->gcs[e->gm++] = newConstraint(cl, v, 0);
+ e->gcs[e->gm++] = newConstraint(v, cr, 0);
+ }
+ }
+ /* fprintf(stderr," containment constraints... done: \n"); */
}
- e->m=0;
- e->cs=NULL;
- if(e->gm>0) {
- e->vpsc = newIncVPSC(n+e->ndv,e->vs,e->gm,e->gcs);
- e->m=e->gm;
- e->cs=e->gcs;
+ e->m = 0;
+ e->cs = NULL;
+ if (e->gm > 0) {
+ e->vpsc = newIncVPSC(n + e->ndv, e->vs, e->gm, e->gcs);
+ e->m = e->gm;
+ e->cs = e->gcs;
}
- if(packedMat!=NULL) {
- e->A = unpackMatrix(packedMat,n);
+ if (packedMat != NULL) {
+ e->A = unpackMatrix(packedMat, n);
}
#ifdef MOSEK
- e->mosekEnv=NULL;
- if(opt->mosek) {
- e->mosekEnv = mosek_init_sep(e->packedMat,n,e->ndv,e->gcs,e->gm);
+ e->mosekEnv = NULL;
+ if (opt->mosek) {
+ e->mosekEnv =
+ mosek_init_sep(e->packedMat, n, e->ndv, e->gcs, e->gm);
}
#endif
- e->fArray1 = N_GNEW(n,float);
- e->fArray2 = N_GNEW(n,float);
- e->fArray3 = N_GNEW(n,float);
- if(Verbose)
- fprintf(stderr," initCMajVPSC done: %d global constraints generated.\n",e->m);
- return e;
+ e->fArray1 = N_GNEW(n, float);
+ e->fArray2 = N_GNEW(n, float);
+ e->fArray3 = N_GNEW(n, float);
+ if (Verbose)
+ fprintf(stderr,
+ " initCMajVPSC done: %d global constraints generated.\n",
+ e->m);
+ return e;
}
-void
-deleteCMajEnvVPSC(CMajEnvVPSC *e)
+
+void deleteCMajEnvVPSC(CMajEnvVPSC * e)
{
int i;
- if(e->A!=NULL) {
- free (e->A[0]);
- free (e->A);
+ if (e->A != NULL) {
+ free(e->A[0]);
+ free(e->A);
}
- if(e->m>0) {
- deleteVPSC(e->vpsc);
- if(e->cs!=e->gcs && e->gcs!=NULL) deleteConstraints(0,e->gcs);
- deleteConstraints(e->m,e->cs);
- for(i=0;i<e->nv+e->nldv+e->ndv;i++) {
- deleteVariable(e->vs[i]);
- }
- free (e->vs);
+ if (e->m > 0) {
+ deleteVPSC(e->vpsc);
+ if (e->cs != e->gcs && e->gcs != NULL)
+ deleteConstraints(0, e->gcs);
+ deleteConstraints(e->m, e->cs);
+ for (i = 0; i < e->nv + e->nldv + e->ndv; i++) {
+ deleteVariable(e->vs[i]);
+ }
+ free(e->vs);
}
- free (e->fArray1);
- free (e->fArray2);
- free (e->fArray3);
+ free(e->fArray1);
+ free(e->fArray2);
+ free(e->fArray3);
#ifdef MOSEK
- if(e->mosekEnv) {
- mosek_delete(e->mosekEnv);
+ if (e->mosekEnv) {
+ mosek_delete(e->mosekEnv);
}
-#endif //MOSEK
- free (e);
+#endif /* MOSEK */
+ free(e);
}
-// generate non-overlap constraints inside each cluster, including dummy
-// nodes at bounds of cluster
-// generate constraints again for top level nodes and clusters treating
-// clusters as rectangles of dim (l,r,b,t)
-// for each cluster map in-constraints to l out-constraints to r
-//
-// For now, we'll keep the global containment constraints already
-// generated for each cluster, and simply generate non-overlap constraints
-// for all nodes and then an additional set of non-overlap constraints for
-// clusters that we'll map back to the dummy vars as above.
-void generateNonoverlapConstraints(
- CMajEnvVPSC* e,
- float nsizeScale,
- float** coords,
- int k,
- bool transitiveClosure,
- ipsep_options* opt
-) {
+/* generate non-overlap constraints inside each cluster, including dummy
+ * nodes at bounds of cluster
+ * generate constraints again for top level nodes and clusters treating
+ * clusters as rectangles of dim (l,r,b,t)
+ * for each cluster map in-constraints to l out-constraints to r
+ *
+ * For now, we'll keep the global containment constraints already
+ * generated for each cluster, and simply generate non-overlap constraints
+ * for all nodes and then an additional set of non-overlap constraints for
+ * clusters that we'll map back to the dummy vars as above.
+ */
+void generateNonoverlapConstraints(CMajEnvVPSC * e,
+ float nsizeScale,
+ float **coords,
+ int k,
+ bool transitiveClosure,
+ ipsep_options * opt)
+{
Constraint **csol, **csolptr;
- int i,j,mol=0;
- int n=e->nv+e->nldv;
+ int i, j, mol = 0;
+ int n = e->nv + e->nldv;
boxf bb[n];
- bool genclusters=opt->clusters->nclusters>0;
- if(genclusters) {
- // n is the number of real variables, not dummy cluster vars
- n-=2*opt->clusters->nclusters;
+ bool genclusters = opt->clusters->nclusters > 0;
+ if (genclusters) {
+ /* n is the number of real variables, not dummy cluster vars */
+ n -= 2 * opt->clusters->nclusters;
}
- if(k==0) {
- // grow a bit in the x dimension, so that if overlap is resolved
- // horizontally then it won't be considered overlapping vertically
- nsizeScale*=1.0001;
+ if (k == 0) {
+ /* grow a bit in the x dimension, so that if overlap is resolved
+ * horizontally then it won't be considered overlapping vertically
+ */
+ nsizeScale *= 1.0001;
}
- for(i=0;i<n;i++) {
- bb[i].LL.x=coords[0][i]-nsizeScale*opt->nsize[i].x/2.0-opt->gap.x/2.0;
- bb[i].UR.x=coords[0][i]+nsizeScale*opt->nsize[i].x/2.0+opt->gap.x/2.0;
- bb[i].LL.y=coords[1][i]-nsizeScale*opt->nsize[i].y/2.0-opt->gap.y/2.0;
- bb[i].UR.y=coords[1][i]+nsizeScale*opt->nsize[i].y/2.0+opt->gap.y/2.0;
+ for (i = 0; i < n; i++) {
+ bb[i].LL.x =
+ coords[0][i] - nsizeScale * opt->nsize[i].x / 2.0 -
+ opt->gap.x / 2.0;
+ bb[i].UR.x =
+ coords[0][i] + nsizeScale * opt->nsize[i].x / 2.0 +
+ opt->gap.x / 2.0;
+ bb[i].LL.y =
+ coords[1][i] - nsizeScale * opt->nsize[i].y / 2.0 -
+ opt->gap.y / 2.0;
+ bb[i].UR.y =
+ coords[1][i] + nsizeScale * opt->nsize[i].y / 2.0 +
+ opt->gap.y / 2.0;
}
- if(genclusters) {
- Constraint **cscl[opt->clusters->nclusters+1];
- int cm[opt->clusters->nclusters+1];
- for(i=0;i<opt->clusters->nclusters;i++) {
- int cn=opt->clusters->clustersizes[i];
- Variable* cvs[cn+2];
- boxf cbb[cn+2];
- // compute cluster bounding bb
- boxf container;
- container.LL.x=container.LL.y=DBL_MAX;
- container.UR.x=container.UR.y=-DBL_MAX;
- for(j=0;j<cn;j++) {
- int iv=opt->clusters->clusters[i][j];
- cvs[j]=e->vs[iv];
- B2BF(bb[iv],cbb[j]);
- EXPANDBB(container,bb[iv]);
- }
- B2BF(container,opt->clusters->bb[i]);
- cvs[cn]=e->vs[n+2*i];
- cvs[cn+1]=e->vs[n+2*i+1];
- B2BF(container,cbb[cn]);
- B2BF(container,cbb[cn+1]);
- if(k==0) {
- cbb[cn].UR.x=container.LL.x+0.0001;
- cbb[cn+1].LL.x=container.UR.x-0.0001;
- cm[i] = genXConstraints(cn+2,cbb,cvs,&cscl[i],transitiveClosure);
- } else {
- cbb[cn].UR.y=container.LL.y+0.0001;
- cbb[cn+1].LL.y=container.UR.y-0.0001;
- cm[i] = genYConstraints(cn+2,cbb,cvs,&cscl[i]);
- }
- mol+=cm[i];
- }
- // generate top level constraints
- {
- int cn=opt->clusters->ntoplevel+opt->clusters->nclusters;
- Variable* cvs[cn];
- boxf cbb[cn];
- for(i=0;i<opt->clusters->ntoplevel;i++) {
- int iv=opt->clusters->toplevel[i];
- cvs[i]=e->vs[iv];
- B2BF(bb[iv],cbb[i]);
- }
- // make dummy variables for clusters
- for(i=opt->clusters->ntoplevel;i<cn;i++) {
- cvs[i]=newVariable(123+i,1,1);
- j=i-opt->clusters->ntoplevel;
- B2BF(opt->clusters->bb[j],cbb[i]);
- }
- i=opt->clusters->nclusters;
- if(k==0) {
- cm[i] = genXConstraints(cn,cbb,cvs,&cscl[i],transitiveClosure);
- } else {
- cm[i] = genYConstraints(cn,cbb,cvs,&cscl[i]);
- }
- // remap constraints from tmp dummy vars to cluster l and r vars
- for(i=opt->clusters->ntoplevel;i<cn;i++) {
- j=i-opt->clusters->ntoplevel;
- // dgap is the change in required constraint gap.
- // since we are going from a source rectangle the size
- // of the cluster bounding box to a zero width (in x dim,
- // zero height in y dim) rectangle, the change will be
- // half the bb width.
- double dgap;
- if(k==0) {
- dgap=-(cbb[i].UR.x-cbb[i].LL.x)/2.0;
- } else {
- dgap=-(cbb[i].UR.y-cbb[i].LL.y)/2.0;
- }
- remapInConstraints(cvs[i],e->vs[n+2*j],dgap);
- remapOutConstraints(cvs[i],e->vs[n+2*j+1],dgap);
- // there may be problems with cycles between
- // cluster non-overlap and diredge constraints,
- // to resolve:
- //
- // for each constraint c:v->cvs[i]:
- // if exists diredge constraint u->v where u in c:
- // remap v->cl to cr->v (gap = height(v)/2)
- //
- // in=getInConstraints(cvs[i])
- // for(c : in) {
- // assert(c.right==cvs[i]);
- // vin=getOutConstraints(v=c.left)
- // for(d : vin) {
- // if(d.left.cluster==i):
- // tmp=d.left
- // d.left=d.right
- // d.right=tmp
- // d.gap==height(d.right)/2
- // }
- // }
- //
- deleteVariable(cvs[i]);
- }
- mol+=cm[opt->clusters->nclusters];
- }
- csolptr=csol=newConstraints(mol);
- for(i=0;i<opt->clusters->nclusters+1;i++) {
- // copy constraints into csol
- for(j=0;j<cm[i];j++) {
- *csolptr++=cscl[i][j];
- }
- deleteConstraints(0,cscl[i]);
- }
+ if (genclusters) {
+ Constraint **cscl[opt->clusters->nclusters + 1];
+ int cm[opt->clusters->nclusters + 1];
+ for (i = 0; i < opt->clusters->nclusters; i++) {
+ int cn = opt->clusters->clustersizes[i];
+ Variable *cvs[cn + 2];
+ boxf cbb[cn + 2];
+ /* compute cluster bounding bb */
+ boxf container;
+ container.LL.x = container.LL.y = DBL_MAX;
+ container.UR.x = container.UR.y = -DBL_MAX;
+ for (j = 0; j < cn; j++) {
+ int iv = opt->clusters->clusters[i][j];
+ cvs[j] = e->vs[iv];
+ B2BF(bb[iv], cbb[j]);
+ EXPANDBB(container, bb[iv]);
+ }
+ B2BF(container, opt->clusters->bb[i]);
+ cvs[cn] = e->vs[n + 2 * i];
+ cvs[cn + 1] = e->vs[n + 2 * i + 1];
+ B2BF(container, cbb[cn]);
+ B2BF(container, cbb[cn + 1]);
+ if (k == 0) {
+ cbb[cn].UR.x = container.LL.x + 0.0001;
+ cbb[cn + 1].LL.x = container.UR.x - 0.0001;
+ cm[i] =
+ genXConstraints(cn + 2, cbb, cvs, &cscl[i],
+ transitiveClosure);
+ } else {
+ cbb[cn].UR.y = container.LL.y + 0.0001;
+ cbb[cn + 1].LL.y = container.UR.y - 0.0001;
+ cm[i] = genYConstraints(cn + 2, cbb, cvs, &cscl[i]);
+ }
+ mol += cm[i];
+ }
+ /* generate top level constraints */
+ {
+ int cn = opt->clusters->ntoplevel + opt->clusters->nclusters;
+ Variable *cvs[cn];
+ boxf cbb[cn];
+ for (i = 0; i < opt->clusters->ntoplevel; i++) {
+ int iv = opt->clusters->toplevel[i];
+ cvs[i] = e->vs[iv];
+ B2BF(bb[iv], cbb[i]);
+ }
+ /* make dummy variables for clusters */
+ for (i = opt->clusters->ntoplevel; i < cn; i++) {
+ cvs[i] = newVariable(123 + i, 1, 1);
+ j = i - opt->clusters->ntoplevel;
+ B2BF(opt->clusters->bb[j], cbb[i]);
+ }
+ i = opt->clusters->nclusters;
+ if (k == 0) {
+ cm[i] =
+ genXConstraints(cn, cbb, cvs, &cscl[i],
+ transitiveClosure);
+ } else {
+ cm[i] = genYConstraints(cn, cbb, cvs, &cscl[i]);
+ }
+ /* remap constraints from tmp dummy vars to cluster l and r vars */
+ for (i = opt->clusters->ntoplevel; i < cn; i++) {
+ double dgap;
+ j = i - opt->clusters->ntoplevel;
+ /* dgap is the change in required constraint gap.
+ * since we are going from a source rectangle the size
+ * of the cluster bounding box to a zero width (in x dim,
+ * zero height in y dim) rectangle, the change will be
+ * half the bb width.
+ */
+ if (k == 0) {
+ dgap = -(cbb[i].UR.x - cbb[i].LL.x) / 2.0;
+ } else {
+ dgap = -(cbb[i].UR.y - cbb[i].LL.y) / 2.0;
+ }
+ remapInConstraints(cvs[i], e->vs[n + 2 * j], dgap);
+ remapOutConstraints(cvs[i], e->vs[n + 2 * j + 1], dgap);
+ /* there may be problems with cycles between
+ * cluster non-overlap and diredge constraints,
+ * to resolve:
+ *
+ * for each constraint c:v->cvs[i]:
+ * if exists diredge constraint u->v where u in c:
+ * remap v->cl to cr->v (gap = height(v)/2)
+ *
+ * in = getInConstraints(cvs[i])
+ * for(c : in) {
+ * assert(c.right==cvs[i]);
+ * vin = getOutConstraints(v=c.left)
+ * for(d : vin) {
+ * if(d.left.cluster==i):
+ * tmp = d.left
+ * d.left = d.right
+ * d.right = tmp
+ * d.gap = height(d.right)/2
+ * }
+ * }
+ *
+ */
+ deleteVariable(cvs[i]);
+ }
+ mol += cm[opt->clusters->nclusters];
+ }
+ csolptr = csol = newConstraints(mol);
+ for (i = 0; i < opt->clusters->nclusters + 1; i++) {
+ /* copy constraints into csol */
+ for (j = 0; j < cm[i]; j++) {
+ *csolptr++ = cscl[i][j];
+ }
+ deleteConstraints(0, cscl[i]);
+ }
} else {
- if(k==0) {
- mol = genXConstraints(n,bb,e->vs,&csol,transitiveClosure);
- } else {
- mol = genYConstraints(n,bb,e->vs,&csol);
- }
+ if (k == 0) {
+ mol = genXConstraints(n, bb, e->vs, &csol, transitiveClosure);
+ } else {
+ mol = genYConstraints(n, bb, e->vs, &csol);
+ }
}
- // remove constraints from previous iteration
- if(e->m>0) {
- // can't reuse instance of VPSC when constraints change!
- deleteVPSC(e->vpsc);
- for(i=e->gm==0?0:e->gm;i<e->m;i++) {
- // delete previous overlap constraints
- deleteConstraint(e->cs[i]);
- }
- // just delete the array, not the elements
- if(e->cs!=e->gcs) deleteConstraints(0,e->cs);
+ /* remove constraints from previous iteration */
+ if (e->m > 0) {
+ /* can't reuse instance of VPSC when constraints change! */
+ deleteVPSC(e->vpsc);
+ for (i = e->gm == 0 ? 0 : e->gm; i < e->m; i++) {
+ /* delete previous overlap constraints */
+ deleteConstraint(e->cs[i]);
+ }
+ /* just delete the array, not the elements */
+ if (e->cs != e->gcs)
+ deleteConstraints(0, e->cs);
}
- // if we have no global constraints then the overlap constraints
- // are all we have to worry about.
- // Otherwise, we have to copy the global and overlap constraints
- // into the one array
- if(e->gm==0) {
- e->m = mol;
- e->cs = csol;
+ /* if we have no global constraints then the overlap constraints
+ * are all we have to worry about.
+ * Otherwise, we have to copy the global and overlap constraints
+ * into the one array
+ */
+ if (e->gm == 0) {
+ e->m = mol;
+ e->cs = csol;
} else {
- e->m = mol + e->gm;
- e->cs = newConstraints(e->m);
- for(i=0;i<e->m;i++) {
- if(i<e->gm) {
- e->cs[i]=e->gcs[i];
- } else {
- e->cs[i]=csol[i-e->gm];
- }
- }
- // just delete the array, not the elements
- deleteConstraints(0,csol);
+ e->m = mol + e->gm;
+ e->cs = newConstraints(e->m);
+ for (i = 0; i < e->m; i++) {
+ if (i < e->gm) {
+ e->cs[i] = e->gcs[i];
+ } else {
+ e->cs[i] = csol[i - e->gm];
+ }
+ }
+ /* just delete the array, not the elements */
+ deleteConstraints(0, csol);
}
- if(Verbose)
- fprintf(stderr," generated %d constraints\n",e->m);
- e->vpsc = newIncVPSC(e->nv+e->nldv+e->ndv,e->vs,e->m,e->cs);
+ if (Verbose)
+ fprintf(stderr, " generated %d constraints\n", e->m);
+ e->vpsc = newIncVPSC(e->nv + e->nldv + e->ndv, e->vs, e->m, e->cs);
#ifdef MOSEK
- if(opt->mosek) {
- if(e->mosekEnv!=NULL) {
- mosek_delete(e->mosekEnv);
- }
- e->mosekEnv = mosek_init_sep(e->packedMat,e->nv+e->nldv,e->ndv,e->cs,e->m);
+ if (opt->mosek) {
+ if (e->mosekEnv != NULL) {
+ mosek_delete(e->mosekEnv);
+ }
+ e->mosekEnv =
+ mosek_init_sep(e->packedMat, e->nv + e->nldv, e->ndv, e->cs,
+ e->m);
}
#endif
}
* Statically remove overlaps, that is remove all overlaps by moving each node as
* little as possible.
*/
-void removeoverlaps(int n,float** coords, ipsep_options* opt) {
+void removeoverlaps(int n, float **coords, ipsep_options * opt)
+{
int i;
- CMajEnvVPSC *e = initCMajVPSC(n,NULL,NULL,opt,0);
- generateNonoverlapConstraints(e,1.0,coords,0,true,opt);
+ CMajEnvVPSC *e = initCMajVPSC(n, NULL, NULL, opt, 0);
+ generateNonoverlapConstraints(e, 1.0, coords, 0, true, opt);
solveVPSC(e->vpsc);
- for(i=0;i<n;i++) {
- coords[0][i]=getVariablePos(e->vs[i]);
+ for (i = 0; i < n; i++) {
+ coords[0][i] = getVariablePos(e->vs[i]);
}
- generateNonoverlapConstraints(e,1.0,coords,1,false,opt);
+ generateNonoverlapConstraints(e, 1.0, coords, 1, false, opt);
solveVPSC(e->vpsc);
- for(i=0;i<n;i++) {
- coords[1][i]=getVariablePos(e->vs[i]);
+ for (i = 0; i < n; i++) {
+ coords[1][i] = getVariablePos(e->vs[i]);
}
deleteCMajEnvVPSC(e);
}
/*
unpack the "ordering" array into an array of DigColaLevel
*/
-DigColaLevel* assign_digcola_levels(int *ordering, int n, int *level_inds, int num_divisions) {
- int i,j;
- DigColaLevel *l=N_GNEW(num_divisions+1,DigColaLevel);
- // first level
- l[0].num_nodes=level_inds[0];
- l[0].nodes=N_GNEW(l[0].num_nodes,int);
- for(i=0;i<l[0].num_nodes;i++) {
- l[0].nodes[i]=ordering[i];
- }
- // second through second last level
- for(i=1;i<num_divisions;i++) {
- l[i].num_nodes=level_inds[i]-level_inds[i-1];
- l[i].nodes=N_GNEW(l[i].num_nodes,int);
- for(j=0;j<l[i].num_nodes;j++) {
- l[i].nodes[j]=ordering[level_inds[i-1]+j];
- }
+DigColaLevel *assign_digcola_levels(int *ordering, int n, int *level_inds,
+ int num_divisions)
+{
+ int i, j;
+ DigColaLevel *l = N_GNEW(num_divisions + 1, DigColaLevel);
+ /* first level */
+ l[0].num_nodes = level_inds[0];
+ l[0].nodes = N_GNEW(l[0].num_nodes, int);
+ for (i = 0; i < l[0].num_nodes; i++) {
+ l[0].nodes[i] = ordering[i];
+ }
+ /* second through second last level */
+ for (i = 1; i < num_divisions; i++) {
+ l[i].num_nodes = level_inds[i] - level_inds[i - 1];
+ l[i].nodes = N_GNEW(l[i].num_nodes, int);
+ for (j = 0; j < l[i].num_nodes; j++) {
+ l[i].nodes[j] = ordering[level_inds[i - 1] + j];
}
- // last level
- if (num_divisions>0) {
- l[num_divisions].num_nodes=n-level_inds[num_divisions-1];
- l[num_divisions].nodes=N_GNEW(l[num_divisions].num_nodes,int);
- for(i=0;i<l[num_divisions].num_nodes;i++) {
- l[num_divisions].nodes[i]=ordering[level_inds[num_divisions-1]+i];
- }
+ }
+ /* last level */
+ if (num_divisions > 0) {
+ l[num_divisions].num_nodes = n - level_inds[num_divisions - 1];
+ l[num_divisions].nodes = N_GNEW(l[num_divisions].num_nodes, int);
+ for (i = 0; i < l[num_divisions].num_nodes; i++) {
+ l[num_divisions].nodes[i] =
+ ordering[level_inds[num_divisions - 1] + i];
}
- return l;
+ }
+ return l;
}
-void delete_digcola_levels(DigColaLevel *l, int num_levels) {
- int i;
- for(i=0;i<num_levels;i++) {
- free(l[i].nodes);
- }
- free(l);
+void delete_digcola_levels(DigColaLevel * l, int num_levels)
+{
+ int i;
+ for (i = 0; i < num_levels; i++) {
+ free(l[i].nodes);
+ }
+ free(l);
}
-void print_digcola_levels(FILE *logfile, DigColaLevel *levels, int num_levels) {
- fprintf(logfile,"levels:\n");
- int i,j;
- for(i=0;i<num_levels;i++) {
- fprintf(logfile," l[%d]:",i);
- for(j=0;j<levels[i].num_nodes;j++) {
- fprintf(logfile,"%d ",levels[i].nodes[j]);
- }
- fprintf(logfile,"\n");
+void print_digcola_levels(FILE * logfile, DigColaLevel * levels,
+ int num_levels)
+{
+ int i, j;
+ fprintf(logfile, "levels:\n");
+ for (i = 0; i < num_levels; i++) {
+ fprintf(logfile, " l[%d]:", i);
+ for (j = 0; j < levels[i].num_nodes; j++) {
+ fprintf(logfile, "%d ", levels[i].nodes[j]);
}
+ fprintf(logfile, "\n");
+ }
}
+
/*********************
get number of separation constraints based on the number of nodes in each level
ie, num_sep_constraints = sum_i^{num_levels-1} (|L[i]|+|L[i+1]|)
**********************/
-int get_num_digcola_constraints(DigColaLevel *levels, int num_levels) {
- int i,nc=0;
- for(i=1;i<num_levels;i++) {
- nc+=levels[i].num_nodes+levels[i-1].num_nodes;
- }
- nc+=levels[0].num_nodes+levels[num_levels-1].num_nodes;
- return nc;
+int get_num_digcola_constraints(DigColaLevel * levels, int num_levels)
+{
+ int i, nc = 0;
+ for (i = 1; i < num_levels; i++) {
+ nc += levels[i].num_nodes + levels[i - 1].num_nodes;
+ }
+ nc += levels[0].num_nodes + levels[num_levels - 1].num_nodes;
+ return nc;
}
-#endif /* DIGCOLA */
-
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
+#endif /* IPSEPCOLA */
+/* vim:set shiftwidth=4 ts=8: */
+
+/**
+ * Authors:
+ * Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * This version is released under the CPL (Common Public License) with
+ * the Graphviz distribution.
+ * A version is also available under the LGPL as part of the Adaptagrams
+ * project: http://sourceforge.net/projects/adaptagrams.
+ * If you make improvements or bug fixes to this code it would be much
+ * appreciated if you could also contribute those changes back to the
+ * Adaptagrams repository.
+ */
+
/**********************************************************
* Written by Tim Dwyer for the graphviz package *
* http://www.graphviz.org/ *
#include "digcola.h"
#ifdef MOSEK
#include "mosek_quad_solve.h"
-#endif //MOSEK
+#endif /* MOSEK */
typedef struct CMajEnvVPSC {
float **A;
float *packedMat;
- int nv; // number of actual vars
- int nldv; // number of dummy nodes included in lap matrix
- int ndv; // number of dummy nodes not included in lap matrix
+ int nv; /* number of actual vars */
+ int nldv; /* number of dummy nodes included in lap matrix */
+ int ndv; /* number of dummy nodes not included in lap matrix */
Variable **vs;
int m; /* total number of constraints for next iteration */
int gm; /* number of global constraints */
float *fArray3;
#ifdef MOSEK
MosekEnv *mosekEnv;
-#endif // MOSEK
+#endif /* MOSEK */
} CMajEnvVPSC;
extern CMajEnvVPSC* initCMajVPSC(int n, float *packedMat, vtx_data* graph, ipsep_options *opt, int diredges);
} DigColaLevel;
/*
- unpack the "ordering" array into an array of DigColaLevel (as defined above)
-*/
+ * unpack the "ordering" array into an array of DigColaLevel (as defined above)
+ */
extern DigColaLevel* assign_digcola_levels(int *ordering, int n, int *level_inds, int num_divisions);
extern void delete_digcola_levels(DigColaLevel *l, int num_levels);
extern void print_digcola_levels(FILE* logfile, DigColaLevel *levels, int num_levels);
void IMDS_given_dim(vtx_data* graph, int n, double* given_coords,
double* new_coords, double conj_tol)
{
- int iterations1,iterations2, mat_mult_count1;
+ int iterations2;
int i,j;
DistType** Dij;
float* f_storage = NULL;
double b;
bool converged;
+#if 0
iterations1=mat_mult_count1=0; /* We don't compute the x-axis at all. */
+#endif
Dij = compute_apsp(graph, n);
projects / graphviz / commitdiff