//Graph Utilities
#include "graph.h"

uint graph::FindEdge(long Label)
{
	for (uint i=1; long(i)<=Edges.TopIndex(); i++)
		if (Edges[i].Label == Label || Edges[i].Label == -Label) return i;
	return 0;
}

uint graph::FindVertex(uint Label)
{
	for (uint i=1; long(i)<=Vertices.TopIndex(); i++) if (Vertices[i].Label == Label) return i;
	return 0;
}

long graph::Derivative(long Label)
{
	uint Index = FindEdge(Label);
	if (!Edges[Index].Image.TopIndex()) return 0;
	if (Label>0) return (Edges[Index].Image[1]);
	return (-Edges[Index].Image[Edges[Index].Image.TopIndex()]);
}

long graph::Derivative(long Label, uint n)
{
	for (uint i=1; i<=n; i++)
	{
		Label = Derivative(Label);
		if (!Label) return 0;
	}
	return Label;
}

long graph::AltDerivative(long Label, uint n)
{
	for (uint i=1; i<=n; i++)
	{
		uint Index = FindEdge(Label);
		intarray Image = Edges[Index].Image;
		if (Label<0) Image.Invert();
		uint j=1;
		while (long(j)<=Image.TopIndex() && IsPeripheral(Image[j])) j++;
		if (long(j)>Image.TopIndex()) return 0;
		Label = Image[j];
	}
	return Label;
}

void graph::Replace(long Label, intarray& L)
{
	edgeiterator I(Edges);
	do
	{
		((I++).Image).Replace(Label, L);
	} while (!I.AtOrigin());
}

void graph::RemoveAll(long Label)
{
	edgeiterator I(Edges);
	do
	{
		((I++).Image).RemoveAll(Label);
	} while (!I.AtOrigin());
}

bool graph::Tighten()
{
	bool Result = false;
	edgeiterator I(Edges);
	do
	{
	  if (((I++).Image).Tighten()) Result = true;
	} while (!I.AtOrigin());
	return Result;
}

decimal graph::Growth()
{
	matrix M(*this);
	return M.GrowthRate();
}

bool graph::IntersectsP(long Label)
{
	uint Index = FindEdge(Label);
	edge& Now = Edges[Index];
	if (Now.Type == Peripheral) return true;
	return (OnP(Now.Start) || OnP(Now.End));
}

uint graph::OnPInd(uint Index)
{
   intiterator I(Vertices[Index].Edges);
   do
   {
   	uint EIndex = FindEdge(I++);
   	if (Edges[EIndex].Type == Peripheral) return Edges[EIndex].Puncture;
  } while (!I.AtOrigin());
   return 0;
}

uint graph::OnP(uint Label)
{
	uint Index = FindVertex(Label);
	intiterator I(Vertices[Index].Edges);
	do
	{
		uint EIndex = FindEdge(I++);
		if (Edges[EIndex].Type == Peripheral) return Edges[EIndex].Puncture;
	} while (!I.AtOrigin());
	return 0;
}

bool graph::IsPeripheral(long Label)
{
	uint Index = FindEdge(Label);
	return (Edges[Index].Type == Peripheral);
}

void graph::FindTypes()
{
	if (!Punctures)
	{
		for (uint i=1; i<=NumberEdges(); i++) Edges[i].Type = Main;
		return;
	}
	edgeiterator I(Edges);
	//First Flag all peripheral edges
	do
	{
		edge& Now = I++;
		Now.Flag = (Now.Type == Peripheral) ? true : false;
	} while (!I.AtOrigin());
	bool Changed = true;
	while (Changed)
	{
		Changed = false;
		do
		{
			edge& Now = I++;
			if (Now.Flag) continue;
			intiterator J(Now.Image);
			bool ShouldFlag = true;
			do
			{
				uint Index = FindEdge(J++);
				if (!Edges[Index].Flag)
				{
					ShouldFlag = false;
					continue;
				}
			} while (!J.AtOrigin());
			if (ShouldFlag)
			{
				Now.Flag = true;
				Now.Type = Preperipheral;
				Changed = true;
			}
		} while (!I.AtOrigin());
	}
}

uint graph::From(long Label)
{
	uint Index = FindEdge(Label);
	return ( (Label>0) ? Edges[Index].Start : Edges[Index].End );
}

void graph::Flush()
{
    uint i;
	for (i=1; i<=NumberEdges(); i++) Edges[i].Image.Flush();
	for (i=1; i<=NumberVertices(); i++) Vertices[i].Edges.Flush();
	Edges.Flush();
	Vertices.Flush();
	Turns.Flush();
	NextVertexLabel = 1; NextEdgeLabel = 1; Punctures = 0;
	Type = Unknown;
}


bool graph::IsProperSubForest(bool* Inset, uint n)
{
	bool Proper = false;      //First check subgraph is proper
    uint i;
	for (i=1; i<=n; i++) if (!Inset[i]) Proper = true;
	if (!Proper) return false;
	uint m = NumberVertices();
	bool *VertSet = new bool[m+1], *Changed = new bool[m+1], *NewChanged = new bool[m+1];
	bool Result = true;
	for (i=1; i<=n; i++) if (Inset[i]) //Look for loop in subgraph containing i
	{
		long Label = -Edges[i].Label; //Don't use this edge in making loop
		uint FromVertex = Edges[i].End; //Start here
		uint AimVertex = Edges[i].Start; //Aim to arrive here
        uint j;
		for (j=1; j<=m; j++) VertSet[j] = Changed[j] = NewChanged[j] = false;
		VertSet[FindVertex(FromVertex)] = Changed[FindVertex(FromVertex)] = true;
		bool SomeChanged = true;
		while (SomeChanged)
		{
			for (j=1; j<=m; j++) if (Changed[j])
			{
				intarray& Now = Vertices[j].Edges;
				intiterator I(Now);
				do
				{
					long NowLabel = (I++);
					uint Index = FindEdge(NowLabel);
					if (NowLabel == Label || !Inset[Index] || Index<i) continue; //Now have legitimate edge
					uint NewVertex = FindVertex(To(NowLabel));
					if (To(NowLabel)==AimVertex)
					{
						Result = false;
						break;
					}
					if (!VertSet[NewVertex]) VertSet[NewVertex] = NewChanged[NewVertex] = true;
				} while (!I.AtOrigin());
				if (!Result) break;
			}
			if (!Result) break;
			SomeChanged = false;
			for (j=1; j<=m; j++)
			{
				Changed[j] = NewChanged[j];
				if (NewChanged[j]) SomeChanged = true;
				NewChanged[j] = false;
			}
		}
		if (!Result) break;
	}
	delete[] VertSet;
	delete[] Changed;
	delete[] NewChanged;
	return Result;
}

void graph::IterateTurn(turn& T)
{
	T.i = Derivative(T.i);
	T.j = Derivative(T.j);
}

turn graph::FindTurns()
{
	Turns.Flush();
	turn Result, Current;
	turnlist T;
	Result.Level = 0;
	edgeiterator I(Edges);
	do
	{
		intarray& Image = (I++).Image;
		for (uint k=1; long(k)<Image.TopIndex(); k++)
		{
			Current.i = -Image[k];
			Current.j = Image[k+1];
			long Found;
			if ((Found = Turns.Find(Current)) != -1)
			{
				Current.Level = Turns[Found].Level;
				if (Current.Level && (!Result.Level || Current.Level<Result.Level)) Result = Current;
				continue;
			}
			T.Flush();
			T[1] = Current;
			bool Finished = false;
			while (!Finished)
			{
				IterateTurn(Current);
				Found = Turns.Find(Current);
				if (Found != -1)
				{
					Finished = true;
					uint FoundLevel = Turns[Found].Level;
					if (!FoundLevel) for (uint l=1; long(l)<=T.TopIndex(); l++) T[l].Level = 0;
					else for (uint l=1; long(l)<=T.TopIndex(); l++) T[l].Level = FoundLevel+T.TopIndex()+1-l;
				}
				else
				{
					if (T.Find(Current) != -1)
					{
						Finished = true;
						for (uint l=1; long(l)<=T.TopIndex(); l++) T[l].Level = 0;
					}
					else
					{
						if (Current.IsDegenerate())
						{
							Finished = true;
							for (uint l=1; long(l)<=T.TopIndex(); l++) T[l].Level = T.TopIndex()+1-l;
						}
						else T.Add(Current);
					}
				}
			}
			Turns.Append(T);
			if (T[1].Level && (!Result.Level || T[1].Level<Result.Level)) Result = T[1];
		}
	} while (!I.AtOrigin());
	return Result;
}


ostream& operator<<(ostream& Out, turn T)
{
	Out << "Turn " << T.i << ", " << T.j << " at level " << T.Level;
	return (Out);
}

turnplace graph::LocateTurn(turn& T)
{
	turnplace Result;
	edgeiterator I(Edges);
	do
	{
		intarray& Now = (I.Now()).Image;
		for (uint i=1; long(i)<Now.TopIndex(); i++)
			if ( (Now[i]==-T.i && Now[i+1]==T.j) || (Now[i]==-T.j && Now[i+1]==T.i) )
			{
				Result.Label = (I.Now()).Label;
				Result.Position = i;
				return Result;
			}
		I++;
	} while (!I.AtOrigin());
	THROW("Locating non-existent turn",1);
}







bool graph::HasIrreducibleMatrix()
{
	matrix M(*this);
	return M.IsIrreducible();
}

bool graph::RetractsOntoP(bool* Inset, uint n)
{
// Need to test whether we can get from one peripheral loop to another in Inset.
//First test whether Inset minus one edge from each peripheral loop is a proper subforest.
	bool *PunctureDone = new bool[Punctures+1];
    uint i;
	for (i=1; i<=Punctures; i++) PunctureDone[i] = false;
	for (i=1; i<=n; i++) if (IsPeripheral(i) && !PunctureDone[Edges[i].Puncture])
	{
		Inset[i] = false;
		PunctureDone[Edges[i].Puncture] = true;
	}
	delete[] PunctureDone;
	if (!IsProperSubForest(Inset, n)) return false;
	for (i=1; i<=n; i++) if (IsPeripheral(i)) Inset[i] = true; //Restore peripheral subgraph
	uint m = NumberVertices();
	bool *VertSet = new bool[m+1], *Changed = new bool[m+1], *NewChanged = new bool[m+1], Result = true;
	for (i=1; i<=m; i++) if (OnP(Vertices[i].Label))
	{
	//Determine set of vertices connected to Vertices[i]
		uint FromPuncture = OnP(Vertices[i].Label); //Start from this peripheral loop
        uint j;
		for (j=1; j<=m; j++) VertSet[j] = Changed[j] = NewChanged[j] = false;
		VertSet[i] = Changed[i] = true;
		bool SomeChanged = true;
		while (SomeChanged)
		{
			for (j=1; j<=m; j++) if (Changed[j])
			{
				intarray& Now = Vertices[j].Edges;
				intiterator I(Now);
				do
				{
					long NowLabel = (I++);
					uint Index = FindEdge(NowLabel);
					if (!Inset[Index]) continue;
					uint NewVertex = To(NowLabel);
					if (OnP(NewVertex) && (OnP(NewVertex) != FromPuncture) )
					{
						Result = false;
						break;
					}
					uint NewVertIndex = FindVertex(NewVertex);
					if (!VertSet[NewVertIndex]) VertSet[NewVertIndex] = NewChanged[NewVertIndex] = true;
				} while (!I.AtOrigin());
				if (!Result) break;
			}
			if (!Result) break;
			SomeChanged = false;
			for (j=1; j<=m; j++)
			{
				Changed[j] = NewChanged[j];
				if (NewChanged[j]) SomeChanged = true;
				NewChanged[j] = false;
			}
		}
		if (!Result) break;
	}
	delete[] VertSet; delete[] Changed; delete[] NewChanged;
	return Result;
}

bool graph::NeedToAbsorb()
{
	edgeiterator I(Edges);
	do
	{
		edge& Now = I++;
		if (Now.Type == Peripheral) continue;
		if (OnP(Now.Start) && IsPeripheral(Now.Image[1])) return true;
		if (OnP(Now.End) && IsPeripheral(Now.Image[Now.Image.TopIndex()])) return true;
	} while (!I.AtOrigin());
	return false;
}

bool graph::Collapses(intarray& L)
{
	intarray M;
	for (uint i=1; i<=2*NumberEdges(); i++) //INEFFICIENT
	{
		M.Flush();
		intarray Image = Edges[FindEdge(L[1])].Image;
		if (L[1]<0) Image.Invert();
		uint j = Image.TopIndex(); while (IsPeripheral(Image[j])) j--;
		while (long(j)<=Image.TopIndex()) M[M.TopIndex()+1] = Image[j++];
		for (j=2; long(j)<L.TopIndex(); j++) //Peripheral edges
		{
			Image = Edges[FindEdge(L[j])].Image;
			if (L[j]<0) Image.Invert();
			M.Append(Image);
		}
		Image = Edges[FindEdge(L[L.TopIndex()])].Image;
		if (L[L.TopIndex()]<0) Image.Invert();
		j=1; while (IsPeripheral(Image[j])) j++;
		for (uint k=1; k<=j; k++) M[M.TopIndex()+1] = Image[k];
		M.Tighten();
		if (!M.TopIndex()) return true; 
		L = M;
	}
	return false;
}

void graph::FindReduction()
{
	if (!Type == Reducible1) THROW("Trying to find reduction with no invariant subgraph",1);
	Reduction.Flush();
	bool* Indices = new bool[NumberEdges()];
	matrix M(*this);
	M.IsIrreducible(Indices);
	uint i=0;
	for (uint j=1; j<=NumberEdges(); j++)
		if (Edges[j].Type == Main)
			if (Indices[i++]) Reduction.SureAdd(Edges[j].Label);
	delete[] Indices;
}

void graph::FindGates()
{
	if (Type != pA_or_red) THROW("Trying to find gates with non-efficient graph map",1);
	Reduction.Flush(); //Will hold gates in format (0,vert label,0,gate1,0,gate2,0,..,gaten,00,vertlabel..)
	Turns.Flush(); //Will hold infinitesimal edges i,j give gates, Level gives vertex label
	//Determine gates at each vertex in turn
    uint i;
	for (i=1; i<=NumberVertices(); i++)
	{
		Reduction.SureAdd(0); Reduction.SureAdd(0);
		intarray& Now = Vertices[i].Edges;
		intarray Class;
        uint j;
		for (j=1; long(j)<=Now.TopIndex(); j++) Class[j] = Derivative(Now[j], 2*NumberEdges());
		Reduction.SureAdd(Vertices[i].Label); Reduction.SureAdd(0);
		uint Start = 2; while (Class[Start] == Class[1]) Start++;
		j=Start;
		do
		{
			Reduction.SureAdd(Now[j]);
			uint Next = (long(j)==Now.TopIndex()) ? 1 : j+1;
			if (Class[j] != Class[Next]) Reduction.SureAdd(0);
			j = Next;
		} while (j!= Start);
	}
	Reduction.SureAdd(0);
	//Determine all infinitesimal edges connecting gates
	for (i=1; i<=NumberEdges(); i++)
	{
		intarray& Now = Edges[i].Image;
		if (Now.TopIndex()==1) continue;
		for (uint j=1; long(j)< Now.TopIndex(); j++)
		{
			bool Continue = true;
			long Label1 = -Now[j], Label2 = Now[j+1];
			while (Continue)
			{
				uint Vertex = From(Label2);
				//Find Vertex in gate list
				uint k=3; while ((Reduction[k] != long(Vertex)) || (Reduction[k-1]) || (Reduction[k-2])) k++;
				//Find canonical labels of  gates
				uint l=k+2; while (Reduction[l] != Label1) l++; while (Reduction[l-1]) l--;
				Label1 = Reduction[l];
				l=k+2; while (Reduction[l] != Label2) l++; while (Reduction[l-1]) l--;
				Label2 = Reduction[l];
				if (Label2<Label1)
				{
					long Temp = Label2; Label2 = Label1; Label1 = Temp;
				}
				turn T; T.i = Label1; T.j = Label2; T.Level = Vertex;
				Continue = Turns.Add(T);
				if (Continue)
				{
					Label1 = Derivative(Label1); Label2 = Derivative(Label2);
				}
			}
		}
	}
}

void graph::FindTrack()
{
   if (Type == pA || Type == Reducible2) return; //Should already have calculated
	FindGates();
	//Determine connectivity at each vertex
	Type = pA; //Change to reducible2 if any vertex is not connected
	for (uint i=1; i<=NumberVertices(); i++)
	{
		uint Count = 0;
		for (uint j=1; long(j)<=Turns.TopIndex(); j++) if (Turns[j].Level == Vertices[i].Label) Count++;
		//Find Vertex in gate list
		uint k=3; while (Reduction[k] != long(Vertices[i].Label) || Reduction[k-1]  || Reduction[k-2]) k++;
		//Count gates at vertex
		uint NumberGates = 0;
		k+=3;
		while (true)
		{
			if (!Reduction[k])
			{
				NumberGates++;
				if (!Reduction[++k]) break;
			}
			k++;
		}
		if (Count < NumberGates-1)
		{
			Type = Reducible2;
			Vertices[i].Flag = false;
		}
		else Vertices[i].Flag = true;
	}
}

/*void graph::TampDown()
{
   bool Finished = false;
   bool Touched = false;
   int Count = 0;
   while (!Finished)
   {
      Finished = true;
      Count++;
      if (Count > 10000) THROW("Tamping not converging", 4);
      //Ensure all vertices non-flagged
      for (uint i=1; i<=NumberVertices(); i++) Vertices[i].Flag = false;
      // Look for a non-peripheral vertex with two adjacent edges which have same derivative
      for (uint i=1; i<=NumberVertices(); i++)
      {
         if (OnPInd(i)) continue;
         intarray& E = Vertices[i].Edges;
         for (uint j=1; long(j)<=E.TopIndex(); j++)
         {
            uint k = (long(j)==E.TopIndex()) ? 1 : j+1;
            if (Derivative(E[j]) == Derivative(E[k]))
            {
               Finished = false;
               Touched = true;
               FoldAsMuchAsPossible(E[j],E[k]);
               AbsorbIntoP();
               PullTight();
               break;
            }
         }
         if (!Finished) break;
      }
   }
   if (Touched)
   {
      //Just to be sure...
      FindTrainTrack();
      FindTrack();
   }
}*/


void graph::TampDown()
{
   bool Finished = false;
   bool Touched = false;
   int Count = 0;
   while (!Finished)
   {
      Finished = true;
      Count++;
      if (Count > 10000) THROW("Tamping not converging", 4);
      //Ensure all vertices non-flagged
      for (uint i=1; i<=NumberVertices(); i++) Vertices[i].Flag = false;
      // Look for a non-peripheral vertex with two adjacent edges which have same derivative
      // Or for a peripheral vertex with three adjacent edges which have same derivative
      long Label1, Label2; uint Steps = 0;
      for (uint i=1; i<=NumberVertices(); i++)
      {
         if (OnPInd(i))
         {
 /*           //Look for 3 adjacent. Make sure Label1 and Label2 bracket the third
            intarray& E = Vertices[i].Edges;
            if (E.TopIndex()<5) continue;
            for (uint j=1; long(j)<=E.TopIndex(); j++)
            {
               uint k = j+2; if (long(k)>E.TopIndex()) k-=E.TopIndex();
               long D1 = Derivative(E[j],2*NumberEdges()), D2 = Derivative(E[k],2*NumberEdges());
               if (D1 == D2)
               {
                  uint l = (long(j)==E.TopIndex()) ? 1 : j+1;
                  if (D1 == Derivative(E[l],2*NumberEdges()))
                  {
                     Finished = false;
                     Touched = true;
                     long L1 = E[j], L2 = E[k], L3 = E[l];
                     uint s = 1;
                     while ((Derivative(L1)!=Derivative(L2)) || (Derivative(L1)!=Derivative(L3))) {s++; L1=Derivative(L1); L2=Derivative(L2); L3=Derivative(L3);}
                     if (Steps == 0 || s<Steps)
                     {
                        Steps=s; Label1=L1; Label2=L2;
                     }
                  }
               }
            }  */
         }
         else
         {
            intarray& E = Vertices[i].Edges;
            for (uint j=1; long(j)<=E.TopIndex(); j++)
            {
               uint k = (long(j)==E.TopIndex()) ? 1 : j+1;
               if (Derivative(E[j],2*NumberEdges()) == Derivative(E[k],2*NumberEdges()))
               {
                  Finished = false;
                  Touched = true;
                  long L1 = E[j], L2 = E[k];
                  uint s = 1;  while (Derivative(L1)!=Derivative(L2)) {s++; L1=Derivative(L1); L2=Derivative(L2);}
                  if (Steps == 0 || s<Steps)
                  {
                     Steps = s; Label1 = L1; Label2 = L2;
                  }
               }
            }
         }
      }
      if (!Finished)
      {
         FoldAsMuchAsPossible(Label1, Label2, false);
         AbsorbIntoP();
         PullTight();
      }
   }
   if (Touched)
   {
      //Just to be sure...
      FindTrainTrack();
      FindTrack();
   }
}