Graphviz  2.41.20170921.2350
quad_prog_solve.c
Go to the documentation of this file.
1 /* $Id$ $Revision$ */
2 /* vim:set shiftwidth=4 ts=8: */
3 
4 /*************************************************************************
5  * Copyright (c) 2011 AT&T Intellectual Property
6  * All rights reserved. This program and the accompanying materials
7  * are made available under the terms of the Eclipse Public License v1.0
8  * which accompanies this distribution, and is available at
9  * http://www.eclipse.org/legal/epl-v10.html
10  *
11  * Contributors: See CVS logs. Details at http://www.graphviz.org/
12  *************************************************************************/
13 
14 #include "digcola.h"
15 #ifdef DIGCOLA
16 #include <math.h>
17 #include <stdlib.h>
18 #include <time.h>
19 #include <stdio.h>
20 #include <float.h>
21 #include <assert.h>
22 #include "matrix_ops.h"
23 #include "kkutils.h"
24 #include "quad_prog_solver.h"
25 
26 #define quad_prog_tol 1e-2
27 
28 float **unpackMatrix(float *packedMat, int n)
29 {
30  float **mat;
31  int i, j, k;
32 
33  mat = N_GNEW(n, float *);
34  mat[0] = N_GNEW(n * n, float);
35  set_vector_valf(n * n, 0, mat[0]);
36  for (i = 1; i < n; i++) {
37  mat[i] = mat[0] + i * n;
38  }
39 
40  for (i = 0, k = 0; i < n; i++) {
41  for (j = i; j < n; j++, k++) {
42  mat[j][i] = mat[i][j] = packedMat[k];
43  }
44  }
45  return mat;
46 }
47 
48 static void ensureMonotonicOrdering(float *place, int n, int *ordering)
49 {
50  /* ensure that 'ordering' is monotonically increasing by 'place', */
51  /* this also implies that levels are separated in the initial layout */
52  int i, node;
53  float lower_bound = place[ordering[0]];
54  for (i = 1; i < n; i++) {
55  node = ordering[i];
56  if (place[node] < lower_bound) {
57  place[node] = lower_bound;
58  }
59  lower_bound = place[node];
60  }
61 }
62 
63 static void
64 ensureMonotonicOrderingWithGaps(float *place, int n, int *ordering,
65  int *levels, int num_levels,
66  float levels_gap)
67 {
68  /* ensure that levels are separated in the initial layout and that
69  * places are monotonic increasing within layer
70  */
71 
72  int i;
73  int node, level, max_in_level;
74  float lower_bound = (float) -1e9;
75 
76  level = -1;
77  max_in_level = 0;
78  for (i = 0; i < n; i++) {
79  if (i >= max_in_level) {
80  /* we are entering a new level */
81  level++;
82  if (level == num_levels) {
83  /* last_level */
84  max_in_level = n;
85  } else {
86  max_in_level = levels[level];
87  }
88  lower_bound =
89  i > 0 ? place[ordering[i - 1]] + levels_gap : (float) -1e9;
90  quicksort_placef(place, ordering, i, max_in_level - 1);
91  }
92 
93  node = ordering[i];
94  if (place[node] < lower_bound) {
95  place[node] = lower_bound;
96  }
97  }
98 }
99 
100 static void
101 computeHierarchyBoundaries(float *place, int n, int *ordering, int *levels,
102  int num_levels, float *hierarchy_boundaries)
103 {
104  int i;
105  for (i = 0; i < num_levels; i++) {
106  hierarchy_boundaries[i] = place[ordering[levels[i] - 1]];
107  }
108 }
109 
110 
111 int
112 constrained_majorization_new(CMajEnv * e, float *b, float **coords,
113  int cur_axis, int dims, int max_iterations,
114  float *hierarchy_boundaries, float levels_gap)
115 {
116  int n = e->n;
117  float *place = coords[cur_axis];
118  float **lap = e->A;
119  int *ordering = e->ordering;
120  int *levels = e->levels;
121  int num_levels = e->num_levels;
122  int i, j;
123  float new_place_i;
124  boolean converged = FALSE;
125  float upper_bound, lower_bound;
126  int node;
127  int left, right;
128  float cur_place;
129  float des_place_block;
130  float block_deg;
131  float toBlockConnectivity;
132  float *lap_node;
133  int block_len;
134  int first_next_level;
135  int level = -1, max_in_level = 0;
136  float *desired_place;
137  float *prefix_desired_place;
138  float *suffix_desired_place;
139  int *block;
140  int *lev;
141  int counter;
142 
143  if (max_iterations <= 0) {
144  return 0;
145  }
146  if (levels_gap != 0) {
147  return constrained_majorization_new_with_gaps(e, b, coords,
148  cur_axis, dims,
149  max_iterations,
150  hierarchy_boundaries,
151  levels_gap);
152  }
153 
154  /* ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels); */
155  ensureMonotonicOrdering(place, n, ordering);
156  /* it is important that in 'ordering' nodes are always sorted by layers,
157  * and within a layer by 'place'
158  */
159 
160  /* the desired place of each individual node in the current block */
161  desired_place = e->fArray1;
162  /* the desired place of each prefix of current block */
163  prefix_desired_place = e->fArray2;
164  /* the desired place of each suffix of current block */
165  suffix_desired_place = e->fArray3;
166  /* current block (nodes connected by active constraints) */
167  block = e->iArray1;
168 
169  lev = e->iArray2; /* level of each node */
170  for (i = 0; i < n; i++) {
171  if (i >= max_in_level) {
172  /* we are entering a new level */
173  level++;
174  if (level == num_levels) {
175  /* last_level */
176  max_in_level = n;
177  } else {
178  max_in_level = levels[level];
179  }
180  }
181  node = ordering[i];
182  lev[node] = level;
183  }
184 
185  for (counter = 0; counter < max_iterations && !converged; counter++) {
186  converged = TRUE;
187  lower_bound = -1e9; /* no lower bound for first level */
188  for (left = 0; left < n; left = right) {
189  int best_i;
190  double max_movement;
191  double movement;
192  float prefix_des_place, suffix_des_place;
193  /* compute a block 'ordering[left]...ordering[right-1]' of
194  * nodes with the same coordinate:
195  */
196  cur_place = place[ordering[left]];
197  for (right = left + 1; right < n; right++) {
198  if (place[ordering[right]] != cur_place) {
199  break;
200  }
201  }
202 
203  /* compute desired place of nodes in block: */
204  for (i = left; i < right; i++) {
205  node = ordering[i];
206  new_place_i = -b[node];
207  lap_node = lap[node];
208  for (j = 0; j < n; j++) {
209  if (j == node) {
210  continue;
211  }
212  new_place_i += lap_node[j] * place[j];
213  }
214  desired_place[node] = new_place_i / (-lap_node[node]);
215  }
216 
217  /* reorder block by levels, and within levels by "relaxed" desired position */
218  block_len = 0;
219  first_next_level = 0;
220  for (i = left; i < right; i = first_next_level) {
221  level = lev[ordering[i]];
222  if (level == num_levels) {
223  /* last_level */
224  first_next_level = right;
225  } else {
226  first_next_level = MIN(right, levels[level]);
227  }
228 
229  /* First, collect all nodes with desired places smaller than current place */
230  for (j = i; j < first_next_level; j++) {
231  node = ordering[j];
232  if (desired_place[node] < cur_place) {
233  block[block_len++] = node;
234  }
235  }
236  /* Second, collect all nodes with desired places equal to current place */
237  for (j = i; j < first_next_level; j++) {
238  node = ordering[j];
239  if (desired_place[node] == cur_place) {
240  block[block_len++] = node;
241  }
242  }
243  /* Third, collect all nodes with desired places greater than current place */
244  for (j = i; j < first_next_level; j++) {
245  node = ordering[j];
246  if (desired_place[node] > cur_place) {
247  block[block_len++] = node;
248  }
249  }
250  }
251 
252  /* loop through block and compute desired places of its prefixes */
253  des_place_block = 0;
254  block_deg = 0;
255  for (i = 0; i < block_len; i++) {
256  node = block[i];
257  toBlockConnectivity = 0;
258  lap_node = lap[node];
259  for (j = 0; j < i; j++) {
260  toBlockConnectivity -= lap_node[block[j]];
261  }
262  toBlockConnectivity *= 2;
263  /* update block stats */
264  des_place_block =
265  (block_deg * des_place_block +
266  (-lap_node[node]) * desired_place[node] +
267  toBlockConnectivity * cur_place) / (block_deg -
268  lap_node[node] +
269  toBlockConnectivity);
270  prefix_desired_place[i] = des_place_block;
271  block_deg += toBlockConnectivity - lap_node[node];
272  }
273 
274  /* loop through block and compute desired places of its suffixes */
275  des_place_block = 0;
276  block_deg = 0;
277  for (i = block_len - 1; i >= 0; i--) {
278  node = block[i];
279  toBlockConnectivity = 0;
280  lap_node = lap[node];
281  for (j = i + 1; j < block_len; j++) {
282  toBlockConnectivity -= lap_node[block[j]];
283  }
284  toBlockConnectivity *= 2;
285  /* update block stats */
286  des_place_block =
287  (block_deg * des_place_block +
288  (-lap_node[node]) * desired_place[node] +
289  toBlockConnectivity * cur_place) / (block_deg -
290  lap_node[node] +
291  toBlockConnectivity);
292  suffix_desired_place[i] = des_place_block;
293  block_deg += toBlockConnectivity - lap_node[node];
294  }
295 
296 
297  /* now, find best place to split block */
298  best_i = -1;
299  max_movement = 0;
300  for (i = 0; i < block_len; i++) {
301  suffix_des_place = suffix_desired_place[i];
302  prefix_des_place =
303  i > 0 ? prefix_desired_place[i - 1] : suffix_des_place;
304  /* limit moves to ensure that the prefix is placed before the suffix */
305  if (suffix_des_place < prefix_des_place) {
306  if (suffix_des_place < cur_place) {
307  if (prefix_des_place > cur_place) {
308  prefix_des_place = cur_place;
309  }
310  suffix_des_place = prefix_des_place;
311  } else if (prefix_des_place > cur_place) {
312  prefix_des_place = suffix_des_place;
313  }
314  }
315  movement =
316  (block_len - i) * fabs(suffix_des_place - cur_place) +
317  i * fabs(prefix_des_place - cur_place);
318  if (movement > max_movement) {
319  max_movement = movement;
320  best_i = i;
321  }
322  }
323  /* Actually move prefix and suffix */
324  if (best_i >= 0) {
325  suffix_des_place = suffix_desired_place[best_i];
326  prefix_des_place =
327  best_i >
328  0 ? prefix_desired_place[best_i -
329  1] : suffix_des_place;
330 
331  /* compute right border of feasible move */
332  if (right >= n) {
333  /* no nodes after current block */
334  upper_bound = 1e9;
335  } else {
336  upper_bound = place[ordering[right]];
337  }
338  suffix_des_place = MIN(suffix_des_place, upper_bound);
339  prefix_des_place = MAX(prefix_des_place, lower_bound);
340 
341  /* limit moves to ensure that the prefix is placed before the suffix */
342  if (suffix_des_place < prefix_des_place) {
343  if (suffix_des_place < cur_place) {
344  if (prefix_des_place > cur_place) {
345  prefix_des_place = cur_place;
346  }
347  suffix_des_place = prefix_des_place;
348  } else if (prefix_des_place > cur_place) {
349  prefix_des_place = suffix_des_place;
350  }
351  }
352 
353  /* move prefix: */
354  for (i = 0; i < best_i; i++) {
355  place[block[i]] = prefix_des_place;
356  }
357  /* move suffix: */
358  for (i = best_i; i < block_len; i++) {
359  place[block[i]] = suffix_des_place;
360  }
361 
362  lower_bound = suffix_des_place; /* lower bound for next block */
363 
364  /* reorder 'ordering' to reflect change of places
365  * Note that it is enough to reorder the level where
366  * the split was done
367  */
368 #if 0
369  int max_in_level, min_in_level;
370 
371  level = lev[best_i];
372  if (level == num_levels) {
373  /* last_level */
374  max_in_level = MIN(right, n);
375  } else {
376  max_in_level = MIN(right, levels[level]);
377  }
378  if (level == 0) {
379  /* first level */
380  min_in_level = MAX(left, 0);
381  } else {
382  min_in_level = MAX(left, levels[level - 1]);
383  }
384 #endif
385  for (i = left; i < right; i++) {
386  ordering[i] = block[i - left];
387  }
388  converged = converged
389  && fabs(prefix_des_place - cur_place) < quad_prog_tol
390  && fabs(suffix_des_place - cur_place) < quad_prog_tol;
391 
392  } else {
393  /* no movement */
394  lower_bound = cur_place; /* lower bound for next block */
395  }
396 
397  }
398  }
399 
400  computeHierarchyBoundaries(place, n, ordering, levels, num_levels,
401  hierarchy_boundaries);
402 
403  return counter;
404 }
405 
406 #ifdef IPSEPCOLA
407 static float *place;
408 static int compare_incr(const void *a, const void *b)
409 {
410  if (place[*(int *) a] > place[*(int *) b]) {
411  return 1;
412  } else if (place[*(int *) a] < place[*(int *) b]) {
413  return -1;
414  }
415  return 0;
416 }
417 
418 /*
419 While not converged: move everything towards the optimum, then satisfy constraints with as little displacement as possible.
420 Returns number of iterations before convergence.
421 */
422 int constrained_majorization_gradient_projection(CMajEnv * e,
423  float *b, float **coords,
424  int ndims, int cur_axis,
425  int max_iterations,
426  float
427  *hierarchy_boundaries,
428  float levels_gap)
429 {
430 
431  int i, j, counter;
432  int *ordering = e->ordering;
433  int *levels = e->levels;
434  int num_levels = e->num_levels;
435  boolean converged = FALSE;
436  float *g = e->fArray1;
437  float *old_place = e->fArray2;
438  float *d = e->fArray4;
439  float test = 0, tmptest = 0;
440  float beta;
441 
442  if (max_iterations == 0)
443  return 0;
444 
445  place = coords[cur_axis];
446 #ifdef CONMAJ_LOGGING
447  double prev_stress = 0;
448  static int call_no = 0;
449  for (i = 0; i < e->n; i++) {
450  prev_stress += 2 * b[i] * place[i];
451  for (j = 0; j < e->n; j++) {
452  prev_stress -= e->A[i][j] * place[j] * place[i];
453  }
454  }
455  FILE *logfile = fopen("constrained_majorization_log", "a");
456 
457  fprintf(logfile, "grad proj %d: stress=%f\n", call_no, prev_stress);
458 #endif
459  for (counter = 0; counter < max_iterations && !converged; counter++) {
460  float alpha;
461  float numerator = 0, denominator = 0, r;
462  converged = TRUE;
463  /* find steepest descent direction */
464  for (i = 0; i < e->n; i++) {
465  old_place[i] = place[i];
466  g[i] = 2 * b[i];
467  for (j = 0; j < e->n; j++) {
468  g[i] -= 2 * e->A[i][j] * place[j];
469  }
470  }
471  for (i = 0; i < e->n; i++) {
472  numerator += g[i] * g[i];
473  r = 0;
474  for (j = 0; j < e->n; j++) {
475  r += 2 * e->A[i][j] * g[j];
476  }
477  denominator -= r * g[i];
478  }
479  alpha = numerator / denominator;
480  for (i = 0; i < e->n; i++) {
481  if (alpha > 0 && alpha < 1000) {
482  place[i] -= alpha * g[i];
483  }
484  }
485  if (num_levels)
486  qsort((int *) ordering, (size_t) levels[0], sizeof(int),
487  compare_incr);
488  /* project to constraint boundary */
489  for (i = 0; i < num_levels; i++) {
490  int endOfLevel = i == num_levels - 1 ? e->n : levels[i + 1];
491  int ui, li, u, l;
492 
493  /* ensure monotic increase in position within levels */
494  qsort((int *) ordering + levels[i],
495  (size_t) endOfLevel - levels[i], sizeof(int),
496  compare_incr);
497  /* If there are overlapping levels find offending nodes and place at average position */
498  ui = levels[i]; li = ui - 1;
499  l = ordering[li--]; u = ordering[ui++];
500  if (place[l] + levels_gap > place[u]) {
501  float sum =
502  place[l] + place[u] - levels_gap * (e->lev[l] +
503  e->lev[u]), w = 2;
504  float avgPos = sum / w;
505  float pos;
506  boolean finished;
507  do {
508  finished = TRUE;
509  if (ui < endOfLevel) {
510  u = ordering[ui];
511  pos = place[u] - levels_gap * e->lev[u];
512  if (pos < avgPos) {
513  ui++;
514  w++;
515  sum += pos;
516  avgPos = sum / w;
517  finished = FALSE;
518  }
519  }
520 
521  if (li >= 0) {
522  l = ordering[li];
523  pos = place[l] - levels_gap * e->lev[l];
524  if (pos > avgPos) {
525  li--;
526  w++;
527  sum += pos;
528  avgPos = sum / w;
529  finished = FALSE;
530  }
531  }
532  } while (!finished);
533  for (j = li + 1; j < ui; j++) {
534  place[ordering[j]] =
535  avgPos + levels_gap * e->lev[ordering[j]];
536  }
537  }
538  }
539  /* set place to the intersection of old_place-g and boundary and compute d, the vector from intersection pnt to projection pnt */
540  for (i = 0; i < e->n; i++) {
541  d[i] = place[i] - old_place[i];
542  }
543  /* now compute beta */
544  numerator = 0, denominator = 0;
545  for (i = 0; i < e->n; i++) {
546  numerator += g[i] * d[i];
547  r = 0;
548  for (j = 0; j < e->n; j++) {
549  r += 2 * e->A[i][j] * d[j];
550  }
551  denominator += r * d[i];
552  }
553  beta = numerator / denominator;
554 
555  for (i = 0; i < e->n; i++) {
556  if (beta > 0 && beta < 1.0) {
557  place[i] = old_place[i] + beta * d[i];
558  }
559  tmptest = fabs(place[i] - old_place[i]);
560  if (test < tmptest)
561  test = tmptest;
562  }
563  computeHierarchyBoundaries(place, e->n, ordering, levels,
564  num_levels, hierarchy_boundaries);
565 #if 0
566  if (num_levels)
567  qsort((int *) ordering, (size_t) levels[0], sizeof(int),
568  compare_incr);
569  for (i = 0; i < num_levels; i++) {
570  int endOfLevel = i == num_levels - 1 ? e->n : levels[i + 1];
571  /* ensure monotic increase in position within levels */
572  qsort((int *) ordering + levels[i],
573  (size_t) endOfLevel - levels[i], sizeof(int),
574  compare_incr);
575  /* If there are overlapping levels find offending nodes and place at average position */
576  int l = ordering[levels[i] - 1], u = ordering[levels[i]];
577  /* assert(place[l]+levels_gap<=place[u]+0.00001); */
578  }
579 #endif
580 #ifdef CONMAJ_LOGGING
581  double stress = 0;
582  for (i = 0; i < e->n; i++) {
583  stress += 2 * b[i] * place[i];
584  for (j = 0; j < e->n; j++) {
585  stress -= e->A[i][j] * place[j] * place[i];
586  }
587  }
588  fprintf(logfile, "%d: stress=%f, test=%f, %s\n", call_no, stress,
589  test, (stress >= prev_stress) ? "No Improvement" : "");
590  prev_stress = stress;
591 #endif
592  if (test > quad_prog_tol) {
593  converged = FALSE;
594  }
595  }
596 #ifdef CONMAJ_LOGGING
597  call_no++;
598  fclose(logfile);
599 #endif
600  return counter;
601 }
602 #endif
603 
604 int
605 constrained_majorization_new_with_gaps(CMajEnv * e, float *b,
606  float **coords, int ndims,
607  int cur_axis, int max_iterations,
608  float *hierarchy_boundaries,
609  float levels_gap)
610 {
611  float *place = coords[cur_axis];
612  int i, j;
613  int n = e->n;
614  float **lap = e->A;
615  int *ordering = e->ordering;
616  int *levels = e->levels;
617  int num_levels = e->num_levels;
618  float new_place_i;
619  boolean converged = FALSE;
620  float upper_bound, lower_bound;
621  int node;
622  int left, right;
623  float cur_place;
624  float des_place_block;
625  float block_deg;
626  float toBlockConnectivity;
627  float *lap_node;
628  float *desired_place;
629  float *prefix_desired_place;
630  float *suffix_desired_place;
631  int *block;
632  int block_len;
633  int first_next_level;
634  int *lev;
635  int level = -1, max_in_level = 0;
636  int counter;
637  float *gap;
638  float total_gap, target_place;
639 
640  if (max_iterations <= 0) {
641  return 0;
642  }
643 
644  ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels,
645  levels_gap);
646  /* it is important that in 'ordering' nodes are always sorted by layers,
647  * and within a layer by 'place'
648  */
649 
650  /* the desired place of each individual node in the current block */
651  desired_place = e->fArray1;
652  /* the desired place of each prefix of current block */
653  prefix_desired_place = e->fArray2;
654  /* the desired place of each suffix of current block */
655  suffix_desired_place = e->fArray3;
656  /* current block (nodes connected by active constraints) */
657  block = e->iArray1;
658 
659  lev = e->iArray2; /* level of each node */
660  for (i = 0; i < n; i++) {
661  if (i >= max_in_level) {
662  /* we are entering a new level */
663  level++;
664  if (level == num_levels) {
665  /* last_level */
666  max_in_level = n;
667  } else {
668  max_in_level = levels[level];
669  }
670  }
671  node = ordering[i];
672  lev[node] = level;
673  }
674 
675  /* displacement of block's nodes from block's reference point */
676  gap = e->fArray4;
677 
678  for (counter = 0; counter < max_iterations && !converged; counter++) {
679  converged = TRUE;
680  lower_bound = -1e9; /* no lower bound for first level */
681  for (left = 0; left < n; left = right) {
682  int best_i;
683  double max_movement;
684  double movement;
685  float prefix_des_place, suffix_des_place;
686  /* compute a block 'ordering[left]...ordering[right-1]' of
687  * nodes connected with active constraints
688  */
689  cur_place = place[ordering[left]];
690  total_gap = 0;
691  target_place = cur_place;
692  gap[ordering[left]] = 0;
693  for (right = left + 1; right < n; right++) {
694  if (lev[right] > lev[right - 1]) {
695  /* we are entering a new level */
696  target_place += levels_gap; /* note that 'levels_gap' may be negative */
697  total_gap += levels_gap;
698  }
699  node = ordering[right];
700 #if 0
701  if (place[node] != target_place)
702 #endif
703  /* not sure if this is better than 'place[node]!=target_place' */
704  if (fabs(place[node] - target_place) > 1e-9) {
705  break;
706  }
707  gap[node] = place[node] - cur_place;
708  }
709 
710  /* compute desired place of block's reference point according
711  * to each node in the block:
712  */
713  for (i = left; i < right; i++) {
714  node = ordering[i];
715  new_place_i = -b[node];
716  lap_node = lap[node];
717  for (j = 0; j < n; j++) {
718  if (j == node) {
719  continue;
720  }
721  new_place_i += lap_node[j] * place[j];
722  }
723  desired_place[node] =
724  new_place_i / (-lap_node[node]) - gap[node];
725  }
726 
727  /* reorder block by levels, and within levels
728  * by "relaxed" desired position
729  */
730  block_len = 0;
731  first_next_level = 0;
732  for (i = left; i < right; i = first_next_level) {
733  level = lev[ordering[i]];
734  if (level == num_levels) {
735  /* last_level */
736  first_next_level = right;
737  } else {
738  first_next_level = MIN(right, levels[level]);
739  }
740 
741  /* First, collect all nodes with desired places smaller
742  * than current place
743  */
744  for (j = i; j < first_next_level; j++) {
745  node = ordering[j];
746  if (desired_place[node] < cur_place) {
747  block[block_len++] = node;
748  }
749  }
750  /* Second, collect all nodes with desired places equal
751  * to current place
752  */
753  for (j = i; j < first_next_level; j++) {
754  node = ordering[j];
755  if (desired_place[node] == cur_place) {
756  block[block_len++] = node;
757  }
758  }
759  /* Third, collect all nodes with desired places greater
760  * than current place
761  */
762  for (j = i; j < first_next_level; j++) {
763  node = ordering[j];
764  if (desired_place[node] > cur_place) {
765  block[block_len++] = node;
766  }
767  }
768  }
769 
770  /* loop through block and compute desired places of its prefixes */
771  des_place_block = 0;
772  block_deg = 0;
773  for (i = 0; i < block_len; i++) {
774  node = block[i];
775  toBlockConnectivity = 0;
776  lap_node = lap[node];
777  for (j = 0; j < i; j++) {
778  toBlockConnectivity -= lap_node[block[j]];
779  }
780  toBlockConnectivity *= 2;
781  /* update block stats */
782  des_place_block =
783  (block_deg * des_place_block +
784  (-lap_node[node]) * desired_place[node] +
785  toBlockConnectivity * cur_place) / (block_deg -
786  lap_node[node] +
787  toBlockConnectivity);
788  prefix_desired_place[i] = des_place_block;
789  block_deg += toBlockConnectivity - lap_node[node];
790  }
791 
792  if (block_len == n) {
793  /* fix is needed since denominator was 0 in this case */
794  prefix_desired_place[n - 1] = cur_place; /* a "neutral" value */
795  }
796 
797  /* loop through block and compute desired places of its suffixes */
798  des_place_block = 0;
799  block_deg = 0;
800  for (i = block_len - 1; i >= 0; i--) {
801  node = block[i];
802  toBlockConnectivity = 0;
803  lap_node = lap[node];
804  for (j = i + 1; j < block_len; j++) {
805  toBlockConnectivity -= lap_node[block[j]];
806  }
807  toBlockConnectivity *= 2;
808  /* update block stats */
809  des_place_block =
810  (block_deg * des_place_block +
811  (-lap_node[node]) * desired_place[node] +
812  toBlockConnectivity * cur_place) / (block_deg -
813  lap_node[node] +
814  toBlockConnectivity);
815  suffix_desired_place[i] = des_place_block;
816  block_deg += toBlockConnectivity - lap_node[node];
817  }
818 
819  if (block_len == n) {
820  /* fix is needed since denominator was 0 in this case */
821  suffix_desired_place[0] = cur_place; /* a "neutral" value? */
822  }
823 
824 
825  /* now, find best place to split block */
826  best_i = -1;
827  max_movement = 0;
828  for (i = 0; i < block_len; i++) {
829  suffix_des_place = suffix_desired_place[i];
830  prefix_des_place =
831  i > 0 ? prefix_desired_place[i - 1] : suffix_des_place;
832  /* limit moves to ensure that the prefix is placed before the suffix */
833  if (suffix_des_place < prefix_des_place) {
834  if (suffix_des_place < cur_place) {
835  if (prefix_des_place > cur_place) {
836  prefix_des_place = cur_place;
837  }
838  suffix_des_place = prefix_des_place;
839  } else if (prefix_des_place > cur_place) {
840  prefix_des_place = suffix_des_place;
841  }
842  }
843  movement =
844  (block_len - i) * fabs(suffix_des_place - cur_place) +
845  i * fabs(prefix_des_place - cur_place);
846  if (movement > max_movement) {
847  max_movement = movement;
848  best_i = i;
849  }
850  }
851  /* Actually move prefix and suffix */
852  if (best_i >= 0) {
853  suffix_des_place = suffix_desired_place[best_i];
854  prefix_des_place =
855  best_i >
856  0 ? prefix_desired_place[best_i -
857  1] : suffix_des_place;
858 
859  /* compute right border of feasible move */
860  if (right >= n) {
861  /* no nodes after current block */
862  upper_bound = 1e9;
863  } else {
864  /* notice that we have to deduct 'gap[block[block_len-1]]'
865  * since all computations should be relative to
866  * the block's reference point
867  */
868  if (lev[ordering[right]] > lev[ordering[right - 1]]) {
869  upper_bound =
870  place[ordering[right]] - levels_gap -
871  gap[block[block_len - 1]];
872  } else {
873  upper_bound =
874  place[ordering[right]] -
875  gap[block[block_len - 1]];
876  }
877  }
878  suffix_des_place = MIN(suffix_des_place, upper_bound);
879  prefix_des_place = MAX(prefix_des_place, lower_bound);
880 
881  /* limit moves to ensure that the prefix is placed before the suffix */
882  if (suffix_des_place < prefix_des_place) {
883  if (suffix_des_place < cur_place) {
884  if (prefix_des_place > cur_place) {
885  prefix_des_place = cur_place;
886  }
887  suffix_des_place = prefix_des_place;
888  } else if (prefix_des_place > cur_place) {
889  prefix_des_place = suffix_des_place;
890  }
891  }
892 
893  /* move prefix: */
894  for (i = 0; i < best_i; i++) {
895  place[block[i]] = prefix_des_place + gap[block[i]];
896  }
897  /* move suffix: */
898  for (i = best_i; i < block_len; i++) {
899  place[block[i]] = suffix_des_place + gap[block[i]];
900  }
901 
902 
903  /* compute lower bound for next block */
904  if (right < n
905  && lev[ordering[right]] > lev[ordering[right - 1]]) {
906  lower_bound = place[block[block_len - 1]] + levels_gap;
907  } else {
908  lower_bound = place[block[block_len - 1]];
909  }
910 
911 
912  /* reorder 'ordering' to reflect change of places
913  * Note that it is enough to reorder the level where
914  * the split was done
915  */
916 #if 0
917  int max_in_level, min_in_level;
918 
919  level = lev[best_i];
920  if (level == num_levels) {
921  /* last_level */
922  max_in_level = MIN(right, n);
923  } else {
924  max_in_level = MIN(right, levels[level]);
925  }
926  if (level == 0) {
927  /* first level */
928  min_in_level = MAX(left, 0);
929  } else {
930  min_in_level = MAX(left, levels[level - 1]);
931  }
932 #endif
933  for (i = left; i < right; i++) {
934  ordering[i] = block[i - left];
935  }
936  converged = converged
937  && fabs(prefix_des_place - cur_place) < quad_prog_tol
938  && fabs(suffix_des_place - cur_place) < quad_prog_tol;
939 
940 
941  } else {
942  /* no movement */
943  /* compute lower bound for next block */
944  if (right < n
945  && lev[ordering[right]] > lev[ordering[right - 1]]) {
946  lower_bound = place[block[block_len - 1]] + levels_gap;
947  } else {
948  lower_bound = place[block[block_len - 1]];
949  }
950  }
951  }
952  orthog1f(n, place); /* for numerical stability, keep ||place|| small */
953  computeHierarchyBoundaries(place, n, ordering, levels, num_levels,
954  hierarchy_boundaries);
955  }
956 
957  return counter;
958 }
959 
960 void deleteCMajEnv(CMajEnv * e)
961 {
962  free(e->A[0]);
963  free(e->A);
964  free(e->lev);
965  free(e->fArray1);
966  free(e->fArray2);
967  free(e->fArray3);
968  free(e->fArray4);
969  free(e->iArray1);
970  free(e->iArray2);
971  free(e->iArray3);
972  free(e->iArray4);
973  free(e);
974 }
975 
976 CMajEnv *initConstrainedMajorization(float *packedMat, int n,
977  int *ordering, int *levels,
978  int num_levels)
979 {
980  int i, level = -1, start_of_level_above = 0;
981  CMajEnv *e = GNEW(CMajEnv);
982  e->A = NULL;
983  e->n = n;
984  e->ordering = ordering;
985  e->levels = levels;
986  e->num_levels = num_levels;
987  e->A = unpackMatrix(packedMat, n);
988  e->lev = N_GNEW(n, int);
989  for (i = 0; i < e->n; i++) {
990  if (i >= start_of_level_above) {
991  level++;
992  start_of_level_above =
993  (level == num_levels) ? e->n : levels[level];
994  }
995  e->lev[ordering[i]] = level;
996  }
997  e->fArray1 = N_GNEW(n, float);
998  e->fArray2 = N_GNEW(n, float);
999  e->fArray3 = N_GNEW(n, float);
1000  e->fArray4 = N_GNEW(n, float);
1001  e->iArray1 = N_GNEW(n, int);
1002  e->iArray2 = N_GNEW(n, int);
1003  e->iArray3 = N_GNEW(n, int);
1004  e->iArray4 = N_GNEW(n, int);
1005  return e;
1006 }
1007 #endif /* DIGCOLA */
#define MAX(a, b)
Definition: agerror.c:17
#define MIN(a, b)
Definition: arith.h:35
void set_vector_valf(int n, float val, float *result)
Definition: matrix_ops.c:685
Definition: block.h:30
#define NULL
Definition: logic.h:39
#define GNEW(t)
Definition: memory.h:37
#define right(i)
Definition: closest.c:87
Agnode_t * node(Agraph_t *g, char *name)
Definition: gv.cpp:103
#define alpha
Definition: shapes.c:3902
#define left
Definition: dthdr.h:16
void orthog1f(int n, float *vec)
Definition: matrix_ops.c:544
#define N_GNEW(n, t)
Definition: agxbuf.c:20
#define FALSE
Definition: cgraph.h:35
void quicksort_placef(float *place, int *ordering, int first, int last)
Definition: kkutils.c:194
#define TRUE
Definition: cgraph.h:38