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
28float **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
48static 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
63static void
64ensureMonotonicOrderingWithGaps(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
100static void
101computeHierarchyBoundaries(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
111int
112constrained_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
407static float *place;
408static 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/*
419While not converged: move everything towards the optimum, then satisfy constraints with as little displacement as possible.
420Returns number of iterations before convergence.
421*/
422int 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
604int
605constrained_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
960void 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
976CMajEnv *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 */
1008