-# $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/cdt \
+ -I$(top_srcdir)/lib/libvpsc
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
+ matrix_ops.h pca.h stress.h quad_prog_solver.h digcola.h \
+ quad_prog_vpsc.h mosek_quad_solve.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
+ compute_hierarchy.c \
+ quad_prog_vpsc.c constrained_majorization_ipsep.c \
+ mosek_quad_solve.c
EXTRA_DIST = Makefile.old
--- /dev/null
+/**********************************************************
+ * Based on constrained_majorization.c
+ *
+ * Perform stress majorization subject
+ * to separation constraints, for background see the paper:
+ * "IPSep-CoLa: An Incremental Procedure for Separation Constraint Layout of Graphs"
+ * by Tim Dwyer, Yehuda Koren and Kim Marriott
+ *
+ * Available separation constraints so far are:
+ * o Directed edge constraints
+ * o Node non-overlap constraints
+ * o Cluster containment constraints
+ * o Cluster/node non-overlap constraints
+ *
+ * Tim Dwyer, 2006
+ **********************************************************/
+
+#include "digcola.h"
+#ifdef DIGCOLA
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stdio.h>
+#include <float.h>
+#include "stress.h"
+#include "dijkstra.h"
+#include "bfs.h"
+#include "matrix_ops.h"
+#include "kkutils.h"
+#include "conjgrad.h"
+#include <csolve_VPSC.h>
+#include "quad_prog_vpsc.h"
+#include "quad_prog_solver.h"
+#include "matrix_ops.h"
+
+#define localConstrMajorIterations 1000
+
+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 dim, /* Dimemsionality of layout */
+ int model, /* difference model */
+ int maxi, /* max iterations */
+ ipsep_options* opt
+)
+{
+ int iterations = 0; /* Output: number of iteration of the process */
+
+ /*************************************************
+ ** Computation of full, dense, unrestricted k-D **
+ ** stress minimization by majorization **
+ ** 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;
+
+ initLayout(graph, n, dim, d_coords);
+ 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);
+ }
+ }
+
+ /****************************************************
+ ** Compute the all-pairs-shortest-distances matrix **
+ ****************************************************/
+
+ 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);
+ } 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");
+ }
+ }
+ if (!Dij) {
+ 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];
+ }
+ }
+
+ /* 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 (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 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;
+ 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;
+ }
+ /* compute diagonal entries */
+ 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 (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;
+ }
+ }
+
+ /* compute constant term in stress sum
+ * which is \sum_{i<j} w_{ij}d_{ij}^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];
+ }
+
+ /* 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");
+ }
+ 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) {
+#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 (Verbose) {
+ fprintf(stderr, "\nfinal e = %f %d iterations %.2f sec\n",
+ new_stress, iterations, elapsed_sec());
+ }
+ deleteCMajEnvVPSC(cMajEnvHor);
+ deleteCMajEnvVPSC(cMajEnvVrt);
+
+ 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);
+ }
+ 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 :
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) {
#endif
} vtx_data;
+ typedef 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 */
+ } cluster_data;
+
+
typedef int DistType; /* must be signed!! */
#ifdef UNUSED
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);
+typedef struct ipsep_options {
+ int diredges; /* 1=generate directed edge constraints */
+ /* 2=generate directed hierarchy level constraints (DiG-CoLa) */
+ double edge_gap; /* amount to force vertical separation of */
+ /* start/end nodes */
+ int noverlap; /* 1=generate non-overlap constraints */
+ /* 2=remove overlaps after layout */
+ pointf gap; /* hor and vert gap to enforce when removing overlap*/
+ pointf* nsize; /* node widths and heights */
+ cluster_data* clusters;
+ /* list of node indices for each cluster */
+#ifdef MOSEK
+ int mosek; /* use Mosek as constraint optimization engine */
+#endif //MOSEK
+} ipsep_options;
+
+// stress majorization, for Constraint Layout
+extern int stress_majorization_cola(vtx_data*, int, int, double**, int, int, int, ipsep_options*);
#endif
#endif
--- /dev/null
+/*
+ * Interface to Mosek (www.mosek.com) quadratic programming solver for solving
+ * instances of the "Variable Placement with Separation Constraints" problem,
+ * and also DiG-CoLa style level constraints.
+ *
+ * Tim Dwyer, 2006
+ */
+#ifdef MOSEK
+#include <stdio.h>
+#include <assert.h>
+#include "defs.h"
+#include "mosek_quad_solve.h"
+#include "quad_prog_vpsc.h"
+
+//#define DUMP_CONSTRAINTS
+//#define EQUAL_WIDTH_LEVELS
+
+static FILE *logfile;
+static void MSKAPI printstr(void *handle, char str[])
+{
+ fprintf(logfile,"%s",str);
+} /* printstr */
+
+/**********************
+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
+n: number of nodes (length of coords array)
+ordering: array containing sequences of nodes for each level,
+ ie, ordering[levels[i]] is first node of (i+1)th level
+level_indexes: array of starting node for each level in ordering
+ ie, levels[i] is index to first node of (i+1)th level
+ also, levels[0] is number of nodes in first level
+ and, levels[i]-levels[i-1] is number of nodes in ith level
+ and, n - levels[num_divisions-1] is number of nodes in last level
+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)
+{
+ 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);
+#ifdef DUMP_CONSTRAINTS
+ 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.
+#ifdef EQUAL_WIDTH_LEVELS
+ 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++;
+ }
+ }
+#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,"]\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);
+ }
+
+ /* 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++;
+ }
+#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++;
+ }
+#endif
+ 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);
+ }
+ }
+ }
+ delete_digcola_levels(levels,num_levels);
+ return mskEnv;
+}
+/*
+b: coefficients of linear part of optimisation function
+n: number of nodes
+coords: optimal y* vector, coord[i] is coordinate of node[i]
+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]);
+ }
+#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");
+#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);
+
+#ifdef DUMP_CONSTRAINTS
+ fprintf(logfile,"Primal solution\n");
+#endif
+ coords[0]=0;
+ for(j=0; j<mskEnv->num_variables; ++j) {
+#ifdef DUMP_CONSTRAINTS
+ 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];
+ }
+ }
+ }
+ 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)
+{
+ 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++;
+ }
+ }
+#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,"]\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);
+ }
+
+ /* 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 );
+ }
+ }
+ }
+ 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);
+
+#ifdef DUMP_CONSTRAINTS
+ fprintf(logfile,"Primal solution\n");
+#endif
+ coords[0]=0;
+ for(j=1; j<=n; j++) {
+#ifdef DUMP_CONSTRAINTS
+ fprintf(logfile,"x[%d]: %.2f\n",j,mskEnv->xx[j-1]);
+#endif
+ coords[j]=-mskEnv->xx[j-1];
+ }
+ }
+ 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);
+
+ if(logfile) {
+ fclose(logfile);
+ logfile=NULL;
+ }
+ free(mskEnv->qval);
+ free(mskEnv->qsubi);
+ free(mskEnv->qsubj);
+ free(mskEnv->xx);
+ free(mskEnv);
+}
+#endif /* MOSEK */
--- /dev/null
+#ifdef MOSEK
+#ifndef _QSOLVE_H_
+#define _QSOLVE_H_
+
+#include <mosek.h> /* Include the MOSEK definition file. */
+#include "types.h"
+#include <csolve_VPSC.h>
+
+typedef struct {
+ int r;
+ MSKenv_t env;
+ MSKtask_t task;
+ double *qval;
+ int *qsubi,*qsubj;
+ double *xx;
+ int num_variables;
+} MosekEnv;
+
+MosekEnv* mosek_init_hier(float* lap, int n,int *ordering,int *level_indexes,int num_divisions, float separation);
+void mosek_quad_solve_hier(MosekEnv*,float *b,int n,float* coords, float *hierarchy_boundaries);
+MosekEnv* mosek_init_sep(float* lap, int nv, int ndv, Constraint** cs, int m);
+void mosek_quad_solve_sep(MosekEnv*,int n,float *b,float* coords);
+void mosek_delete(MosekEnv*);
+
+
+#endif // _QSOLVE_H_
+#endif // MOSEK
#define MODE_KK 0
#define MODE_MAJOR 1
#define MODE_HIER 2
+#define MODE_IPSEP 3
#define INIT_SELF 0
#define INIT_REGULAR 1
}
}
+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;
+ edge_t *me;
+ // 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;
+ }
+ cluster_data *cdata = GNEW(cluster_data);
+ cdata->ntoplevel=agnnodes(g);
+ mm = mastergraph->meta_node;
+ mg = mm->graph;
+ for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
+ mn = me->head;
+ subg = agusergraph(mn);
+ if (!strncmp(subg->name, "cluster", 7)) {
+ nclusters++;
+ }
+ }
+ 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);
+ 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));
+ for (n = agfstnode(subg); n; n = agnxtnode(subg, n)) {
+ node_t *gn;
+ 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);
+ *c++=ind;
+ assigned[ind]=true;
+ cdata->ntoplevel--;
+ }
+ }
+ }
+ cdata->bb=N_GNEW(cdata->nclusters,boxf);
+ cdata->toplevel=N_GNEW(cdata->ntoplevel,int);
+ for(i=j=0;i<agnnodes(g);i++) {
+ if(!assigned[i]) {
+ cdata->toplevel[j++]=i;
+ }
+ }
+ assert(cdata->ntoplevel==agnnodes(g)-cdata->nvars);
+ return cdata;
+}
+static void freeClusterData(cluster_data *c) {
+ if(c->nclusters>0) {
+ free(c->clusters[0]);
+ free(c->clusters);
+ free(c->clustersizes);
+ free(c->toplevel);
+ free(c->bb);
+ }
+ free(c);
+}
/* 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;
+ else if (streq(str, "ipsep"))
+ mode = MODE_IPSEP;
#endif
else
agerr(AGWARN,
* dfs for breaking cycles in vtxdata
*/
static void
-dfsCycle (vtx_data* graph, int i)
+dfsCycle (vtx_data* graph, int i,int mode)
{
node_t *np, *hp;
int j, e, f;
j = graph[i].edges[e];
hp = graph[j].np;
if (ND_onstack(hp)) { /* back edge: reverse it */
- graph[i].edists[e] = 1.0;
- for (f = 1; (f < graph[j].nedges) &&(graph[j].edges[f] != i); f++) ;
- assert (f < graph[j].nedges);
- graph[j].edists[f] = -1.0;
+ // 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;
+ for (f = 1; (f < graph[j].nedges) &&(graph[j].edges[f] != i); f++) ;
+ assert (f < graph[j].nedges);
+ graph[j].edists[f] = -1.0;
}
- else if (ND_mark(hp) == FALSE) dfsCycle(graph, j);
+ else if (ND_mark(hp) == FALSE) dfsCycle(graph, j, mode);
}
ND_onstack(np) = FALSE;
* Do a dfs of the vtx_data, looking for cycles, reversing edges.
*/
static void
-acyclic (vtx_data* graph, int nv)
+acyclic (vtx_data* graph, int nv, int mode)
{
int i;
node_t* np;
}
for (i = 0; i < nv; i++) {
if (ND_mark(graph[i].np)) continue;
- dfsCycle (graph, i);
+ dfsCycle (graph, i, mode);
}
}
haveLen = (agindex(g->root->proto->e, "len") >= 0);
haveWt = (E_weight != 0);
}
- if (mode == MODE_HIER)
+ if (mode == MODE_HIER || mode == MODE_IPSEP)
haveDir = TRUE;
else
haveDir = FALSE;
else if (haveDir)
*ewgts++ = 1.0;
#ifdef DIGCOLA
- if (haveDir) {
- *edists++ = (np == ep->head ? 1.0 : -1.0);
- }
+ if (haveDir) {
+ char *s=agget(ep,"dir");
+ if(s&&!strncmp(s,"none",4)) {
+ *edists++ = 0;
+ } else {
+ *edists++ = (np == ep->head ? 1.0 : -1.0);
+ }
+ }
#endif
i_nedges++;
}
#ifdef DIGCOLA
if (haveDir) {
/* Make graph acyclic */
- acyclic (graph, nv);
+ acyclic (graph, nv, mode);
}
#endif
* weight?
*/
static void
-majorization(graph_t * g, int nv, int mode, int model, int dim, int steps)
+majorization(graph_t *mg, graph_t * g, int nv, int mode, int model, int dim, int steps)
{
double **coords;
int ne;
int i;
node_t *v;
vtx_data *gp;
+ cluster_data *cs;
+
int init;
init = checkStart(g, nv, (mode == MODE_HIER ? INIT_SELF : INIT_RANDOM));
start_timer();
}
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) {
- double lgap = late_double(g, agfindattr(g, "levelsgap"), 0.0, -MAXDOUBLE);
- stress_majorization_with_hierarchy(gp, nv, ne, coords, Ndim,
- (init == INIT_SELF), model, MaxIter, lgap);
+ if (mode == MODE_HIER || mode == MODE_IPSEP) {
+ 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 {
+ char* str;
+ ipsep_options opt;
+ pointf nsize[nv];
+ 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;
+ str = agget(g, "diredgeconstraints");
+ if(str && !strncmp(str,"true",4)) {
+ opt.diredges = 1;
+ if(Verbose)
+ fprintf(stderr,"Generating Edge Constraints...\n");
+ } else if(str && !strncmp(str,"hier",4)) {
+ opt.diredges = 2;
+ if(Verbose)
+ fprintf(stderr,"Generating DiG-CoLa Edge Constraints...\n");
+ }
+ str = agget(g, "overlapconstraints");
+ if(str && !strncmp(str,"true",4)) {
+ opt.noverlap = 1;
+ if(Verbose)
+ fprintf(stderr,"Generating Non-overlap Constraints...\n");
+ } else if(str && !strncmp(str,"post",4)) {
+ opt.noverlap = 2;
+ if(Verbose)
+ fprintf(stderr,"Removing overlaps as postprocess...\n");
+ }
+#ifdef MOSEK
+ str = agget(g, "mosek");
+ if(str && !strncmp(str,"true",4)) {
+ opt.mosek = 1;
+ 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);
+ if(Verbose)
+ fprintf(stderr,"gap=%f,%f\n",opt.gap.x,opt.gap.y);
+ } else {
+ opt.gap.x =opt.gap.y = atof(str);
+ }
+ }
+ for (i=0, v = agfstnode(g); v; v = agnxtnode(g, v),i++) {
+ nsize[i].x=ND_width(v);
+ nsize[i].y=ND_height(v);
+ }
+
+ stress_majorization_cola(gp, nv, ne, coords, Ndim, model, MaxIter, &opt);
+ }
}
else
#endif
}
}
freeGraphData(gp);
+ freeClusterData(cs);
free(coords[0]);
free(coords);
}
/* neatoLayout:
* Use stress optimization to layout a single component
*/
-void neatoLayout(Agraph_t * g, int layoutMode, int layoutModel)
+void neatoLayout(Agraph_t * mg, Agraph_t * g, int layoutMode, int layoutModel)
{
int nG;
char *str;
if (!nG)
return;
if (layoutMode)
- majorization(g, nG, layoutMode, layoutModel, Ndim, MaxIter);
+ majorization(mg, g, nG, layoutMode, layoutModel, Ndim, MaxIter);
else
kkNeato(g, nG, layoutModel);
}
for (i = 0; i < n_cc; i++) {
gc = cc[i];
nodeInduce(gc);
- neatoLayout(gc, layoutMode, model);
+ neatoLayout(g, gc, layoutMode, model);
adjustNodes(gc);
}
if (n_cc > 1) {
free_scan_graph(gc);
agdelete(g, gc);
}
+ {
+ graph_t *mg, *subg;
+ node_t *mm, *mn;
+ edge_t *me;
+ mm = g->meta_node;
+ mg = mm->graph;
+ for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
+ mn = me->head;
+ subg = agusergraph(mn);
+ if (!strncmp(subg->name, "cluster", 7)) {
+ add_cluster(g,subg);
+ compute_bb(subg);
+ }
+ }
+ }
} else {
- neatoLayout(g, layoutMode, model);
+ neatoLayout(g, g, layoutMode, model);
adjustNodes(g);
addZ (g);
spline_edges(g);
extern void neato_init_edge(Agedge_t *);
extern void neato_init_node(Agnode_t *);
extern void neato_layout(Agraph_t * g);
- extern void neatoLayout(Agraph_t * g, int layoutMode, int layoutModel);
+ extern void neatoLayout(Agraph_t *mg, Agraph_t * g, int layoutMode, int layoutModel);
extern void neato_init_graph(graph_t * g);
extern void neato_init_graphn(graph_t * g, int);
extern int Plegal_arrangement(Ppoly_t ** polys, int n_polys);
return counter;
}
+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;
+}
+/*
+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;
+
+#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];
+ }
+ }
+ FILE *logfile = fopen("constrained_majorization_log","a");
+
+ 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]];
+ }
+ }
+ }
+ // 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];
+ }
+ 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);
+ }
+ */
+#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;
+#endif
+ if(test>quad_prog_tol) {
+ converged=false;
+ }
+ }
+#ifdef CONMAJ_LOGGING
+ call_no++;
+ fclose(logfile);
+#endif
+ return counter;
+}
+
int
constrained_majorization_new_with_gaps(CMajEnv *e, float* b, float** coords,
int ndims, int cur_axis, int max_iterations, float* hierarchy_boundaries,
extern int constrained_majorization_new(CMajEnv*, float*, float**,
int, int, int, float*, float);
+extern int constrained_majorization_new_with_gaps(CMajEnv*, float*, float**,
+ int, int, int, float*, float);
+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);
extern void deleteCMajEnv(CMajEnv *e);
--- /dev/null
+/**********************************************************
+ *
+ * Solve a quadratic function f(X) = X' e->A X + b X
+ * subject to a set of separation constraints e->cs
+ *
+ * Tim Dwyer, 2006
+ **********************************************************/
+
+#include "digcola.h"
+#ifdef DIGCOLA
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stdio.h>
+#include <float.h>
+#include <assert.h>
+#include "matrix_ops.h"
+#include "kkutils.h"
+#include <csolve_VPSC.h>
+#include "quad_prog_vpsc.h"
+#include "quad_prog_solver.h"
+
+//#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 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;
+#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");
+ }
+#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];
+ }
+ }
+ FILE *logfile = fopen("constrained_majorization_log","a");
+
+ //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]);
+ }
+#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;
+#endif
+ if(test>quad_prog_tol) {
+ converged=false;
+ }
+ }
+#ifdef CONMAJ_LOGGING
+ call_no++;
+ fclose(logfile);
+#endif
+ return counter;
+}
+
+/*
+ * Set up environment and global constraints (dir-edge constraints, containment constraints
+ * etc).
+ *
+ * diredges: 0=no dir edge constraints
+ * 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)
+{
+ 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;
+ /* if we have clusters then we'll need two constraints for each var in
+ * a cluster */
+ 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->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);
+ }
+ }
+ //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;
+ }
+ 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);
+ }
+#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;
+}
+void
+deleteCMajEnvVPSC(CMajEnvVPSC *e)
+{
+ int i;
+ 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);
+ }
+ free (e->fArray1);
+ free (e->fArray2);
+ free (e->fArray3);
+#ifdef MOSEK
+ if(e->mosekEnv) {
+ mosek_delete(e->mosekEnv);
+ }
+#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
+) {
+ Constraint **csol, **csolptr;
+ 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;
+ }
+ 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;
+ }
+ 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]);
+ }
+ } else {
+ 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);
+ }
+ // 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);
+ }
+ 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);
+ }
+#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) {
+ int i;
+ 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]);
+ }
+ generateNonoverlapConstraints(e,1.0,coords,1,false,opt);
+ solveVPSC(e->vpsc);
+ 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];
+ }
+ }
+ // 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;
+}
+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");
+ }
+}
+/*********************
+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;
+}
+
+#endif /* DIGCOLA */
+
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
--- /dev/null
+/**********************************************************
+* Written by Tim Dwyer for the graphviz package *
+* http://www.graphviz.org/ *
+* *
+**********************************************************/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#ifndef _QUAD_PROG_VPSC_H_
+#define _QUAD_PROG_VPSC_H_
+
+#ifdef DIGCOLA
+
+#include "defs.h"
+#include "digcola.h"
+#ifdef MOSEK
+#include "mosek_quad_solve.h"
+#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
+ Variable **vs;
+ int m; /* total number of constraints for next iteration */
+ int gm; /* number of global constraints */
+ Constraint **cs;
+ /* global constraints are persistant throughout optimisation process */
+ Constraint **gcs;
+ VPSC *vpsc;
+ float *fArray1; /* utility arrays - reusable memory */
+ float *fArray2;
+ float *fArray3;
+#ifdef MOSEK
+ MosekEnv *mosekEnv;
+#endif // MOSEK
+} CMajEnvVPSC;
+
+extern CMajEnvVPSC* initCMajVPSC(int n, float *packedMat, vtx_data* graph, ipsep_options *opt, int diredges);
+
+extern int constrained_majorization_vpsc(CMajEnvVPSC*, float*, float*, int);
+
+extern void deleteCMajEnvVPSC(CMajEnvVPSC *e);
+extern void generateNonoverlapConstraints(
+ CMajEnvVPSC* e,
+ float nsizeScale,
+ float** coords,
+ int k,
+ bool transitiveClosure,
+ ipsep_options* opt
+);
+
+extern void removeoverlaps(int,float**,ipsep_options*);
+
+typedef struct {
+ int *nodes;
+ int num_nodes;
+} DigColaLevel;
+
+/*
+ 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);
+int get_num_digcola_constraints(DigColaLevel *levels, int num_levels);
+#endif
+
+#endif /* _QUAD_PROG_VPSC_H_ */
+
+#ifdef __cplusplus
+}
+#endif
projects / graphviz / commitdiff