// IncidenceMatrix.H
// JR 9/8/95,9/16/95,9/17/95
//	10/29/95  Reviewed in light of changes to lower-level classes.
//	11/3/95   Compiled.
//      11/4/95   Tested constructor.
//	11/8/95   Testing eigenstuff.
//		  One example (5x5) tried. Got answer agreeing
//		  with Mathematica to 4 digits when tol=.00001.	  


#ifndef INCIDENCEMATRIX_H
#define INCIDENCEMATRIX_H

#include <assert.h>
#include "Map.H"
#include "Set.H"
#include "LinList.H"
#include "MatrixD.H"
#include "MatrixI.H"
#include "Boolean.H"

class IncidenceMatrix 
{
   public:
        IncidenceMatrix()	// Use not permitted.
	{  assert(0);	} 
       	IncidenceMatrix( const Set<Name>& names, const Map& m );
	IncidenceMatrix(const IncidenceMatrix& m);  // copy constructor       
       ~IncidenceMatrix()
	{}
	int size() const { return mi.rows(); }
	const MatrixI& matrix() { return mi; }
        Boolean reducible() const;
	Boolean isPermutation() const;    
        double eigenvectorComponent(const Name& n,double tol,int maxIts=0);
	double  growthRate(double tol,int maxIts=0);

   private:
        LinList<Name> sources;
	MatrixI mi;
	MatrixD v;	// eigenvector
	double  eigenvalue;
	double  currentTol;

	void compute_eigenstuff(double tol,int maxIts);
	Boolean addDestinations(Set<int>& s) const;
};

           						  // constructor
IncidenceMatrix::IncidenceMatrix(const Set<Name>& names,const Map& m)
: currentTol(1.0e30) // a huge number
{
				// make a list from the set of names
        SetBrowser<Name> i(names);
        for(i.init();i.ok();i.advance()) sources.append(i.current());
        int n=sources.length();
        mi = MatrixI(n);
       
        LinListBrowser<Name> js(sources),jc(sources);

	int njs;
        for(js.init(),njs=0;js.ok();++js,++njs)   // for each edge name 
		  		      // make a column of M_H       
        {       
                DELListBrowser k(m.imageOf(DE(js.ahead(),0)));
		int njc;
		for(jc.init(),njc=0;jc.ok();++jc,++njc)
	                for(k.init();k.ok();++k) 
        	                if(k.ahead().name()==jc.ahead())
                	                mi(njc,njs)+=1;  // using what indexing??
        }
//cout << mi << endl;

 }

IncidenceMatrix::IncidenceMatrix(const IncidenceMatrix& m)  // copy constructor
:  sources(m.sources), mi(m.mi)
{}

double IncidenceMatrix::eigenvectorComponent
				(const Name& n,double tol,int maxIts)
{
     	// Return component of positive eigenvector approximation.

	assert( tol > 0.0 );
	if( tol < currentTol)  // need to recompute eigenstuff
	{
		compute_eigenstuff(tol,maxIts);
		currentTol = tol;
	}
				// find index, jname, of Name n
	LinListBrowser<Name> i(sources);
	int jname=0;
	for(i.init();i.ok();++i,jname++) 
		if(i.ahead()==n)  return v(jname,0);

	assert(0);  // should never get to here 
}

double IncidenceMatrix::growthRate(double tol,int maxIts)
{
     	// Return approximation of the eigenvalue of positive eigenvector.
	assert( tol > 0.0 );
	if( tol < currentTol )
	{
		compute_eigenstuff(tol,maxIts);
		currentTol = tol;
	}
	return eigenvalue;	
}

void IncidenceMatrix::compute_eigenstuff(double tol,int maxIts)
{
        // Uses the Power Method with a shift of +1.
	// A positive integer matrix has a unique positive eigenvector
	// with eigenvalue equal to the spectral radius.
	// All other eigenvalues have strictly smaller modulus, but
	// here we shift by +1 to enhance the dominance of the largest.

// Is irreducibility necessary for the above?
// Do we need to 

	if(reducible()) cout << "WARNING: MH is reducible." << endl;
//	assert(!reducible());

// ? See Bestvina&Handel bottom of p110

	v = MatrixD(mi.rows(),1,1.0);
				 // Note: this cannot be orthogonal
			    	 // to the (positive) eigenvector.
	MatrixD vOld(mi.rows(),1,1.0e30);
	MatrixD mis(mi);
	for(int i=0;i<mi.rows();i++) mis(i,i) += 1.0;  // add identity
       	int nIterations=0;
	double eval;
       	while(     ((v-vOld).largestMag() > tol) 
		&& ( nIterations < maxIts))
       	{
		nIterations++;
               	vOld=v;
               	v = mis*v;
		eval=v.largestMag();
		v/=eval;
//cout << v << endl;
       	}
//cout << "  Finished eigenvector estimation in " << nIterations << " iterations." << endl;
//cout << v << endl;
	if(nIterations>=maxIts)
	{
cout << "WARNING: Failed to converge on eigenvector." << endl;
	}
	eigenvalue = eval - 1.0;	// compensate for +1 shift
}

Boolean IncidenceMatrix::isPermutation() const
{
    // See if there is exactly one 1 in each row and column.
    int m=mi.rows();
    for(int i=0;i<m;i++)
    {
       int count=0;
       for(int j=0;j<m;j++)
       {
          if(mi(i,j)>1) return 0;
          if(mi(i,j)==1) count++;
       }
       if(count!=1) return 0;

       count=0;
       for(j=0;j<m;j++)
       {
          if(mi(j,i)>1) return 0;
          if(mi(j,i)==1) count++;
       }
       if(count!=1) return 0;
    }
    return 1;	    // Matrix is a permutation matrix.
}


Boolean IncidenceMatrix::reducible() const
{
   // Decide if the incidence matrix, is reducible or irreducible
   // (under symmetric permutations). 
//cout << "Entered IncidenceMatrix::reducible(), size of MH is " << mi.rows()<< endl;
   const Boolean REDUCIBLE=1;
   const Boolean IRREDUCIBLE=0;

   for(int i=0;i<mi.rows();i++)  // Determine "strong component" of edge i.
   {   	Set<int> s(i);
       	while(addDestinations(s))
	{}
//cout << "Edge number " << i << " has strong component " << s << endl;	
       	if(s.cardinality()<mi.rows()) 
	{	
cout << "  M_H is reducible:" << endl;
cout << mi << endl;
		return REDUCIBLE;
	}
   }
   
//cout << "  M_H is irreducible." << endl;
   return IRREDUCIBLE;
}

Boolean IncidenceMatrix::addDestinations(Set<int>& s) const
{
   Boolean somethingAdded=0;
   SetBrowser<int> si(s);
   for(si.init();si.ok();si.advance())
   {   int j=si.current();
       for(int i=0;i<mi.rows();i++)
          if( mi(i,j) ) if(s.add(i)) somethingAdded=1;
   }
   return somethingAdded;
}

double Homeo::growthRate(double tol,int maxIts) const
{
	Set<Name> H(theNameManager.edges());
	Set<Name> P(theNameManager.periphery());
	Set<Name> prePee(preP());
	H -= P;
	H -= prePee;
cout << "  Computing growth rate." << endl;
cout << "    pre-P is " << prePee << "." << endl;
cout << "    H is " << H << "." << endl;
	IncidenceMatrix MH( H, theMap );
	return MH.growthRate(tol,maxIts);
}													



#endif