// This file contains the member definitions of the master class
#pragma warning(disable: 4786)

#include <iostream>
#include <fstream>
#include <map>
#include <vector>
#include <math.h>
#include <time.h>
#include <string>

using namespace std;

#include "ebGraph2d.h"

// constructor -- initializes the schedule variables (as in 
// graph constructor)

eb_graph2d::eb_graph2d ( char *wel_file )
{
		  
	  // MPI parameters
	  myid = 0;
	  num_procs = 1;

	  // initial annealing parameters
	  STAGE = 0;
	  iterations = 0;
	  temperature = 2000;
	  attraction = 10;
	  damping_mult = 1.0;
	  min_edges = 20;
	  first_add = fine_first_add = true;
	  fineDensity = false;

	  // WES
	  current_id = 0;
	  
	  // Brian's original Vx schedule
	  liquid.iterations = 200;
	  liquid.temperature = 2000;
	  liquid.attraction = 2;
	  liquid.damping_mult = 1.0;
	  liquid.time_elapsed = 0;

	  expansion.iterations = 200;
	  expansion.temperature = 2000;
	  expansion.attraction = 10;
	  expansion.damping_mult = 1.0;
	  expansion.time_elapsed = 0;

	  cooldown.iterations = 200;
	  cooldown.temperature = 2000;
	  cooldown.attraction = 1;
	  cooldown.damping_mult = .1f;
	  cooldown.time_elapsed = 0;

	  crunch.iterations = 50;
	  crunch.temperature = 250;
	  crunch.attraction = 1;
	  crunch. damping_mult = .25;
	  crunch.time_elapsed = 0;

	  simmer.iterations = 100;
	  simmer.temperature = 250;
	  simmer.attraction = .5;
	  simmer.damping_mult = 0.0;
	  simmer.time_elapsed = 0;

	  // scan .wel file for node info
	  scan_wel ( wel_file );
	  
	  // populate node positions and ids
	  positions.reserve ( num_nodes );
	  map < int, int >::iterator cat_iter;
	  for ( cat_iter = id_catalog.begin();
		    cat_iter != id_catalog.end();
			cat_iter++ )
		positions.push_back ( Node2d( cat_iter->first ) );
	  
	  /*
	  // output positions .ids for debugging
	  for ( int id = 0; id < num_nodes; id++ )
		cout << positions[id].id << endl;
	  */
	  
	  // read .wel file for graph info
	  read_wel ( wel_file );
	  
	  // initialize density server
	  density_server.Init();
}

// The following subroutine scans the .wel file for the following
// information: number nodes, node ids, and highest similarity.  The
// corresponding graph globals are populated: num_nodes, id_catalog,
// and highest_sim.

void eb_graph2d::scan_wel ( char *filename )
{

  cout << "Proc. " << myid << " scanning .wel file ..." << endl;
  
  // Open (sim) File
  ifstream fp ( filename );
  if ( !fp )
  {
	cout << "Error: could not open " << filename << ".  Program terminated." << endl;
	  exit (1);
  }	
  
  // Read file, parse, and add into data structure
  int id1, id2;
  string label1, label2, dummy;
  float edge_weight;
  highest_sim = -1.0;

  // scan number of vertices and edges 
  fp >> num_vertices; fp >> num_edges;
  while ( !fp.eof () )
	{
	  fp >> label1 >> label2 >> edge_weight;
	  
	  // ignore negative weights!
	  if ( edge_weight <= 0 )
	  {
		  edge_weight *= -1;
	   }

	   if ( highest_sim < edge_weight )
	      highest_sim = edge_weight;
	
	   if (vertex_ids.find(label1) != vertex_ids.end())
	   {
		   id1 = vertex_ids[label1];
	   }
	   else
	   {
		   id1 = current_id;
		   vertex_ids[label1] = id1;
		   vertex_labels[id1] = label1;
		   current_id++;
	   }
	   if (vertex_ids.find(label2) != vertex_ids.end())
	   {
		   id2 = vertex_ids[label2];
	   }
	   else
	   {
		   id2 = current_id;
		   vertex_ids[label2] = id2;
		   vertex_labels[id2] = label2;
		   current_id++;
	   }
	   id_catalog[id1] = 1;
	   id_catalog[id2] = 1;
	}

  fp.close();

  if ( id_catalog.size() == 0 )
  {
    cout << "Error: Proc. " << myid << ": " << filename << " is empty.  Program terminated." << endl;
	  exit (1);
  }
  
  // label nodes with sequential integers starting at 0
  map< int, int>::iterator cat_iter;
  int id_label;
  for ( cat_iter = id_catalog.begin(), id_label = 0;
	    cat_iter != id_catalog.end(); cat_iter++, id_label++ )
    cat_iter->second = id_label;

  num_nodes = id_catalog.size();  
}


// init_parms -- this subroutine initializes the edge_cut variables
// used in the original VxOrd starting with the edge_cut parameter.
// In our version, edge_cut = 0 means no cutting, 1 = maximum cut.
// We also set the random seed here.

void eb_graph2d::init_parms ( int rand_seed, float edge_cut, float real_parm )
{

	// first we translate edge_cut the former tcl sliding scale
	//CUT_END = cut_length_end = 39000.0 * (1.0 - edge_cut) + 1000.0;
	CUT_END = cut_length_end = 40000.f * (1.f - edge_cut);
	
	// cut_length_end cannot actually be 0
	if ( cut_length_end <= 1.0 )
	  cut_length_end = 1.0;
	  
	float cut_length_start = 4.f * cut_length_end;
	
	// now we set the parameters used by ReCompute
	cut_off_length = cut_length_start;
	cut_rate = ( cut_length_start - cut_length_end ) / 400.f;
	
	// finally set the number of iterations to leave .real coords fixed
	int full_comp_iters;
	full_comp_iters = liquid.iterations + expansion.iterations +
					  cooldown.iterations + crunch.iterations + 3;
    
	// adjust real parm to iterations (do not enter simmer halfway)
	if ( real_parm < 0 )
		real_iterations = (int)real_parm;
	else if ( real_parm == 1)
		real_iterations = full_comp_iters + simmer.iterations + 100;
	else
		real_iterations = (int)(real_parm*full_comp_iters);
		
	tot_iterations = 0;
	if ( real_iterations > 0 )
		real_fixed = true;
	else
		real_fixed = false;
		
	// set random seed
	srand ( rand_seed );
	
}


// The read_part_int subroutine reads the .wel
// file produced by convert_sim and gathers the nodes and their
// neighbors in the range start_ind to end_ind.

void eb_graph2d::read_wel ( char *file_name )
{

	ifstream wel_file;
	
	wel_file.open ( file_name );
	if ( !wel_file )
	{
		cout << "Error (worker process " << myid << "): could not open .wel file." << endl;
		  exit (1);
	}
	
	cout << "Processor " << myid << " reading .wel file ..." << endl;
	
	string node_1, node_2, dummy;
	float weight;
	
	// bypass first two number of wel file
	wel_file >> dummy; wel_file >> dummy;
    while ( !wel_file.eof() )
	{
		weight = 0;		// all weights should be >= 0
		wel_file >> node_1 >> node_2 >> weight;
		weight = (float)fabs(weight);
		if ( weight )		// otherwise we are at end of file
								// or it is a self-connected node
		{
			    // normalization from original vxord
			    weight /= highest_sim;
				weight = weight*(float)fabs(weight);
				
				// initialize graph
				if ( ( vertex_ids[node_1] % num_procs ) == myid )
					(neighbors[id_catalog[vertex_ids[node_1]]])[id_catalog[vertex_ids[node_2]]] = weight;
				if ( ( vertex_ids[node_2] % num_procs ) == myid )
					(neighbors[id_catalog[vertex_ids[node_2]]])[id_catalog[vertex_ids[node_1]]] = weight;
		}
	}
	wel_file.close();
}

/*********************************************
 * Function: ReCompute				         *
 * Description: Compute the graph locations	 *
 * Modified from original code by B. Wylie   *
 ********************************************/

int eb_graph2d::ReCompute( ) 
{

  // carryover from original VxOrd
  int MIN = 1;
  
  // energy information is across processors
  float tot_energy;
  
  /* Compute Energies for individual nodes */
  update_nodes ();
  
  // check to see if we need to free fixed nodes
  tot_iterations++;
  if ( tot_iterations >= real_iterations )
	 real_fixed = false;


  	// ****************************************
	// AUTOMATIC CONTROL SECTION
	// ****************************************

	// STAGE 0: LIQUID
	if (STAGE == 0) {

		if ( iterations == 0 )
		{
			start_time = time( NULL );
			if ( myid == 0 )
				cout << "Entering liquid stage ...";
		}

		if (iterations < liquid.iterations) {
			temperature = liquid.temperature;
			attraction = liquid.attraction;
			damping_mult = liquid.damping_mult;
			iterations++;
			if ( myid == 0 )
			{
				int prev_percent = 
					(int)((iterations - 1) 
					/ (float)liquid.iterations * 100);
				int percent = 
					(int)(iterations 
					/ (float)liquid.iterations * 100);
					
				if (iterations != 1)
				{
					if (prev_percent < 10)
					{
						cout << "\b\b";
					}
					else
					{
						cout << "\b\b\b";
					}
				}
				cout << percent << "%" << flush;
			}
		
		} else {

			stop_time = time( NULL );
			liquid.time_elapsed = liquid.time_elapsed + (stop_time - start_time);
			temperature = expansion.temperature;
			attraction = expansion.attraction;
			damping_mult = expansion.damping_mult;
			iterations = 0;

		    tot_energy = get_tot_energy ();
			if ( myid == 0 )
				cout << endl << "Liquid stage completed in " << liquid.time_elapsed
					 << " seconds, total energy = " << tot_energy << "." << endl;
			
			// go to next stage
			STAGE = 1;
			start_time = time( NULL );
			
			if ( myid == 0 )
				cout << "Entering expansion stage ...";
		}
	}

	// STAGE 1: EXPANSION
	if (STAGE == 1) {

		if (iterations < expansion.iterations) {
				
			// Play with vars
			if (attraction > 1) attraction -= .05f;
			if (min_edges > 12) min_edges -= .05f;
			cut_off_length -= cut_rate;
			if (damping_mult > .1) damping_mult -= .005f;
			iterations++;
			if ( myid == 0 ) 
			{
				int prev_percent = 
					(int)((iterations - 1) 
					/ (float)expansion.iterations * 100);
				int percent = 
					(int)(iterations 
					/ (float)expansion.iterations * 100);
					
				if (iterations != 1)
				{
					if (prev_percent < 10)
					{
						cout << "\b\b";
					}
					else
					{
						cout << "\b\b\b";
					}
				}
				cout << percent << "%" << flush;
			}
	
		} else {

			stop_time = time( NULL );
			expansion.time_elapsed = expansion.time_elapsed + (stop_time - start_time);
		  	min_edges = 12;
			damping_mult = cooldown.damping_mult;
			
			tot_energy = get_tot_energy ();
			if ( myid == 0 )
				cout << endl << "Finished expansion stage in " << expansion.time_elapsed
					 << " seconds, total energy = " << tot_energy << "." << endl;
			
			STAGE = 2;
			attraction = cooldown.attraction;
			temperature = cooldown.temperature;
			iterations = 0;
			start_time = time( NULL );
			
			if ( myid == 0 )
				cout << "Entering cool-down stage ...";
		}
	}

	// STAGE 2: Cool down and cluster
	else if(STAGE==2) {

		if (iterations < cooldown.iterations) {

			// Reduce temperature
			if (temperature > 50) temperature -= 10;

			// Reduce cut length
			if (cut_off_length > cut_length_end) cut_off_length -= cut_rate*2;
			if (min_edges > MIN) min_edges -= .2f;
			//min_edges = 99;
			iterations++;
			if ( myid == 0 )
			{
				int prev_percent = 
					(int)((iterations - 1) 
					/ (float)cooldown.iterations * 100);
				int percent = 
					(int)(iterations 
					/ (float)cooldown.iterations * 100);
					
				if (iterations != 1)
				{
					if (prev_percent < 10)
					{
						cout << "\b\b";
					}
					else
					{
						cout << "\b\b\b";
					}
				}
				cout << percent << "%" << flush;
			}
			
		} else {

			stop_time = time( NULL );
			cooldown.time_elapsed = cooldown.time_elapsed + (stop_time - start_time);
			cut_off_length = cut_length_end;
			temperature = crunch.temperature;
			damping_mult = crunch.damping_mult;
			min_edges = (float)MIN;
			//min_edges = 99; // In other words: no more cutting
			
			tot_energy = get_tot_energy ();
			if ( myid == 0 )
				cout << endl << "Completed cool-down stage in " << cooldown.time_elapsed
					 << " seconds, total energy = " << tot_energy << "." << endl;
			
			STAGE = 3;
			iterations = 0;
			attraction = crunch.attraction;
			start_time = time( NULL );
			
			if ( myid == 0 )
				cout << "Entering crunch stage ...";
		}  
	}

	// STAGE 3: Crunch
	else if(STAGE==3) {
	
		if (iterations < crunch.iterations)
		{
			iterations++;
			if ( myid == 0 ) 
			{
				int prev_percent = 
					(int)((iterations - 1) 
					/ (float)crunch.iterations * 100);
				int percent = 
					(int)(iterations 
					/ (float)crunch.iterations * 100);
					
				if (iterations != 1)
				{
					if (prev_percent < 10)
					{
						cout << "\b\b";
					}
					else
					{
						cout << "\b\b\b";
					}
				}
				cout << percent << "%" << flush;
			}
		}
		else {
		
			stop_time = time( NULL );
			crunch.time_elapsed = crunch.time_elapsed + (stop_time - start_time);
			iterations = 0;
			temperature = simmer.temperature;
			attraction = simmer.attraction;
			damping_mult = simmer.damping_mult;
			min_edges = 99;
			fineDensity = true;
			
			tot_energy = get_tot_energy ();
			if ( myid == 0 )
				cout << endl << "Finished crunch stage in " << crunch.time_elapsed
					 << " seconds, total energy = " << tot_energy << "." << endl;
			
			STAGE = 5;
			start_time = time( NULL );
				
			if ( myid == 0 )
				cout << "Entering simmer stage ...";
		}  
	}

	// STAGE 5: Simmer
	else if( STAGE==5 ) {

		if (iterations < simmer.iterations) {
			if (temperature > 50) temperature -= 2;
			iterations++;
			if ( myid == 0 )
			{
				int prev_percent = 
					(int)((iterations - 1) 
					/ (float)simmer.iterations * 100);
				int percent = 
					(int)(iterations 
					/ (float)simmer.iterations * 100);
					
				if (iterations != 1)
				{
					if (prev_percent < 10)
					{
						cout << "\b\b";
					}
					else
					{
						cout << "\b\b\b";
					}
				}
				cout << percent << "%" << flush;
			}
		
		} else {
			stop_time = time( NULL );
			simmer.time_elapsed = simmer.time_elapsed + (stop_time - start_time);
			
			tot_energy = get_tot_energy();
			if ( myid == 0 )
				cout << endl << "Finished simmer stage in " << simmer.time_elapsed
					 << " seconds, total energy = "<< tot_energy << "." << endl;

			STAGE = 6;
			
			if ( myid == 0 )
				cout << "Layout calculation completed in " <<
				  ( liquid.time_elapsed + expansion.time_elapsed +
				    cooldown.time_elapsed + crunch.time_elapsed +
				    simmer.time_elapsed )
				     << " seconds (not including I/O)." 
				     << endl;
		}
	}

	// STAGE 6: All Done!
	else if ( STAGE == 6)
	{
	  
	  /*
	  // output parameters (for debugging)
	  cout << "ReCompute is using the following parameters: "<< endl;
      cout << "STAGE: " << STAGE << ", iter: " << iterations << ", temp = " << temperature
           << ", attract = " << attraction << ", damping_mult = " << damping_mult
	       << ", min_edges = " << min_edges << ", cut_off_length = " << cut_off_length
	       << ", fineDensity = " << fineDensity << endl; 
	  */
	  
	  return 0;
	}
	
	// ****************************************
	// END AUTOMATIC CONTROL SECTION
	// ****************************************

	// Still need more recomputation
	return 1;

}

// update_nodes -- this function will complete the primary node update
// loop in layout's recompute routine.  It follows exactly the same
// sequence to ensure similarity of parallel layout to the standard layout

void eb_graph2d::update_nodes ( )
{
	int i, j;
   unsigned int k;
	vector<int> node_indices;			// node list of nodes currently being updated
	float old_positions[3];	// positions before update
	float new_positions[3];	// positions after update
    
	bool all_fixed;						// check if all nodes are fixed
	
	// initial node list consists of 0,1,...,num_procs
	for ( i = 0; i < num_procs; i++ )
	  node_indices.push_back( i );

	// next we calculate the number of nodes there would be if the
	// num_nodes by num_procs schedule grid were perfectly square
	int square_num_nodes = (int)(num_procs + num_procs*floor ((float)(num_nodes-1)/(float)num_procs ));

	for ( i = myid; i < square_num_nodes; i += num_procs )
	{
	
		// get old positions
		get_positions ( node_indices, old_positions );
		
		// default new position is old position
		get_positions ( node_indices, new_positions );
		
		if ( i < num_nodes )
		{

		  // advance random sequence according to myid
		  for ( j = 0; j < 3*myid; j++ )
		    rand();

		  // calculate node energy possibilities
		  if ( !(positions[i].fixed && real_fixed) )
			update_node_pos ( i, old_positions, new_positions );

		  // advance random sequence for next iteration
		  for ( j = 3*myid; j < (int)(3*(node_indices.size()-1)); j++ )
		    rand();

		}
		else
		{
		  // advance random sequence according to use by
		  // the other processors
		  for ( j = 0; j < (int)(3*(node_indices.size())); j++ )
		    rand();
		}
		
		// check if anything was actually updated (e.g. everything was fixed)
		all_fixed = true;
		for ( k = 0; k < node_indices.size (); k++ )
		  if ( !(positions [ node_indices[k] ].fixed && real_fixed) )
		    all_fixed = false;
		  
		// update positions across processors (if not all fixed)
		if ( !all_fixed )
		{
		  // update positions (old to new)
		  update_density ( node_indices, old_positions, new_positions );
		}
		
		// compute node list for next update
		for ( k = 0; k < node_indices.size(); k++ )
		  node_indices [k] += num_procs;
		
		while ( node_indices [ node_indices.size()-1] >= num_nodes )
		  node_indices.pop_back ( );
			
	}
	
	// update first_add and fine_first_add
	first_add = false;
	if ( fineDensity ) fine_first_add = false;
	
}

// The get_positions function takes the node_indices list
// and returns the corresponding positions in an array.

void eb_graph2d::get_positions ( vector<int> &node_indices,
			    float return_positions[3]  )
{
	
	// fill positions
	for(unsigned int i=0; i < node_indices.size(); i++)
	{
		return_positions[3*i] = positions[ node_indices[i] ].x;
		return_positions[3*i+1] = positions[ node_indices[i] ].y;
		//return_positions[3*i+2] = positions[ node_indices[i] ].z;
		return_positions[3*i+2] = 0;
	}
	
}

// update_node_pos -- this subroutine does the actual work of computing
// the new position of a given node.  num_act_proc gives the number
// of active processes at this level for use by the random number
// generators.

void eb_graph2d::update_node_pos ( int node_ind,
			      float old_positions[3],
			      float new_positions[3] )
{	

		double energies[2];			// node energies for possible positions
		float updated_pos[2][3];	// possible positions
		float pos_x, pos_y;
		
		// old VxOrd parameter
		float jump_length = .01f * temperature;
		
		// subtract old node
		density_server.Subtract ( positions[node_ind], first_add, fine_first_add, fineDensity );

		// compute node energy for old solution
		energies[0] = Compute_Node_Energy ( node_ind );

	        // move node to centroid position
		Solve_Analytic ( node_ind, pos_x, pos_y );
		positions[node_ind].x = updated_pos[0][0] = pos_x;
		positions[node_ind].y = updated_pos[0][1] = pos_y;
//		positions[node_ind].z = updated_pos[0][2] = pos_z;

		// Do random method (RAND_MAX is C++ maximum random number)
		updated_pos[1][0] = updated_pos[0][0] + (.5f - rand()/(float)RAND_MAX) * jump_length;
		updated_pos[1][1] = updated_pos[0][1] + (.5f - rand()/(float)RAND_MAX) * jump_length;
//		updated_pos[1][2] = updated_pos[0][2] + (.5f - rand()/(float)RAND_MAX) * jump_length;
		updated_pos[1][2] = 0;
		
		// compute node energy for random position
		positions[node_ind].x = updated_pos[1][0];
		positions[node_ind].y = updated_pos[1][1];
//		positions[node_ind].z = updated_pos[1][2];
		energies[1] = Compute_Node_Energy ( node_ind );
		
		// add back old position
		positions[node_ind].x = old_positions[3*myid];
		positions[node_ind].y = old_positions[3*myid+1];
		if ( !fineDensity && !first_add )
		{
			density_server.Add ( positions[node_ind], fineDensity );
		}
		else if ( !fine_first_add )
		{
			density_server.Add ( positions[node_ind], fineDensity );
		}
		
		// choose updated node position with lowest energy
		if ( energies[0] < energies[1] )
		{
			new_positions[3*myid] = updated_pos[0][0];
			new_positions[3*myid+1] = updated_pos[0][1];
//			new_positions[3*myid+2] = updated_pos[0][2];
			new_positions[3*myid+2] = 0;
			positions[node_ind].energy = (float)energies[0];
		}
		else
		{
			new_positions[3*myid] = updated_pos[1][0];
			new_positions[3*myid+1] = updated_pos[1][1];
//			new_positions[3*myid+2] = updated_pos[1][2];
			new_positions[3*myid+2] = 0;
			positions[node_ind].energy = (float)energies[1];
		}
		
}

// update_density takes a sequence of node_indices and their positions and
// updates the positions by subtracting the old positions and adding the
// new positions to the density grid.

void eb_graph2d::update_density ( vector<int> &node_indices,
			     float old_positions[3],
			     float new_positions[3] )
{
	
	// go through each node and subtract old position from
	// density grid before adding new position
	for ( unsigned int i = 0; i < node_indices.size(); i++ )
	{
		positions[node_indices[i]].x = old_positions[3*i];
		positions[node_indices[i]].y = old_positions[3*i+1];
//		positions[node_indices[i]].z = old_positions[3*i+2];
		density_server.Subtract ( positions[node_indices[i]],
					  first_add, fine_first_add, fineDensity );
		
		positions[node_indices[i]].x = new_positions[3*i];
		positions[node_indices[i]].y = new_positions[3*i+1];
//		positions[node_indices[i]].z = new_positions[3*i+2];
		density_server.Add ( positions[node_indices[i]], fineDensity );
	}	

}

/********************************************
* Function: Compute_Node_Energy			    *
* Description: Compute the node energy		*
* This code has been modified from the      *
* original code by B. Wylie.                *
*********************************************/

double eb_graph2d::Compute_Node_Energy( int node_ind )
{
	
	/* Want to expand 4th power range of attraction */
	float attraction_factor = (float)(attraction*attraction*
			attraction*attraction*2e-2);
	
	map <int,float>::iterator EI;
	double x_dis,y_dis;
	double energy_distance, weight;
	double node_energy=0;
	
	// Add up all connection energies
	for(EI = neighbors[node_ind].begin(); EI != neighbors[node_ind].end(); ++EI) {

		// Get edge weight
		weight = EI->second;
				
		// Compute x,y,z distance
		x_dis = positions[ node_ind ].x - positions[ EI->first ].x;
		y_dis = positions[ node_ind ].y - positions[ EI->first ].y;
		
		// Energy Distance
		energy_distance = x_dis*x_dis + y_dis*y_dis;
		if (STAGE<2) energy_distance *= energy_distance;

		// In the liquid phase we want to discourage long link distances
		if (STAGE==0) energy_distance *= energy_distance;

		node_energy += weight * attraction_factor * energy_distance;
	}

	// add density
	node_energy += density_server.GetDensity ( positions[ node_ind ].x, 
											   positions[ node_ind ].y,
											   fineDensity );

if (node_energy < 0) exit(1);	
	// return computated energy
	return node_energy;
}


/*********************************************
* Function: Solve_Analytic			         *
* Description: Compute the node position     *
* This is a modified version of the function *
* originally written by B. Wylie		     *
*********************************************/

void eb_graph2d::Solve_Analytic( int node_ind, float &pos_x, float &pos_y )
{

   map <int,float>::iterator EI;
   float total_weight = 0;
   float x_dis, y_dis, x_cen=0, y_cen=0;
   float x=0,y=0,dis;
   float damping,weight;

   // Sum up all connections
   for(EI = neighbors[node_ind].begin(); EI != neighbors[node_ind].end(); ++EI) {
		weight = EI->second;
		total_weight += weight;
		x +=  weight * positions[ EI->first ].x;  
		y +=  weight * positions[ EI->first ].y;
   }

   // Now set node position
   if (total_weight > 0) {

		// Compute centriod
		x_cen = x/total_weight;
		y_cen = y/total_weight;
		damping = 1.f - damping_mult;
		pos_x = (float)(damping*positions[ node_ind ].x + (1.0-damping) * x_cen);
		pos_y = (float)(damping*positions[ node_ind ].y + (1.0-damping) * y_cen);
   }
   
   // No cut edge flag (?)
   if (min_edges == 99) return;

   // Don't cut at end of scale
   if ( CUT_END >= 39500 ) return;

   float num_connections = (float)sqrtf( (float)neighbors[node_ind].size() );
   float maxLength = 0;

   map<int, float>::iterator maxIndex;

   // Go through nodes edges... cutting if necessary
   for(EI = maxIndex = neighbors[node_ind].begin();
	   EI !=neighbors[node_ind].end(); ++EI) {

		// Check for at least min edges
		if (neighbors[node_ind].size() < min_edges) continue;

		x_dis = x_cen - positions[ EI->first ].x;
		y_dis = y_cen - positions[ EI->first ].y;
		dis = x_dis*x_dis+y_dis*y_dis;
		dis *= num_connections;

		// Store maximum edge
		if (dis > maxLength) {maxLength = dis; maxIndex=EI;}
   }

   // If max length greater than cut_length then cut
   if (maxLength > cut_off_length) neighbors[ node_ind ].erase( maxIndex ); 
   
}


// write_gvl writes out the gvl file of the final solutions

void eb_graph2d::write_gvl( const char *file_name, const char *wel_file_name, float thresh, int absValFlag, int attractFlag, float grav, int springAlgo)
{
  unsigned int i;
  ofstream gvlOUT( file_name );
  if ( !gvlOUT )
  {
	cout << "Could not open " << file_name << ".  Program terminated." << endl;
	  exit (1);
  }
  
  cout << "Writing out solution to " << file_name << " ..." << endl;
  
  float max_pos = 0;
  for (i = 0; i < positions.size(); i++) {
	  if ((float)fabs(positions[i].x) > max_pos)
	  {
		  max_pos = (float)fabs(positions[i].x);
	  }
	  if ((float)fabs(positions[i].y) > max_pos)
	  {
		  max_pos = (float)fabs(positions[i].y);
	  }
  }
  
  float d_val = 1; 
  if (max_pos != 0)
  {
	  d_val = max_pos / 29999;
  }

  // place metadata in gvl
  gvlOUT << wel_file_name << "\t" << thresh << "\t" << absValFlag << "\t" << attractFlag
	  << "\t" << grav << "\t" << springAlgo << endl;
  gvlOUT << num_vertices << "\t" << num_edges << endl;
  
  for (i = 0; i < positions.size(); i++) {
	  // WES
	  positions[i].x /= d_val;
	  positions[i].y /= d_val;
	  
    gvlOUT << vertex_labels[positions[i].id] << "\t" << (int)positions[i].x << "\t" << (int)positions[i].y << "\t" << 0.0 << endl;
  }
  gvlOUT.close();
  
}

// get_tot_energy adds up the energy for each node to give an estimate of the
// quality of the minimization.

float eb_graph2d::get_tot_energy ( )
{

	float my_tot_energy, tot_energy;
	my_tot_energy = 0;
	for ( int i = myid; i < num_nodes; i += num_procs )
	  my_tot_energy += positions[i].energy;
	  
	tot_energy = my_tot_energy;
	
	return tot_energy;
	
}


// The following subroutine draws the graph with possible intermediate
// output (int_out is set to 0 if not proc. 0).  int_out is the parameter
// passed by the user, and gvl_file is the .gvl file.

void eb_graph2d::draw_graph ( int int_out, char *gvl_file )
{
	
	// layout graph (with possible intermediate output)
	int count_iter = 0, count_file = 1;
	//char int_gvl_file [254];
	while ( ReCompute( ) )
	{}
	/*
		if ( (int_out > 0) && (count_iter == int_out) )
		{
			// output intermediate solution
			sprintf ( int_gvl_file, "%s.%d", gvl_file, count_file );
			write_gvl ( int_gvl_file );
			
			count_iter = 0;
			count_file++;
		}
		else
			count_iter++;*/
	
}
