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 */ |
1008 | |