Graphviz  2.41.20171026.1811
mosek_quad_solve.c
Go to the documentation of this file.
1 /* $Id$ $Revision$ */
2 /* vim:set shiftwidth=4 ts=8: */
3 
19 /*
20  * Interface to Mosek (www.mosek.com) quadratic programming solver for solving
21  * instances of the "Variable Placement with Separation Constraints" problem,
22  * and also DiG-CoLa style level constraints.
23  *
24  * Tim Dwyer, 2006
25  */
26 #ifdef MOSEK
27 #include <stdio.h>
28 #include <assert.h>
29 #include "defs.h"
30 #include "mosek_quad_solve.h"
31 #include "quad_prog_vpsc.h"
32 
33 /* #define DUMP_CONSTRAINTS */
34 /* #define EQUAL_WIDTH_LEVELS */
35 
36 static FILE *logfile;
37 static void MSKAPI printstr(void *handle, char str[])
38 {
39  fprintf(logfile, "%s", str);
40 }
41 
42 #define INIT_sub_val(a,b) \
43  MSKidxt subi[2]; \
44  double vali[2]; \
45  subi[0] = a; \
46  subi[1] = b; \
47  vali[0] = 1.0; \
48  vali[1] = -1.0;
49 
50 #define INIT_sub_val3(a,b,c) \
51  MSKidxt subi[3]; \
52  double vali[3]; \
53  subi[0] = a; \
54  subi[1] = b; \
55  subi[2] = c; \
56  vali[0] = 1.0; \
57  vali[1] = -2.0; \
58  vali[2] = 1.0;
59 
60 /**********************
61 lap: the upper RHS of the symmetric graph laplacian matrix which will be transformed
62  to the hessian of the non-linear part of the optimisation function
63 n: number of nodes (length of coords array)
64 ordering: array containing sequences of nodes for each level,
65  ie, ordering[levels[i]] is first node of (i+1)th level
66 level_indexes: array of starting node for each level in ordering
67  ie, levels[i] is index to first node of (i+1)th level
68  also, levels[0] is number of nodes in first level
69  and, levels[i]-levels[i-1] is number of nodes in ith level
70  and, n - levels[num_divisions-1] is number of nodes in last level
71 num_divisions: number of divisions between levels, ie number of levels - 1
72 separation: the minimum separation between nodes on different levels
73 ***********************/
74 MosekEnv *mosek_init_hier(float *lap, int n, int *ordering,
75  int *level_indexes, int num_divisions,
76  float separation)
77 {
78  int count = 0;
79  int i, j, num_levels = num_divisions + 1;
80  int num_constraints;
81  MosekEnv *mskEnv = GNEW(MosekEnv);
82  DigColaLevel *levels;
83  int nonzero_lapsize = (n * (n - 1)) / 2;
84  /* vars for nodes (except x0) + dummy nodes between levels
85  * x0 is fixed at 0, and therefore is not included in the opt problem
86  * add 2 more vars for top and bottom constraints
87  */
88  mskEnv->num_variables = n + num_divisions + 1;
89 
90  logfile = fopen("quad_solve_log", "w");
91  levels = assign_digcola_levels(ordering, n, level_indexes, num_divisions);
92 #ifdef DUMP_CONSTRAINTS
93  print_digcola_levels(logfile, levels, num_levels);
94 #endif
95 
96  /* nonlinear coefficients matrix of objective function */
97  /* int lapsize=mskEnv->num_variables+(mskEnv->num_variables*(mskEnv->num_variables-1))/2; */
98  mskEnv->qval = N_GNEW(nonzero_lapsize, double);
99  mskEnv->qsubi = N_GNEW(nonzero_lapsize, int);
100  mskEnv->qsubj = N_GNEW(nonzero_lapsize, int);
101 
102  /* solution vector */
103  mskEnv->xx = N_GNEW(mskEnv->num_variables, double);
104 
105  /* constraint matrix */
106  separation /= 2.0; /* separation between each node and it's adjacent constraint */
107  num_constraints = get_num_digcola_constraints(levels,
108  num_levels) + num_divisions + 1;
109  /* constraints of the form x_i - x_j >= sep so 2 non-zero entries per constraint in LHS matrix
110  * except x_0 (fixed at 0) constraints which have 1 nz val each.
111  */
112 #ifdef EQUAL_WIDTH_LEVELS
113  num_constraints += num_divisions;
114 #endif
115  /* pointer to beginning of nonzero sequence in a column */
116 
117  for (i = 0; i < n - 1; i++) {
118  for (j = i; j < n - 1; j++) {
119  mskEnv->qval[count] = -2 * lap[count + n];
120  assert(mskEnv->qval[count] != 0);
121  mskEnv->qsubi[count] = j;
122  mskEnv->qsubj[count] = i;
123  count++;
124  }
125  }
126 #ifdef DUMP_CONSTRAINTS
127  fprintf(logfile, "Q=[");
128  int lapcntr = n;
129  for (i = 0; i < mskEnv->num_variables; i++) {
130  if (i != 0)
131  fprintf(logfile, ";");
132  for (j = 0; j < mskEnv->num_variables; j++) {
133  if (j < i || i >= n - 1 || j >= n - 1) {
134  fprintf(logfile, "0 ");
135  } else {
136  fprintf(logfile, "%f ", -2 * lap[lapcntr++]);
137  }
138  }
139  }
140  fprintf(logfile, "]\nQ=Q-diag(diag(Q))+Q'\n");
141 #endif
142  fprintf(logfile, "\n");
143  /* Make the mosek environment. */
144  mskEnv->r = MSK_makeenv(&mskEnv->env, NULL, NULL, NULL, NULL);
145 
146  /* Check whether the return code is ok. */
147  if (mskEnv->r == MSK_RES_OK) {
148  /* Directs the log stream to the user
149  * specified procedure 'printstr'.
150  */
151  MSK_linkfunctoenvstream(mskEnv->env, MSK_STREAM_LOG, NULL,
152  printstr);
153  }
154 
155  /* Initialize the environment. */
156  mskEnv->r = MSK_initenv(mskEnv->env);
157  if (mskEnv->r == MSK_RES_OK) {
158  /* Make the optimization task. */
159  mskEnv->r =
160  MSK_maketask(mskEnv->env, num_constraints,
161  mskEnv->num_variables, &mskEnv->task);
162 
163  if (mskEnv->r == MSK_RES_OK) {
164  int c_ind = 0;
165  int c_var = n - 1;
166  mskEnv->r =
167  MSK_linkfunctotaskstream(mskEnv->task, MSK_STREAM_LOG,
168  NULL, printstr);
169  /* Resize the task. */
170  if (mskEnv->r == MSK_RES_OK)
171  mskEnv->r = MSK_resizetask(mskEnv->task, num_constraints, mskEnv->num_variables, 0, /* no cones!! */
172  /* each constraint applies to 2 vars */
173  2 * num_constraints +
174  num_divisions, nonzero_lapsize);
175 
176  /* Append the constraints. */
177  if (mskEnv->r == MSK_RES_OK)
178  mskEnv->r = MSK_append(mskEnv->task, 1, num_constraints);
179 
180  /* Append the variables. */
181  if (mskEnv->r == MSK_RES_OK)
182  mskEnv->r =
183  MSK_append(mskEnv->task, 0, mskEnv->num_variables);
184  /* Put variable bounds. */
185  for (j = 0;
186  j < mskEnv->num_variables && mskEnv->r == MSK_RES_OK; ++j)
187  mskEnv->r =
188  MSK_putbound(mskEnv->task, 0, j, MSK_BK_RA,
189  -MSK_INFINITY, MSK_INFINITY);
190  for (j = 0; j < levels[0].num_nodes && mskEnv->r == MSK_RES_OK;
191  j++) {
192  int node = levels[0].nodes[j] - 1;
193  if (node >= 0) {
194  INIT_sub_val(c_var,node);
195  mskEnv->r =
196  MSK_putavec(mskEnv->task, 1, c_ind, 2, subi, vali);
197  } else {
198  /* constraint for y0 (fixed at 0) */
199  mskEnv->r =
200  MSK_putaij(mskEnv->task, c_ind, c_var, 1.0);
201  }
202  mskEnv->r =
203  MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
204  separation, MSK_INFINITY);
205  c_ind++;
206  }
207  for (i = 0; i < num_divisions && mskEnv->r == MSK_RES_OK; i++) {
208  c_var = n + i;
209  for (j = 0;
210  j < levels[i].num_nodes && mskEnv->r == MSK_RES_OK;
211  j++) {
212  /* create separation constraint a>=b+separation */
213  int node = levels[i].nodes[j] - 1;
214  if (node >= 0) { /* no constraint for fixed node */
215  INIT_sub_val(node,c_var);
216  mskEnv->r =
217  MSK_putavec(mskEnv->task, 1, c_ind, 2, subi,
218  vali);
219  } else {
220  /* constraint for y0 (fixed at 0) */
221  mskEnv->r =
222  MSK_putaij(mskEnv->task, c_ind, c_var, -1.0);
223  }
224  mskEnv->r =
225  MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
226  separation, MSK_INFINITY);
227  c_ind++;
228  }
229  for (j = 0;
230  j < levels[i + 1].num_nodes
231  && mskEnv->r == MSK_RES_OK; j++) {
232  int node = levels[i + 1].nodes[j] - 1;
233  if (node >= 0) {
234  INIT_sub_val(c_var,node);
235  mskEnv->r =
236  MSK_putavec(mskEnv->task, 1, c_ind, 2, subi,
237  vali);
238  } else {
239  /* constraint for y0 (fixed at 0) */
240  mskEnv->r =
241  MSK_putaij(mskEnv->task, c_ind, c_var, 1.0);
242  }
243  mskEnv->r =
244  MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
245  separation, MSK_INFINITY);
246  c_ind++;
247  }
248  }
249  c_var = n + i;
250  for (j = 0; j < levels[i].num_nodes && mskEnv->r == MSK_RES_OK;
251  j++) {
252  /* create separation constraint a>=b+separation */
253  int node = levels[i].nodes[j] - 1;
254  if (node >= 0) { /* no constraint for fixed node */
255  INIT_sub_val(node,c_var);
256  mskEnv->r =
257  MSK_putavec(mskEnv->task, 1, c_ind, 2, subi, vali);
258  } else {
259  /* constraint for y0 (fixed at 0) */
260  mskEnv->r =
261  MSK_putaij(mskEnv->task, c_ind, c_var, -1.0);
262  }
263  mskEnv->r =
264  MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO,
265  separation, MSK_INFINITY);
266  c_ind++;
267  }
268  /* create constraints preserving the order of dummy vars */
269  for (i = 0; i < num_divisions + 1 && mskEnv->r == MSK_RES_OK;
270  i++) {
271  int c_var = n - 1 + i, c_var2 = c_var + 1;
272  INIT_sub_val(c_var,c_var2);
273  mskEnv->r =
274  MSK_putavec(mskEnv->task, 1, c_ind, 2, subi, vali);
275  mskEnv->r =
276  MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_LO, 0,
277  MSK_INFINITY);
278  c_ind++;
279  }
280 #ifdef EQUAL_WIDTH_LEVELS
281  for (i = 1; i < num_divisions + 1 && mskEnv->r == MSK_RES_OK;
282  i++) {
283  int c_var = n - 1 + i, c_var_lo = c_var - 1, c_var_hi =
284  c_var + 1;
285  INIT_sub_val3(c_var_lo, c_var, c_var_h);
286  mskEnv->r =
287  MSK_putavec(mskEnv->task, 1, c_ind, 3, subi, vali);
288  mskEnv->r =
289  MSK_putbound(mskEnv->task, 1, c_ind, MSK_BK_FX, 0, 0);
290  c_ind++;
291  }
292 #endif
293  assert(c_ind == num_constraints);
294 #ifdef DUMP_CONSTRAINTS
295  fprintf(logfile, "A=[");
296  for (i = 0; i < num_constraints; i++) {
297  if (i != 0)
298  fprintf(logfile, ";");
299  for (j = 0; j < mskEnv->num_variables; j++) {
300  double aij;
301  MSK_getaij(mskEnv->task, i, j, &aij);
302  fprintf(logfile, "%f ", aij);
303  }
304  }
305  fprintf(logfile, "]\n");
306  fprintf(logfile, "b=[");
307  for (i = 0; i < num_constraints; i++) {
308  fprintf(logfile, "%f ", separation);
309  }
310  fprintf(logfile, "]\n");
311 #endif
312  if (mskEnv->r == MSK_RES_OK) {
313  /*
314  * The lower triangular part of the Q
315  * matrix in the objective is specified.
316  */
317  mskEnv->r =
318  MSK_putqobj(mskEnv->task, nonzero_lapsize,
319  mskEnv->qsubi, mskEnv->qsubj,
320  mskEnv->qval);
321  }
322  }
323  }
324  delete_digcola_levels(levels, num_levels);
325  return mskEnv;
326 }
327 
328 /*
329 b: coefficients of linear part of optimisation function
330 n: number of nodes
331 coords: optimal y* vector, coord[i] is coordinate of node[i]
332 hierarchy_boundaries: y coord of boundaries between levels
333  (ie, solution values for the dummy variables used in constraints)
334 */
335 void mosek_quad_solve_hier(MosekEnv * mskEnv, float *b, int n,
336  float *coords, float *hierarchy_boundaries)
337 {
338  int i, j;
339  for (i = 1; i < n && mskEnv->r == MSK_RES_OK; i++) {
340  mskEnv->r = MSK_putcj(mskEnv->task, i - 1, -2 * b[i]);
341  }
342 #ifdef DUMP_CONSTRAINTS
343  fprintf(logfile, "x0=[");
344  for (j = 0; j < mskEnv->num_variables; j++) {
345  fprintf(logfile, "%f ", j < n ? b[j] : 0);
346  }
347  fprintf(logfile, "]\n");
348  fprintf(logfile, "f=[");
349  double *c = N_GNEW(mskEnv->num_variables, double);
350  MSK_getc(mskEnv->task, c);
351  for (j = 0; j < mskEnv->num_variables; j++) {
352  fprintf(logfile, "%f ", c[j]);
353  }
354  free(c);
355  fprintf(logfile, "]\n");
356 #endif
357  if (mskEnv->r == MSK_RES_OK)
358  mskEnv->r = MSK_optimize(mskEnv->task);
359 
360  if (mskEnv->r == MSK_RES_OK) {
361  MSK_getsolutionslice(mskEnv->task,
362  MSK_SOL_ITR,
363  MSK_SOL_ITEM_XX,
364  0, mskEnv->num_variables, mskEnv->xx);
365 
366 #ifdef DUMP_CONSTRAINTS
367  fprintf(logfile, "Primal solution\n");
368 #endif
369  coords[0] = 0;
370  for (j = 0; j < mskEnv->num_variables; ++j) {
371 #ifdef DUMP_CONSTRAINTS
372  fprintf(logfile, "x[%d]: %.2f\n", j, mskEnv->xx[j]);
373 #endif
374  if (j < n - 1) {
375  coords[j + 1] = -mskEnv->xx[j];
376  } else if (j >= n && j < mskEnv->num_variables - 1) {
377  hierarchy_boundaries[j - n] = -mskEnv->xx[j];
378  }
379  }
380  }
381  fprintf(logfile, "Return code: %d\n", mskEnv->r);
382 }
383 
384 /**********************
385 lap: the upper RHS of the symmetric graph laplacian matrix which will be transformed
386  to the hessian of the non-linear part of the optimisation function
387  has dimensions num_variables, dummy vars do not have entries in lap
388 cs: array of pointers to separation constraints
389 ***********************/
390 MosekEnv *mosek_init_sep(float *lap, int num_variables, int num_dummy_vars,
391  Constraint ** cs, int num_constraints)
392 {
393  int i, j;
394  MosekEnv *mskEnv = GNEW(MosekEnv);
395  int count = 0;
396  int nonzero_lapsize = num_variables * (num_variables - 1) / 2;
397  /* fix var 0 */
398  mskEnv->num_variables = num_variables + num_dummy_vars - 1;
399 
400  fprintf(stderr, "MOSEK!\n");
401  logfile = fopen("quad_solve_log", "w");
402 
403  /* nonlinear coefficients matrix of objective function */
404  mskEnv->qval = N_GNEW(nonzero_lapsize, double);
405  mskEnv->qsubi = N_GNEW(nonzero_lapsize, int);
406  mskEnv->qsubj = N_GNEW(nonzero_lapsize, int);
407 
408  /* solution vector */
409  mskEnv->xx = N_GNEW(mskEnv->num_variables, double);
410 
411  /* pointer to beginning of nonzero sequence in a column */
412 
413  for (i = 0; i < num_variables - 1; i++) {
414  for (j = i; j < num_variables - 1; j++) {
415  mskEnv->qval[count] = -2 * lap[count + num_variables];
416  /* assert(mskEnv->qval[count]!=0); */
417  mskEnv->qsubi[count] = j;
418  mskEnv->qsubj[count] = i;
419  count++;
420  }
421  }
422 #ifdef DUMP_CONSTRAINTS
423  fprintf(logfile, "Q=[");
424  count = 0;
425  for (i = 0; i < num_variables - 1; i++) {
426  if (i != 0)
427  fprintf(logfile, ";");
428  for (j = 0; j < num_variables - 1; j++) {
429  if (j < i) {
430  fprintf(logfile, "0 ");
431  } else {
432  fprintf(logfile, "%f ", -2 * lap[num_variables + count++]);
433  }
434  }
435  }
436  fprintf(logfile, "]\nQ=Q-diag(diag(Q))+Q'\n");
437 #endif
438  /* Make the mosek environment. */
439  mskEnv->r = MSK_makeenv(&mskEnv->env, NULL, NULL, NULL, NULL);
440 
441  /* Check whether the return code is ok. */
442  if (mskEnv->r == MSK_RES_OK) {
443  /* Directs the log stream to the user
444  specified procedure 'printstr'. */
445  MSK_linkfunctoenvstream(mskEnv->env, MSK_STREAM_LOG, NULL,
446  printstr);
447  }
448 
449  /* Initialize the environment. */
450  mskEnv->r = MSK_initenv(mskEnv->env);
451  if (mskEnv->r == MSK_RES_OK) {
452  /* Make the optimization task. */
453  mskEnv->r =
454  MSK_maketask(mskEnv->env, num_constraints,
455  mskEnv->num_variables, &mskEnv->task);
456 
457  if (mskEnv->r == MSK_RES_OK) {
458  mskEnv->r =
459  MSK_linkfunctotaskstream(mskEnv->task, MSK_STREAM_LOG,
460  NULL, printstr);
461  /* Resize the task. */
462  if (mskEnv->r == MSK_RES_OK)
463  mskEnv->r = MSK_resizetask(mskEnv->task, num_constraints, mskEnv->num_variables, 0, /* no cones!! */
464  /* number of non-zero constraint matrix entries:
465  * each constraint applies to 2 vars
466  */
467  2 * num_constraints,
468  nonzero_lapsize);
469 
470  /* Append the constraints. */
471  if (mskEnv->r == MSK_RES_OK)
472  mskEnv->r = MSK_append(mskEnv->task, 1, num_constraints);
473 
474  /* Append the variables. */
475  if (mskEnv->r == MSK_RES_OK)
476  mskEnv->r =
477  MSK_append(mskEnv->task, 0, mskEnv->num_variables);
478  /* Put variable bounds. */
479  for (j = 0;
480  j < mskEnv->num_variables && mskEnv->r == MSK_RES_OK; j++)
481  mskEnv->r =
482  MSK_putbound(mskEnv->task, 0, j, MSK_BK_RA,
483  -MSK_INFINITY, MSK_INFINITY);
484  for (i = 0; i < num_constraints; i++) {
485  int u = getLeftVarID(cs[i]) - 1;
486  int v = getRightVarID(cs[i]) - 1;
487  double separation = getSeparation(cs[i]);
488  if (u < 0) {
489  mskEnv->r =
490  MSK_putbound(mskEnv->task, 0, v, MSK_BK_RA,
491  -MSK_INFINITY, -separation);
492  assert(mskEnv->r == MSK_RES_OK);
493  } else if (v < 0) {
494  mskEnv->r =
495  MSK_putbound(mskEnv->task, 0, u, MSK_BK_RA,
496  separation, MSK_INFINITY);
497  assert(mskEnv->r == MSK_RES_OK);
498  } else {
499  /* fprintf(stderr,"u=%d,v=%d,sep=%f\n",u,v,separation); */
500  INIT_sub_val(u,v);
501  mskEnv->r =
502  MSK_putavec(mskEnv->task, 1, i, 2, subi, vali);
503  assert(mskEnv->r == MSK_RES_OK);
504  mskEnv->r =
505  MSK_putbound(mskEnv->task, 1, i, MSK_BK_LO,
506  separation, MSK_INFINITY);
507  assert(mskEnv->r == MSK_RES_OK);
508  }
509  }
510  if (mskEnv->r == MSK_RES_OK) {
511  /*
512  * The lower triangular part of the Q
513  * matrix in the objective is specified.
514  */
515  mskEnv->r =
516  MSK_putqobj(mskEnv->task, nonzero_lapsize,
517  mskEnv->qsubi, mskEnv->qsubj,
518  mskEnv->qval);
519  assert(mskEnv->r == MSK_RES_OK);
520  }
521  }
522  }
523  return mskEnv;
524 }
525 
526 /*
527 n: size of b and coords, may be smaller than mskEnv->num_variables if we
528 have dummy vars
529 b: coefficients of linear part of optimisation function
530 coords: optimal y* vector, coord[i] is coordinate of node[i]
531 */
532 void mosek_quad_solve_sep(MosekEnv * mskEnv, int n, float *b,
533  float *coords)
534 {
535  int i, j;
536  assert(n <= mskEnv->num_variables + 1);
537  for (i = 0; i < n - 1 && mskEnv->r == MSK_RES_OK; i++) {
538  mskEnv->r = MSK_putcj(mskEnv->task, i, -2 * b[i + 1]);
539  }
540  if (mskEnv->r == MSK_RES_OK)
541  mskEnv->r = MSK_optimize(mskEnv->task);
542 
543  if (mskEnv->r == MSK_RES_OK) {
544  MSK_getsolutionslice(mskEnv->task,
545  MSK_SOL_ITR,
546  MSK_SOL_ITEM_XX,
547  0, mskEnv->num_variables, mskEnv->xx);
548 
549 #ifdef DUMP_CONSTRAINTS
550  fprintf(logfile, "Primal solution\n");
551 #endif
552  coords[0] = 0;
553  for (j = 1; j <= n; j++) {
554 #ifdef DUMP_CONSTRAINTS
555  fprintf(logfile, "x[%d]: %.2f\n", j, mskEnv->xx[j - 1]);
556 #endif
557  coords[j] = -mskEnv->xx[j - 1];
558  }
559  }
560  fprintf(logfile, "Return code: %d\n", mskEnv->r);
561 }
562 
563 /*
564 please call to clean up
565 */
566 void mosek_delete(MosekEnv * mskEnv)
567 {
568  MSK_deletetask(&mskEnv->task);
569  MSK_deleteenv(&mskEnv->env);
570 
571  if (logfile) {
572  fclose(logfile);
573  logfile = NULL;
574  }
575  free(mskEnv->qval);
576  free(mskEnv->qsubi);
577  free(mskEnv->qsubj);
578  free(mskEnv->xx);
579  free(mskEnv);
580 }
581 #endif /* MOSEK */
assert
#define assert(x)
Definition: cghdr.h:47
quad_prog_vpsc.h
mosek_quad_solve.h
NULL
#define NULL
Definition: logic.h:39
GNEW
#define GNEW(t)
Definition: memory.h:37
node
Agnode_t * node(Agraph_t *g, char *name)
Definition: gv.cpp:103
str
agxbuf * str
Definition: htmlparse.c:85
N_GNEW
#define N_GNEW(n, t)
Definition: agxbuf.c:20
Generated by   doxygen 1.8.5