1 /*
2 * Copyright 2007 Sandia Corporation. Under the terms of Contract
3 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
4 * certain rights in this software.
5 *
6 * All rights reserved.
7 *
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted provided that the following conditions are
10 * met:
11 *
12 * * Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * * Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * * Neither the name of Sandia National Laboratories nor the names of
18 * its contributors may be used to endorse or promote products derived from
19 * this software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
27 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33 // This file contains the member definitions of the master class
34
35
36 #include <map>
37 #include <vector>
38 #include <cmath>
39
40 using namespace std;
41
42 #include "drl_graph.h"
43 #include "igraph_random.h"
44 #include "igraph_interface.h"
45 #include "igraph_progress.h"
46 #include "igraph_interrupt_internal.h"
47 #ifdef MUSE_MPI
48 #include <mpi.h>
49 #endif
50
51 namespace drl {
52
53 // constructor -- initializes the schedule variables (as in
54 // graph constructor)
55
56 // graph::graph ( int proc_id, int tot_procs, char *int_file )
57 // {
58
59 // // MPI parameters
60 // myid = proc_id;
61 // num_procs = tot_procs;
62
63 // // initial annealing parameters
64 // STAGE = 0;
65 // iterations = 0;
66 // temperature = 2000;
67 // attraction = 10;
68 // damping_mult = 1.0;
69 // min_edges = 20;
70 // first_add = fine_first_add = true;
71 // fineDensity = false;
72
73 // // Brian's original Vx schedule
74 // liquid.iterations = 200;
75 // liquid.temperature = 2000;
76 // liquid.attraction = 2;
77 // liquid.damping_mult = 1.0;
78 // liquid.time_elapsed = 0;
79
80 // expansion.iterations = 200;
81 // expansion.temperature = 2000;
82 // expansion.attraction = 10;
83 // expansion.damping_mult = 1.0;
84 // expansion.time_elapsed = 0;
85
86 // cooldown.iterations = 200;
87 // cooldown.temperature = 2000;
88 // cooldown.attraction = 1;
89 // cooldown.damping_mult = .1;
90 // cooldown.time_elapsed = 0;
91
92 // crunch.iterations = 50;
93 // crunch.temperature = 250;
94 // crunch.attraction = 1;
95 // crunch. damping_mult = .25;
96 // crunch.time_elapsed = 0;
97
98 // simmer.iterations = 100;
99 // simmer.temperature = 250;
100 // simmer.attraction = .5;
101 // simmer.damping_mult = 0.0;
102 // simmer.time_elapsed = 0;
103
104 // // scan .int file for node info
105 // scan_int ( int_file );
106
107 // // populate node positions and ids
108 // positions.reserve ( num_nodes );
109 // map < int, int >::iterator cat_iter;
110 // for ( cat_iter = id_catalog.begin();
111 // cat_iter != id_catalog.end();
112 // cat_iter++ )
113 // positions.push_back ( Node( cat_iter->first ) );
114
115 // /*
116 // // output positions .ids for debugging
117 // for ( int id = 0; id < num_nodes; id++ )
118 // cout << positions[id].id << endl;
119 // */
120
121 // // read .int file for graph info
122 // read_int ( int_file );
123
124 // // initialize density server
125 // density_server.Init();
126
127 // }
128
graph(const igraph_t * igraph,const igraph_layout_drl_options_t * options,const igraph_vector_t * weights)129 graph::graph(const igraph_t *igraph,
130 const igraph_layout_drl_options_t *options,
131 const igraph_vector_t *weights) {
132 myid = 0;
133 num_procs = 1;
134
135 STAGE = 0;
136 iterations = options->init_iterations;
137 temperature = options->init_temperature;
138 attraction = options->init_attraction;
139 damping_mult = options->init_damping_mult;
140 min_edges = 20;
141 first_add = fine_first_add = true;
142 fineDensity = false;
143
144 // Brian's original Vx schedule
145 liquid.iterations = options->liquid_iterations;
146 liquid.temperature = options->liquid_temperature;
147 liquid.attraction = options->liquid_attraction;
148 liquid.damping_mult = options->liquid_damping_mult;
149 liquid.time_elapsed = 0;
150
151 expansion.iterations = options->expansion_iterations;
152 expansion.temperature = options->expansion_temperature;
153 expansion.attraction = options->expansion_attraction;
154 expansion.damping_mult = options->expansion_damping_mult;
155 expansion.time_elapsed = 0;
156
157 cooldown.iterations = options->cooldown_iterations;
158 cooldown.temperature = options->cooldown_temperature;
159 cooldown.attraction = options->cooldown_attraction;
160 cooldown.damping_mult = options->cooldown_damping_mult;
161 cooldown.time_elapsed = 0;
162
163 crunch.iterations = options->crunch_iterations;
164 crunch.temperature = options->crunch_temperature;
165 crunch.attraction = options->crunch_attraction;
166 crunch.damping_mult = options->crunch_damping_mult;
167 crunch.time_elapsed = 0;
168
169 simmer.iterations = options->simmer_iterations;
170 simmer.temperature = options->simmer_temperature;
171 simmer.attraction = options->simmer_attraction;
172 simmer.damping_mult = options->simmer_damping_mult;
173 simmer.time_elapsed = 0;
174
175 // scan .int file for node info
176 highest_sim = 1.0;
177 num_nodes = igraph_vcount(igraph);
178 long int no_of_edges = igraph_ecount(igraph);
179 for (long int i = 0; i < num_nodes; i++) {
180 id_catalog[i] = 1;
181 }
182 map< int, int>::iterator cat_iter;
183 for ( cat_iter = id_catalog.begin();
184 cat_iter != id_catalog.end(); cat_iter++) {
185 cat_iter->second = cat_iter->first;
186 }
187
188 // populate node positions and ids
189 positions.reserve ( num_nodes );
190 for ( cat_iter = id_catalog.begin();
191 cat_iter != id_catalog.end();
192 cat_iter++ ) {
193 positions.push_back ( Node( cat_iter->first ) );
194 }
195
196 // read .int file for graph info
197 long int node_1, node_2;
198 double weight;
199 for (long int i = 0; i < no_of_edges; i++) {
200 node_1 = IGRAPH_FROM(igraph, i);
201 node_2 = IGRAPH_TO(igraph, i);
202 weight = weights ? VECTOR(*weights)[i] : 1.0 ;
203 (neighbors[id_catalog[node_1]])[id_catalog[node_2]] = weight;
204 (neighbors[id_catalog[node_2]])[id_catalog[node_1]] = weight;
205 }
206
207 // initialize density server
208 density_server.Init();
209
210 }
211
212 // The following subroutine scans the .int file for the following
213 // information: number nodes, node ids, and highest similarity. The
214 // corresponding graph globals are populated: num_nodes, id_catalog,
215 // and highest_sim.
216
217 // void graph::scan_int ( char *filename )
218 // {
219
220 // cout << "Proc. " << myid << " scanning .int file ..." << endl;
221
222 // // Open (sim) File
223 // ifstream fp ( filename );
224 // if ( !fp )
225 // {
226 // cout << "Error: could not open " << filename << ". Program terminated." << endl;
227 // #ifdef MUSE_MPI
228 // MPI_Abort ( MPI_COMM_WORLD, 1 );
229 // #else
230 // exit (1);
231 // #endif
232 // }
233
234 // // Read file, parse, and add into data structure
235 // int id1, id2;
236 // float edge_weight;
237 // highest_sim = -1.0;
238 // while ( !fp.eof () )
239 // {
240 // fp >> id1 >> id2 >> edge_weight;
241
242 // // ignore negative weights!
243 // if ( edge_weight <= 0 )
244 // {
245 // cout << "Error: found negative edge weight in " << filename << ". Program stopped." << endl;
246 // #ifdef MUSE_MPI
247 // MPI_Abort ( MPI_COMM_WORLD, 1 );
248 // #else
249 // exit (1);
250 // #endif
251 // }
252
253 // if ( highest_sim < edge_weight )
254 // highest_sim = edge_weight;
255
256 // id_catalog[id1] = 1;
257 // id_catalog[id2] = 1;
258 // }
259
260 // fp.close();
261
262 // if ( id_catalog.size() == 0 )
263 // {
264 // cout << "Error: Proc. " << myid << ": " << filename << " is empty. Program terminated." << endl;
265 // #ifdef MUSE_MPI
266 // MPI_Abort ( MPI_COMM_WORLD, 1 );
267 // #else
268 // exit (1);
269 // #endif
270 // }
271
272 // // label nodes with sequential integers starting at 0
273 // map< int, int>::iterator cat_iter;
274 // int id_label;
275 // for ( cat_iter = id_catalog.begin(), id_label = 0;
276 // cat_iter != id_catalog.end(); cat_iter++, id_label++ )
277 // cat_iter->second = id_label;
278
279 // /*
280 // // output id_catalog for debugging:
281 // for ( cat_iter = id_catalog.begin();
282 // cat_iter != id_catalog.end();
283 // cat_iter++ )
284 // cout << cat_iter->first << "\t" << cat_iter->second << endl;
285 // */
286
287 // num_nodes = id_catalog.size();
288 // }
289
290 // read in .parms file, if present
291
292 /*
293 void graph::read_parms ( char *parms_file )
294 {
295
296 // read from .parms file
297 ifstream parms_in ( parms_file );
298 if ( !parms_in )
299 {
300 cout << "Error: could not open .parms file! Program stopped." << endl;
301 #ifdef MUSE_MPI
302 MPI_Abort ( MPI_COMM_WORLD, 1 );
303 #else
304 exit (1);
305 #endif
306 }
307
308 cout << "Processor " << myid << " reading .parms file." << endl;
309
310 // read in stage parameters
311 string parm_label; // this is ignored in the .parms file
312
313 // initial parameters
314 parms_in >> parm_label >> iterations;
315 parms_in >> parm_label >> temperature;
316 parms_in >> parm_label >> attraction;
317 parms_in >> parm_label >> damping_mult;
318
319 // liquid stage
320 parms_in >> parm_label >> liquid.iterations;
321 parms_in >> parm_label >> liquid.temperature;
322 parms_in >> parm_label >> liquid.attraction;
323 parms_in >> parm_label >> liquid.damping_mult;
324
325 // expansion stage
326 parms_in >> parm_label >> expansion.iterations;
327 parms_in >> parm_label >> expansion.temperature;
328 parms_in >> parm_label >> expansion.attraction;
329 parms_in >> parm_label >> expansion.damping_mult;
330
331 // cooldown stage
332 parms_in >> parm_label >> cooldown.iterations;
333 parms_in >> parm_label >> cooldown.temperature;
334 parms_in >> parm_label >> cooldown.attraction;
335 parms_in >> parm_label >> cooldown.damping_mult;
336
337 // crunch stage
338 parms_in >> parm_label >> crunch.iterations;
339 parms_in >> parm_label >> crunch.temperature;
340 parms_in >> parm_label >> crunch.attraction;
341 parms_in >> parm_label >> crunch.damping_mult;
342
343 // simmer stage
344 parms_in >> parm_label >> simmer.iterations;
345 parms_in >> parm_label >> simmer.temperature;
346 parms_in >> parm_label >> simmer.attraction;
347 parms_in >> parm_label >> simmer.damping_mult;
348
349 parms_in.close();
350
351 // print out parameters for double checking
352 if ( myid == 0 )
353 {
354 cout << "Processor 0 reports the following inputs:" << endl;
355 cout << "inital.iterations = " << iterations << endl;
356 cout << "initial.temperature = " << temperature << endl;
357 cout << "initial.attraction = " << attraction << endl;
358 cout << "initial.damping_mult = " << damping_mult << endl;
359 cout << " ..." << endl;
360 cout << "liquid.iterations = " << liquid.iterations << endl;
361 cout << "liquid.temperature = " << liquid.temperature << endl;
362 cout << "liquid.attraction = " << liquid.attraction << endl;
363 cout << "liquid.damping_mult = " << liquid.damping_mult << endl;
364 cout << " ..." << endl;
365 cout << "simmer.iterations = " << simmer.iterations << endl;
366 cout << "simmer.temperature = " << simmer.temperature << endl;
367 cout << "simmer.attraction = " << simmer.attraction << endl;
368 cout << "simmer.damping_mult = " << simmer.damping_mult << endl;
369 }
370
371 }
372 */
373
374 // init_parms -- this subroutine initializes the edge_cut variables
375 // used in the original VxOrd starting with the edge_cut parameter.
376 // In our version, edge_cut = 0 means no cutting, 1 = maximum cut.
377 // We also set the random seed here.
378
init_parms(int rand_seed,float edge_cut,float real_parm)379 void graph::init_parms ( int rand_seed, float edge_cut, float real_parm ) {
380 IGRAPH_UNUSED(rand_seed);
381
382 // first we translate edge_cut the former tcl sliding scale
383 //CUT_END = cut_length_end = 39000.0 * (1.0 - edge_cut) + 1000.0;
384 CUT_END = cut_length_end = 40000.0 * (1.0 - edge_cut);
385
386 // cut_length_end cannot actually be 0
387 if ( cut_length_end <= 1.0 ) {
388 cut_length_end = 1.0;
389 }
390
391 float cut_length_start = 4.0 * cut_length_end;
392
393 // now we set the parameters used by ReCompute
394 cut_off_length = cut_length_start;
395 cut_rate = ( cut_length_start - cut_length_end ) / 400.0;
396
397 // finally set the number of iterations to leave .real coords fixed
398 int full_comp_iters;
399 full_comp_iters = liquid.iterations + expansion.iterations +
400 cooldown.iterations + crunch.iterations + 3;
401
402 // adjust real parm to iterations (do not enter simmer halfway)
403 if ( real_parm < 0 ) {
404 real_iterations = (int)real_parm;
405 } else if ( real_parm == 1) {
406 real_iterations = full_comp_iters + simmer.iterations + 100;
407 } else {
408 real_iterations = (int)(real_parm * full_comp_iters);
409 }
410
411 tot_iterations = 0;
412 if ( real_iterations > 0 ) {
413 real_fixed = true;
414 } else {
415 real_fixed = false;
416 }
417
418 // calculate total expected iterations (for progress bar display)
419 tot_expected_iterations = liquid.iterations +
420 expansion.iterations + cooldown.iterations +
421 crunch.iterations + simmer.iterations;
422
423 /*
424 // output edge_cutting parms (for debugging)
425 cout << "Processor " << myid << ": "
426 << "cut_length_end = CUT_END = " << cut_length_end
427 << ", cut_length_start = " << cut_length_start
428 << ", cut_rate = " << cut_rate << endl;
429 */
430
431 // set random seed
432 // srand ( rand_seed ); // Don't need this in igraph
433
434 }
435
init_parms(const igraph_layout_drl_options_t * options)436 void graph::init_parms(const igraph_layout_drl_options_t *options) {
437 double rand_seed = 0.0;
438 double real_in = -1.0;
439 init_parms(rand_seed, options->edge_cut, real_in);
440 }
441
442 // The following subroutine reads a .real file to obtain initial
443 // coordinates. If a node is missing coordinates the coordinates
444 // are computed
445
446 // void graph::read_real ( char *real_file )
447 // {
448 // cout << "Processor " << myid << " reading .real file ..." << endl;
449
450 // // read in .real file and mark as fixed
451 // ifstream real_in ( real_file );
452 // if ( !real_in )
453 // {
454 // cout << "Error: proc. " << myid << " could not open .real file." << endl;
455 // #ifdef MUSE_MPI
456 // MPI_Abort ( MPI_COMM_WORLD, 1 );
457 // #else
458 // exit (1);
459 // #endif
460 // }
461
462 // int real_id;
463 // float real_x, real_y;
464 // while ( !real_in.eof () )
465 // {
466 // real_id = -1;
467 // real_in >> real_id >> real_x >> real_y;
468 // if ( real_id >= 0 )
469 // {
470 // positions[id_catalog[real_id]].x = real_x;
471 // positions[id_catalog[real_id]].y = real_y;
472 // positions[id_catalog[real_id]].fixed = true;
473
474 // /*
475 // // output positions read (for debugging)
476 // cout << id_catalog[real_id] << " (" << positions[id_catalog[real_id]].x
477 // << ", " << positions[id_catalog[real_id]].y << ") "
478 // << positions[id_catalog[real_id]].fixed << endl;
479 // */
480
481 // // add node to density grid
482 // if ( real_iterations > 0 )
483 // density_server.Add ( positions[id_catalog[real_id]], fineDensity );
484 // }
485
486 // }
487
488 // real_in.close();
489 // }
490
read_real(const igraph_matrix_t * real_mat,const igraph_vector_bool_t * fixed)491 int graph::read_real ( const igraph_matrix_t *real_mat,
492 const igraph_vector_bool_t *fixed) {
493 long int n = igraph_matrix_nrow(real_mat);
494 for (long int i = 0; i < n; i++) {
495 positions[id_catalog[i]].x = MATRIX(*real_mat, i, 0);
496 positions[id_catalog[i]].y = MATRIX(*real_mat, i, 1);
497 positions[id_catalog[i]].fixed = fixed ? VECTOR(*fixed)[i] : false;
498
499 if ( real_iterations > 0 ) {
500 density_server.Add ( positions[id_catalog[i]], fineDensity );
501 }
502 }
503
504 return 0;
505 }
506
507 // The read_part_int subroutine reads the .int
508 // file produced by convert_sim and gathers the nodes and their
509 // neighbors in the range start_ind to end_ind.
510
511 // void graph::read_int ( char *file_name )
512 // {
513
514 // ifstream int_file;
515
516 // int_file.open ( file_name );
517 // if ( !int_file )
518 // {
519 // cout << "Error (worker process " << myid << "): could not open .int file." << endl;
520 // #ifdef MUSE_MPI
521 // MPI_Abort ( MPI_COMM_WORLD, 1 );
522 // #else
523 // exit (1);
524 // #endif
525 // }
526
527 // cout << "Processor " << myid << " reading .int file ..." << endl;
528
529 // int node_1, node_2;
530 // float weight;
531
532 // while ( !int_file.eof() )
533 // {
534 // weight = 0; // all weights should be >= 0
535 // int_file >> node_1 >> node_2 >> weight;
536 // if ( weight ) // otherwise we are at end of file
537 // // or it is a self-connected node
538 // {
539 // // normalization from original vxord
540 // weight /= highest_sim;
541 // weight = weight*fabs(weight);
542
543 // // initialize graph
544 // if ( ( node_1 % num_procs ) == myid )
545 // (neighbors[id_catalog[node_1]])[id_catalog[node_2]] = weight;
546 // if ( ( node_2 % num_procs ) == myid )
547 // (neighbors[id_catalog[node_2]])[id_catalog[node_1]] = weight;
548 // }
549 // }
550 // int_file.close();
551
552 // /*
553 // // the following code outputs the contents of the neighbors structure
554 // // (to be used for debugging)
555
556 // map<int, map<int,float> >::iterator i;
557 // map<int,float>::iterator j;
558
559 // for ( i = neighbors.begin(); i != neighbors.end(); i++ ) {
560 // cout << myid << ": " << i->first << " ";
561 // for (j = (i->second).begin(); j != (i->second).end(); j++ )
562 // cout << j->first << " (" << j->second << ") ";
563 // cout << endl;
564 // }
565 // */
566
567 // }
568
569 /*********************************************
570 * Function: ReCompute *
571 * Description: Compute the graph locations *
572 * Modified from original code by B. Wylie *
573 ********************************************/
574
ReCompute()575 int graph::ReCompute( ) {
576
577 // carryover from original VxOrd
578 int MIN = 1;
579
580 /*
581 // output parameters (for debugging)
582 cout << "ReCompute is using the following parameters: "<< endl;
583 cout << "STAGE: " << STAGE << ", iter: " << iterations << ", temp = " << temperature
584 << ", attract = " << attraction << ", damping_mult = " << damping_mult
585 << ", min_edges = " << min_edges << ", cut_off_length = " << cut_off_length
586 << ", fineDensity = " << fineDensity << endl;
587 */
588
589 /* igraph progress report */
590 float progress = (tot_iterations * 100.0 / tot_expected_iterations);
591
592 switch (STAGE) {
593 case 0:
594 if (iterations == 0) {
595 IGRAPH_PROGRESS("DrL layout (initialization stage)", progress, 0);
596 } else {
597 IGRAPH_PROGRESS("DrL layout (liquid stage)", progress, 0);
598 }
599 break;
600 case 1:
601 IGRAPH_PROGRESS("DrL layout (expansion stage)", progress, 0); break;
602 case 2:
603 IGRAPH_PROGRESS("DrL layout (cooldown and cluster phase)", progress, 0); break;
604 case 3:
605 IGRAPH_PROGRESS("DrL layout (crunch phase)", progress, 0); break;
606 case 5:
607 IGRAPH_PROGRESS("DrL layout (simmer phase)", progress, 0); break;
608 case 6:
609 IGRAPH_PROGRESS("DrL layout (final phase)", 100.0, 0); break;
610 default:
611 IGRAPH_PROGRESS("DrL layout (unknown phase)", 0.0, 0); break;
612 }
613
614 /* Compute Energies for individual nodes */
615 update_nodes ();
616
617 // check to see if we need to free fixed nodes
618 tot_iterations++;
619 if ( tot_iterations >= real_iterations ) {
620 real_fixed = false;
621 }
622
623
624 // ****************************************
625 // AUTOMATIC CONTROL SECTION
626 // ****************************************
627
628 // STAGE 0: LIQUID
629 if (STAGE == 0) {
630
631 if ( iterations == 0 ) {
632 start_time = time( NULL );
633 // if ( myid == 0 )
634 // cout << "Entering liquid stage ...";
635 }
636
637 if (iterations < liquid.iterations) {
638 temperature = liquid.temperature;
639 attraction = liquid.attraction;
640 damping_mult = liquid.damping_mult;
641 iterations++;
642 // if ( myid == 0 )
643 // cout << "." << flush;
644
645 } else {
646
647 stop_time = time( NULL );
648 liquid.time_elapsed = liquid.time_elapsed + (stop_time - start_time);
649 temperature = expansion.temperature;
650 attraction = expansion.attraction;
651 damping_mult = expansion.damping_mult;
652 iterations = 0;
653
654 // go to next stage
655 STAGE = 1;
656 start_time = time( NULL );
657
658 // if ( myid == 0 )
659 // cout << "Entering expansion stage ...";
660 }
661 }
662
663 // STAGE 1: EXPANSION
664 if (STAGE == 1) {
665
666 if (iterations < expansion.iterations) {
667
668 // Play with vars
669 if (attraction > 1) {
670 attraction -= .05;
671 }
672 if (min_edges > 12) {
673 min_edges -= .05;
674 }
675 cut_off_length -= cut_rate;
676 if (damping_mult > .1) {
677 damping_mult -= .005;
678 }
679 iterations++;
680 // if ( myid == 0 ) cout << "." << flush;
681
682 } else {
683
684 stop_time = time( NULL );
685 expansion.time_elapsed = expansion.time_elapsed + (stop_time - start_time);
686 min_edges = 12;
687 damping_mult = cooldown.damping_mult;
688
689 STAGE = 2;
690 attraction = cooldown.attraction;
691 temperature = cooldown.temperature;
692 iterations = 0;
693 start_time = time( NULL );
694
695 // if ( myid == 0 )
696 // cout << "Entering cool-down stage ...";
697 }
698 }
699
700 // STAGE 2: Cool down and cluster
701 else if (STAGE == 2) {
702
703 if (iterations < cooldown.iterations) {
704
705 // Reduce temperature
706 if (temperature > 50) {
707 temperature -= 10;
708 }
709
710 // Reduce cut length
711 if (cut_off_length > cut_length_end) {
712 cut_off_length -= cut_rate * 2;
713 }
714 if (min_edges > MIN) {
715 min_edges -= .2;
716 }
717 //min_edges = 99;
718 iterations++;
719 // if ( myid == 0 )
720 // cout << "." << flush;
721
722 } else {
723
724 stop_time = time( NULL );
725 cooldown.time_elapsed = cooldown.time_elapsed + (stop_time - start_time);
726 cut_off_length = cut_length_end;
727 temperature = crunch.temperature;
728 damping_mult = crunch.damping_mult;
729 min_edges = MIN;
730 //min_edges = 99; // In other words: no more cutting
731
732 STAGE = 3;
733 iterations = 0;
734 attraction = crunch.attraction;
735 start_time = time( NULL );
736
737 // if ( myid == 0 )
738 // cout << "Entering crunch stage ...";
739 }
740 }
741
742 // STAGE 3: Crunch
743 else if (STAGE == 3) {
744
745 if (iterations < crunch.iterations) {
746 iterations++;
747 // if ( myid == 0 ) cout << "." << flush;
748 } else {
749
750 stop_time = time( NULL );
751 crunch.time_elapsed = crunch.time_elapsed + (stop_time - start_time);
752 iterations = 0;
753 temperature = simmer.temperature;
754 attraction = simmer.attraction;
755 damping_mult = simmer.damping_mult;
756 min_edges = 99;
757 fineDensity = true;
758
759 STAGE = 5;
760 start_time = time( NULL );
761
762 // if ( myid == 0 )
763 // cout << "Entering simmer stage ...";
764 }
765 }
766
767 // STAGE 5: Simmer
768 else if ( STAGE == 5 ) {
769
770 if (iterations < simmer.iterations) {
771 if (temperature > 50) {
772 temperature -= 2;
773 }
774 iterations++;
775 // if ( myid == 0 ) cout << "." << flush;
776 } else {
777 stop_time = time( NULL );
778 simmer.time_elapsed = simmer.time_elapsed + (stop_time - start_time);
779
780 STAGE = 6;
781
782 // if ( myid == 0 )
783 // cout << "Layout calculation completed in " <<
784 // ( liquid.time_elapsed + expansion.time_elapsed +
785 // cooldown.time_elapsed + crunch.time_elapsed +
786 // simmer.time_elapsed )
787 // << " seconds (not including I/O)."
788 // << endl;
789 }
790 }
791
792 // STAGE 6: All Done!
793 else if ( STAGE == 6) {
794
795 /*
796 // output parameters (for debugging)
797 cout << "ReCompute is using the following parameters: "<< endl;
798 cout << "STAGE: " << STAGE << ", iter: " << iterations << ", temp = " << temperature
799 << ", attract = " << attraction << ", damping_mult = " << damping_mult
800 << ", min_edges = " << min_edges << ", cut_off_length = " << cut_off_length
801 << ", fineDensity = " << fineDensity << endl;
802 */
803
804 return 0;
805 }
806
807 // ****************************************
808 // END AUTOMATIC CONTROL SECTION
809 // ****************************************
810
811 // Still need more recomputation
812 return 1;
813
814 }
815
816 // update_nodes -- this function will complete the primary node update
817 // loop in layout's recompute routine. It follows exactly the same
818 // sequence to ensure similarity of parallel layout to the standard layout
819
update_nodes()820 void graph::update_nodes ( ) {
821
822 vector<int> node_indices; // node list of nodes currently being updated
823 float old_positions[2 * MAX_PROCS]; // positions before update
824 float new_positions[2 * MAX_PROCS]; // positions after update
825
826 bool all_fixed; // check if all nodes are fixed
827
828 // initial node list consists of 0,1,...,num_procs
829 for ( int i = 0; i < num_procs; i++ ) {
830 node_indices.push_back( i );
831 }
832
833 // next we calculate the number of nodes there would be if the
834 // num_nodes by num_procs schedule grid were perfectly square
835 int square_num_nodes = (int)(num_procs + num_procs * floor ((float)(num_nodes - 1) / (float)num_procs ));
836
837 for ( int i = myid; i < square_num_nodes; i += num_procs ) {
838
839 // get old positions
840 get_positions ( node_indices, old_positions );
841
842 // default new position is old position
843 get_positions ( node_indices, new_positions );
844
845 if ( i < num_nodes ) {
846
847 // advance random sequence according to myid
848 for ( int j = 0; j < 2 * myid; j++ ) {
849 RNG_UNIF01();
850 }
851 // rand();
852
853 // calculate node energy possibilities
854 if ( !(positions[i].fixed && real_fixed) ) {
855 update_node_pos ( i, old_positions, new_positions );
856 }
857
858 // advance random sequence for next iteration
859 for ( unsigned int j = 2 * myid; j < 2 * (node_indices.size() - 1); j++ ) {
860 RNG_UNIF01();
861 }
862 // rand();
863
864 } else {
865 // advance random sequence according to use by
866 // the other processors
867 for ( unsigned int j = 0; j < 2 * (node_indices.size()); j++ ) {
868 RNG_UNIF01();
869 }
870 //rand();
871 }
872
873 // check if anything was actually updated (e.g. everything was fixed)
874 all_fixed = true;
875 for ( unsigned int j = 0; j < node_indices.size (); j++ )
876 if ( !(positions [ node_indices[j] ].fixed && real_fixed) ) {
877 all_fixed = false;
878 }
879
880 // update positions across processors (if not all fixed)
881 if ( !all_fixed ) {
882 #ifdef MUSE_MPI
883 MPI_Allgather ( &new_positions[2 * myid], 2, MPI_FLOAT,
884 new_positions, 2, MPI_FLOAT, MPI_COMM_WORLD );
885 #endif
886
887 // update positions (old to new)
888 update_density ( node_indices, old_positions, new_positions );
889 }
890
891 /*
892 if ( myid == 0 )
893 {
894 // output node list (for debugging)
895 for ( unsigned int j = 0; j < node_indices.size(); j++ )
896 cout << node_indices[j] << " ";
897 cout << endl;
898 }
899 */
900
901 // compute node list for next update
902 for ( unsigned int j = 0; j < node_indices.size(); j++ ) {
903 node_indices [j] += num_procs;
904 }
905
906 while ( !node_indices.empty() && node_indices.back() >= num_nodes ) {
907 node_indices.pop_back ( );
908 }
909
910 }
911
912 // update first_add and fine_first_add
913 first_add = false;
914 if ( fineDensity ) {
915 fine_first_add = false;
916 }
917
918 }
919
920 // The get_positions function takes the node_indices list
921 // and returns the corresponding positions in an array.
922
get_positions(vector<int> & node_indices,float return_positions[2* MAX_PROCS])923 void graph::get_positions ( vector<int> &node_indices,
924 float return_positions[2 * MAX_PROCS] ) {
925
926 // fill positions
927 for (unsigned int i = 0; i < node_indices.size(); i++) {
928 return_positions[2 * i] = positions[ node_indices[i] ].x;
929 return_positions[2 * i + 1] = positions[ node_indices[i] ].y;
930 }
931
932 }
933
934 // update_node_pos -- this subroutine does the actual work of computing
935 // the new position of a given node. num_act_proc gives the number
936 // of active processes at this level for use by the random number
937 // generators.
938
update_node_pos(int node_ind,float old_positions[2* MAX_PROCS],float new_positions[2* MAX_PROCS])939 void graph::update_node_pos ( int node_ind,
940 float old_positions[2 * MAX_PROCS],
941 float new_positions[2 * MAX_PROCS] ) {
942
943 float energies[2]; // node energies for possible positions
944 float updated_pos[2][2]; // possible positions
945 float pos_x, pos_y;
946
947 // old VxOrd parameter
948 float jump_length = .010 * temperature;
949
950 // subtract old node
951 density_server.Subtract ( positions[node_ind], first_add, fine_first_add, fineDensity );
952
953 // compute node energy for old solution
954 energies[0] = Compute_Node_Energy ( node_ind );
955
956 // move node to centroid position
957 Solve_Analytic ( node_ind, pos_x, pos_y );
958 positions[node_ind].x = updated_pos[0][0] = pos_x;
959 positions[node_ind].y = updated_pos[0][1] = pos_y;
960
961 /*
962 // ouput random numbers (for debugging)
963 int rand_0, rand_1;
964 rand_0 = rand();
965 rand_1 = rand();
966 cout << myid << ": " << rand_0 << ", " << rand_1 << endl;
967 */
968
969 // Do random method (RAND_MAX is C++ maximum random number)
970 updated_pos[1][0] = updated_pos[0][0] + (.5 - RNG_UNIF01()) * jump_length;
971 updated_pos[1][1] = updated_pos[0][1] + (.5 - RNG_UNIF01()) * jump_length;
972
973 // compute node energy for random position
974 positions[node_ind].x = updated_pos[1][0];
975 positions[node_ind].y = updated_pos[1][1];
976 energies[1] = Compute_Node_Energy ( node_ind );
977
978 /*
979 // output update possiblities (debugging):
980 cout << node_ind << ": (" << updated_pos[0][0] << "," << updated_pos[0][1]
981 << "), " << energies[0] << "; (" << updated_pos[1][0] << ","
982 << updated_pos[1][1] << "), " << energies[1] << endl;
983 */
984
985 // add back old position
986 positions[node_ind].x = old_positions[2 * myid];
987 positions[node_ind].y = old_positions[2 * myid + 1];
988 if ( !fineDensity && !first_add ) {
989 density_server.Add ( positions[node_ind], fineDensity );
990 } else if ( !fine_first_add ) {
991 density_server.Add ( positions[node_ind], fineDensity );
992 }
993
994 // choose updated node position with lowest energy
995 if ( energies[0] < energies[1] ) {
996 new_positions[2 * myid] = updated_pos[0][0];
997 new_positions[2 * myid + 1] = updated_pos[0][1];
998 positions[node_ind].energy = energies[0];
999 } else {
1000 new_positions[2 * myid] = updated_pos[1][0];
1001 new_positions[2 * myid + 1] = updated_pos[1][1];
1002 positions[node_ind].energy = energies[1];
1003 }
1004
1005 }
1006
1007 // update_density takes a sequence of node_indices and their positions and
1008 // updates the positions by subtracting the old positions and adding the
1009 // new positions to the density grid.
1010
update_density(vector<int> & node_indices,float old_positions[2* MAX_PROCS],float new_positions[2* MAX_PROCS])1011 void graph::update_density ( vector<int> &node_indices,
1012 float old_positions[2 * MAX_PROCS],
1013 float new_positions[2 * MAX_PROCS] ) {
1014
1015 // go through each node and subtract old position from
1016 // density grid before adding new position
1017 for ( unsigned int i = 0; i < node_indices.size(); i++ ) {
1018 positions[node_indices[i]].x = old_positions[2 * i];
1019 positions[node_indices[i]].y = old_positions[2 * i + 1];
1020 density_server.Subtract ( positions[node_indices[i]],
1021 first_add, fine_first_add, fineDensity );
1022
1023 positions[node_indices[i]].x = new_positions[2 * i];
1024 positions[node_indices[i]].y = new_positions[2 * i + 1];
1025 density_server.Add ( positions[node_indices[i]], fineDensity );
1026 }
1027
1028 }
1029
1030 /********************************************
1031 * Function: Compute_Node_Energy *
1032 * Description: Compute the node energy *
1033 * This code has been modified from the *
1034 * original code by B. Wylie. *
1035 *********************************************/
1036
Compute_Node_Energy(int node_ind)1037 float graph::Compute_Node_Energy( int node_ind ) {
1038
1039 /* Want to expand 4th power range of attraction */
1040 float attraction_factor = attraction * attraction *
1041 attraction * attraction * 2e-2;
1042
1043 map <int, float>::iterator EI;
1044 float x_dis, y_dis;
1045 float energy_distance, weight;
1046 float node_energy = 0;
1047
1048 // Add up all connection energies
1049 for (EI = neighbors[node_ind].begin(); EI != neighbors[node_ind].end(); ++EI) {
1050
1051 // Get edge weight
1052 weight = EI->second;
1053
1054 // Compute x,y distance
1055 x_dis = positions[ node_ind ].x - positions[ EI->first ].x;
1056 y_dis = positions[ node_ind ].y - positions[ EI->first ].y;
1057
1058 // Energy Distance
1059 energy_distance = x_dis * x_dis + y_dis * y_dis;
1060 if (STAGE < 2) {
1061 energy_distance *= energy_distance;
1062 }
1063
1064 // In the liquid phase we want to discourage long link distances
1065 if (STAGE == 0) {
1066 energy_distance *= energy_distance;
1067 }
1068
1069 node_energy += weight * attraction_factor * energy_distance;
1070 }
1071
1072 // output effect of density (debugging)
1073 //cout << "[before: " << node_energy;
1074
1075 // add density
1076 node_energy += density_server.GetDensity ( positions[ node_ind ].x, positions[ node_ind ].y,
1077 fineDensity );
1078
1079 // after calling density server (debugging)
1080 //cout << ", after: " << node_energy << "]" << endl;
1081
1082 // return computated energy
1083 return node_energy;
1084 }
1085
1086
1087 /*********************************************
1088 * Function: Solve_Analytic *
1089 * Description: Compute the node position *
1090 * This is a modified version of the function *
1091 * originally written by B. Wylie *
1092 *********************************************/
1093
Solve_Analytic(int node_ind,float & pos_x,float & pos_y)1094 void graph::Solve_Analytic( int node_ind, float &pos_x, float &pos_y ) {
1095
1096 map <int, float>::iterator EI;
1097 float total_weight = 0;
1098 float x_dis, y_dis, x_cen = 0, y_cen = 0;
1099 float x = 0, y = 0, dis;
1100 float damping, weight;
1101
1102 // Sum up all connections
1103 for (EI = neighbors[node_ind].begin(); EI != neighbors[node_ind].end(); ++EI) {
1104 weight = EI->second;
1105 total_weight += weight;
1106 x += weight * positions[ EI->first ].x;
1107 y += weight * positions[ EI->first ].y;
1108 }
1109
1110 // Now set node position
1111 if (total_weight > 0) {
1112
1113 // Compute centriod
1114 x_cen = x / total_weight;
1115 y_cen = y / total_weight;
1116 damping = 1.0 - damping_mult;
1117 pos_x = damping * positions[ node_ind ].x + (1.0 - damping) * x_cen;
1118 pos_y = damping * positions[ node_ind ].y + (1.0 - damping) * y_cen;
1119 } else {
1120 pos_x = positions[ node_ind ].x;
1121 pos_y = positions[ node_ind ].y;
1122 }
1123
1124 // No cut edge flag (?)
1125 if (min_edges == 99) {
1126 return;
1127 }
1128
1129 // Don't cut at end of scale
1130 if ( CUT_END >= 39500 ) {
1131 return;
1132 }
1133
1134 float num_connections = sqrt((double)neighbors[node_ind].size());
1135 float maxLength = 0;
1136
1137 map<int, float>::iterator maxIndex;
1138
1139 // Go through nodes edges... cutting if necessary
1140 for (EI = maxIndex = neighbors[node_ind].begin();
1141 EI != neighbors[node_ind].end(); ++EI) {
1142
1143 // Check for at least min edges
1144 if (neighbors[node_ind].size() < min_edges) {
1145 continue;
1146 }
1147
1148 x_dis = x_cen - positions[ EI->first ].x;
1149 y_dis = y_cen - positions[ EI->first ].y;
1150 dis = x_dis * x_dis + y_dis * y_dis;
1151 dis *= num_connections;
1152
1153 // Store maximum edge
1154 if (dis > maxLength) {
1155 maxLength = dis;
1156 maxIndex = EI;
1157 }
1158 }
1159
1160 // If max length greater than cut_length then cut
1161 if (maxLength > cut_off_length) {
1162 neighbors[ node_ind ].erase( maxIndex );
1163 }
1164
1165 }
1166
1167
1168 // write_coord writes out the coordinate file of the final solutions
1169
1170 // void graph::write_coord( const char *file_name )
1171 // {
1172
1173 // ofstream coordOUT( file_name );
1174 // if ( !coordOUT )
1175 // {
1176 // cout << "Could not open " << file_name << ". Program terminated." << endl;
1177 // #ifdef MUSE_MPI
1178 // MPI_Abort ( MPI_COMM_WORLD, 1 );
1179 // #else
1180 // exit (1);
1181 // #endif
1182 // }
1183
1184 // cout << "Writing out solution to " << file_name << " ..." << endl;
1185
1186 // for (unsigned int i = 0; i < positions.size(); i++) {
1187 // coordOUT << positions[i].id << "\t" << positions[i].x << "\t" << positions[i].y <<endl;
1188 // }
1189 // coordOUT.close();
1190
1191 // }
1192
1193 // write_sim -- outputs .edges file, takes as input .coord filename,
1194 // with .coord extension
1195
1196 /*
1197 void graph::write_sim ( const char *file_name )
1198 {
1199
1200 string prefix_name ( file_name, strlen(file_name)-7 );
1201 prefix_name = prefix_name + ".iedges";
1202
1203 // first we overwrite, then we append
1204 ofstream simOUT;
1205 if ( myid == 0 )
1206 simOUT.open ( prefix_name.c_str() );
1207 else
1208 simOUT.open ( prefix_name.c_str(), ios::app );
1209
1210 if ( !simOUT )
1211 {
1212 cout << "Could not open " << prefix_name << ". Program terminated." << endl;
1213 #ifdef MUSE_MPI
1214 MPI_Abort ( MPI_COMM_WORLD, 1 );
1215 #else
1216 exit (1);
1217 #endif
1218 }
1219
1220
1221 cout << "Proc. " << myid << " writing to " << prefix_name << " ..." << endl;
1222
1223
1224 // the following code outputs the contents of the neighbors structure
1225
1226 map<int, map<int,float> >::iterator i;
1227 map<int,float>::iterator j;
1228
1229 for ( i = neighbors.begin(); i != neighbors.end(); i++ )
1230 for (j = (i->second).begin(); j != (i->second).end(); j++ )
1231 simOUT << positions[i->first].id << "\t"
1232 << positions[j->first].id << "\t"
1233 << j->second << endl;
1234
1235 simOUT.close();
1236
1237 }
1238 */
1239
1240 // get_tot_energy adds up the energy for each node to give an estimate of the
1241 // quality of the minimization.
1242
get_tot_energy()1243 float graph::get_tot_energy ( ) {
1244
1245 float my_tot_energy, tot_energy;
1246 my_tot_energy = 0;
1247 for ( int i = myid; i < num_nodes; i += num_procs ) {
1248 my_tot_energy += positions[i].energy;
1249 }
1250
1251 //vector<Node>::iterator i;
1252 //for ( i = positions.begin(); i != positions.end(); i++ )
1253 // tot_energy += i->energy;
1254
1255 #ifdef MUSE_MPI
1256 MPI_Reduce ( &my_tot_energy, &tot_energy, 1, MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD );
1257 #else
1258 tot_energy = my_tot_energy;
1259 #endif
1260
1261 return tot_energy;
1262
1263 }
1264
1265
1266 // The following subroutine draws the graph with possible intermediate
1267 // output (int_out is set to 0 if not proc. 0). int_out is the parameter
1268 // passed by the user, and coord_file is the .coord file.
1269
1270 // void graph::draw_graph ( int int_out, char *coord_file )
1271 // {
1272
1273 // // layout graph (with possible intermediate output)
1274 // int count_iter = 0, count_file = 1;
1275 // char int_coord_file [MAX_FILE_NAME + MAX_INT_LENGTH];
1276 // while ( ReCompute( ) )
1277 // if ( (int_out > 0) && (count_iter == int_out) )
1278 // {
1279 // // output intermediate solution
1280 // sprintf ( int_coord_file, "%s.%d", coord_file, count_file );
1281 // write_coord ( int_coord_file );
1282
1283 // count_iter = 0;
1284 // count_file++;
1285 // }
1286 // else
1287 // count_iter++;
1288
1289 // }
1290
draw_graph(igraph_matrix_t * res)1291 int graph::draw_graph(igraph_matrix_t *res) {
1292 int count_iter = 0;
1293 while (ReCompute()) {
1294 IGRAPH_ALLOW_INTERRUPTION();
1295 count_iter++;
1296 }
1297 long int n = positions.size();
1298 IGRAPH_CHECK(igraph_matrix_resize(res, n, 2));
1299 for (long int i = 0; i < n; i++) {
1300 MATRIX(*res, i, 0) = positions[i].x;
1301 MATRIX(*res, i, 1) = positions[i].y;
1302 }
1303 return 0;
1304 }
1305
1306 } // namespace drl
1307