// runBH.c (from specified "Homeo").

//#include <iostream.h>
#include <iostream>
using namespace std;
#//	3/7/96	Compiled without any errors or warnings.
//		This ran to completion on one of the homeo*.dat files.
//		After making a BUNCH of changes from call by reference
//		(const DE&) to call by value (DE), the code now crashes
//		the SGI compiler:
//		ugen: internal error writing to /tmp/ctmca002ea: Error 0
//		and gives the following message on the Sparc:
// 		memory allocation exceeded
//		Something is wrong! Go back to testing individual classes.
//		Problem was just that the directory /tmp was getting full
//		during compilation.
//	
//	3/17/96 Code has been running for a week now.
//		Problem has been that the algorithm is going round in a
//		circle: after doing the collapsing after step 7, we
//		get back exactly what we had before step 7.
//		This is not a bug in the code: if I follow the algorithm
//		by hand, I get the same thing.
//	4/20/96 Changed from main() to member function of Homeo.
//	4/22/96	Changed call of step6 from repeated to single.
//	5/23/96 Added second kind of reducibility check: reducibleType2().


#include "Homeo.H"
#include <float.h>	// For DBL_MAX.
#define VERBOSE

const int PERIODIC=0;
const int REDUCIBLE=1;
const int PSEUDOANOSOV=2;

int Homeo::BH()	// Perform the Bestvina-Handel algorithm.
{

    double newGrowth,oldGrowth=DBL_MAX; // a number that should exceed initial growth rate

    if(!componentsOfPeripheryAreCircles())
    {
	  cout << "Cannot handle graph whose peripheral components are not circles." << endl;
	  cout << (*this) << endl;
	  return -1;
    }
    assert(mapIsContinuous());	

newGrowth=growthRate(1.0e-7,100);
cout << "    Growth rate (which was " << oldGrowth <<
") is now " << newGrowth << "." << endl;
//SIZE();
//assert( newGrowth <= oldGrowth );
oldGrowth=newGrowth;

    pullTight();

    cout << endl << (*this) << endl;
    assert(mapIsContinuous());

int counter=0;

newGrowth=growthRate(1.0e-7,100);
cout << "    Growth rate (which was " << oldGrowth <<
") is now " << newGrowth << "." << endl;
//SIZE();
//assert( newGrowth <= oldGrowth );
oldGrowth=newGrowth;


    while(1)
    {
cout << ++counter << endl;
	Set<Name> invariantForest=invariantForestNotIntersectingPInAnEdge();
	if( invariantForest.cardinality()
		==theNameManager.edges().cardinality()
		 - theNameManager.periphery().cardinality() )
	{
		// Collapse to bouquet and stop.
		
cout << "  The entire complement of P, " << invariantForest 
	<< ", is invariant. Collapsing it." << endl; 

		SetIterator<Name> fi(invariantForest);
		for( fi.init(); fi.ok(); fi.advance() ) 
			collapseEdge(DE(fi.current(),0));
cout << endl << (*this) << endl;
cout << "This is periodic." << endl;
		return PERIODIC;
	}



	if( ! collapseInvariantForest())
	{
assert(mapIsContinuous());
	    if( ! removeValenceOneVertices() )
	    {
assert(mapIsContinuous());
	        if(!absorbIntoP())
	        {
assert(mapIsContinuous());
		    if( reducible() )
		    {
cout << endl << (*this) << endl;

assert(mapIsContinuous());
newGrowth=growthRate(1.0e-7,100);
cout << "    Growth rate (which was " << oldGrowth <<
") is now " << newGrowth << "." << endl;
//SIZE();
//assert( newGrowth <= oldGrowth );
oldGrowth=newGrowth;

		    	cout << "This is reducible." << endl;
		    	return REDUCIBLE;
		    }
		    else if ( reducibleType2() )
		    {
cout << endl << (*this) << endl;
		    	cout << "This is reducible (type-2)." << endl;
		    	return REDUCIBLE;
		    }
		    else
		    {
cout << "  M_H is irreducible." << endl;
			if( ! removeAValenceTwoVertex() )
			{
assert(mapIsContinuous());
			    Boolean didAFold=0;
	           	    didAFold = step6();
assert(mapIsContinuous());
			    if( ! didAFold )
			    {
				if( efficient() )
				{

assert(mapIsContinuous());
newGrowth=growthRate(1.0e-7,100);
cout << "    Growth rate (which was " << oldGrowth <<
") is now " << newGrowth << "." << endl;
//SIZE();
//assert( newGrowth <= oldGrowth );
oldGrowth=newGrowth;
cout << endl << (*this) << endl;
		   		    cout << "This is efficient." << endl;
		     		    return PSEUDOANOSOV;
				}
				else
				{
cout << endl << (*this) << endl;
				    step7();
				}
			    }
			}
		    }
		}
	    }
	}

assert(mapIsContinuous());
newGrowth=growthRate(1.0e-7,100);
cout << "    Growth rate (which was " << oldGrowth <<
") is now " << newGrowth << "." << endl;
//SIZE();
//assert( newGrowth <= oldGrowth );
oldGrowth=newGrowth;

	pullTight();

cout << endl << (*this) << endl;
assert(mapIsContinuous());
newGrowth=growthRate(1.0e-7,100);
cout << "    Growth rate (which was " << oldGrowth <<
") is now " << newGrowth << "." << endl;
//SIZE();
//assert( newGrowth <= oldGrowth );
oldGrowth=newGrowth;

    }

}

