package dmrgvisu2;

public class Mps
{
	private static final long serialVersionUID = 1L;
	int N;
	Matrix[][] A;
	Matrix[][] L;
	Matrix[][] R;
	int d;
	int d_max;
	int states;

	public Mps()
	{
		System.out.println("DNU");
		this.d=2;									// Maximale Groesse der Matrizen
		this.d_max=1024;		
		this.N=100000;									// Laenge der Kette
		this.states=2;								// Anzahl der moeglichen Zustaende (Anzahl der Matrizen pro Kettenplatz)
		this.A = new Matrix[this.N][this.states];
		this.R = new Matrix[this.N][this.d_max];
		this.L = new Matrix[this.N][this.d_max];
		for (int i=0; i<this.N/2; i++)
		{
			int d_x = (int)(Math.pow(d,i+1)<d_max?Math.pow(d,i+1):d_max);
			int d_y = (int)(Math.pow(d,i)<d_max?Math.pow(d,i):d_max);
			for (int j=0; j<this.states; j++)
			{
				this.A[i][j]=new Matrix(d_x,d_y);
				this.A[this.N-i-1][j]=new Matrix(d_y,d_x);
			}
		}
		this.FillARandom();
		for (int i=0; i<this.A[N-1][0].GetSize()[1]; i++)
		{
			this.R[N-1][i] = new Matrix(1,1);
			this.R[N-1][i].SetContent(0,0,1);
		}
		for (int i=0; i<this.A[0][0].GetSize()[1]; i++)
		{
			this.L[0][i] = new Matrix(1,1);
			this.L[0][i].SetContent(0,0,1);
		}

	//	this.sweep();
		
	}
	
	public Mps(int N, int d, int d_max, int states, Mpo OpH)
	{
		this.d=d;									// Maximale Groesse der Matrizen
		this.d_max=d_max;		
		this.N=N;									// Laenge der Kette
		this.states=states;								// Anzahl der moeglichen Zustaende (Anzahl der Matrizen pro Kettenplatz)
		this.A = new Matrix[this.N][this.states];
		this.R = new Matrix[this.N][this.d_max];
		this.L = new Matrix[this.N][this.d_max];
		for (int i=0; i<(this.N+1)/2; i++)
		{
			int d_x = (int)(Math.pow(d,i+1)<this.d_max?Math.pow(d,i+1):this.d_max);
			int d_y = (int)(Math.pow(d,i)<this.d_max?Math.pow(d,i):this.d_max);
//			System.out.format("i: %d, y: %d, x: %d%n", i, d_y, d_x);
			for (int j=0; j<this.states; j++)
			{
				this.A[i][j]=new Matrix(d_x,d_y);
				this.A[this.N-i-1][j]=new Matrix(d_y,d_x);
			}
		}
		if (this.N%2==1)
			for (int j=0;j<this.states;j++)
				this.A[this.N/2][j]=new Matrix((int)(Math.pow(d,this.N/2)<d_max?Math.pow(d,this.N/2):d_max),(int)(Math.pow(d,this.N/2)<d_max?Math.pow(d,this.N/2):d_max));
		this.FillARandom();
		for (int i=0; i<this.A[N-1][0].GetSize()[1]; i++)
		{
			this.R[N-1][i] = new Matrix(1,1);
			this.R[N-1][i].SetContent(0,0,1);
		}
		for (int i=0; i<this.A[0][0].GetSize()[1]; i++)
		{
			this.L[0][i] = new Matrix(1,1);
			this.L[0][i].SetContent(0,0,1);
		}

		this.DoLeftSweepFull(OpH);
	}
	
	public Mps(int N, int d, int d_max, int states)
	{
		this.d=d;									// Maximale Groesse der Matrizen
		this.d_max=d_max;		
		this.N=N;									// Laenge der Kette
		this.states=states;								// Anzahl der moeglichen Zustaende (Anzahl der Matrizen pro Kettenplatz)
		this.A = new Matrix[this.N][this.states];
		this.R = new Matrix[this.N][this.d_max];
		this.L = new Matrix[this.N][this.d_max];
		for (int i=0; i<(this.N+1)/2; i++)
		{
			int d_x = (int)(Math.pow(d,i+1)<d_max?Math.pow(d,i+1):d_max);
			int d_y = (int)(Math.pow(d,i)<d_max?Math.pow(d,i):d_max);
			for (int j=0; j<this.states; j++)
			{
				this.A[i][j]=new Matrix(d_x,d_y);
				this.A[this.N-i-1][j]=new Matrix(d_y,d_x);
			}
		}
		if (this.N%2==1)
			for (int j=0;j<this.states;j++)
				this.A[this.N/2][j]=new Matrix((int)(Math.pow(d,this.N/2)<d_max?Math.pow(d,this.N/2):d_max),(int)(Math.pow(d,this.N/2)<d_max?Math.pow(d,this.N/2):d_max));
		
		for (int i=0; i<this.A[N-1][0].GetSize()[1]; i++)
		{
			this.R[N-1][i] = new Matrix(1,1);
			this.R[N-1][i].SetContent(0,0,1);
		}
		for (int i=0; i<this.A[0][0].GetSize()[1]; i++)
		{
			this.L[0][i] = new Matrix(1,1);
			this.L[0][i].SetContent(0,0,1);
		}
	}
	
	private Matrix ReformMatrix(Matrix c)
	{
		Matrix work = new Matrix((int)(c.GetSize()[1]/(2)), (int)(c.GetSize()[0]*(2)));
		for (int i=0; i<work.GetSize()[1]; i++)
			for (int s=0; s<work.GetSize()[0]/2; s++)
				work.SetContent(i,s,c.GetContent(i, s));
		for (int i=0; i<work.GetSize()[1]; i++)
			for (int s=0; s<work.GetSize()[0]/2; s++)
				work.SetContent(i,s+work.GetSize()[0]/2,c.GetContent(i+work.GetSize()[1], s));			
					work.Merges = new int[this.states-1];
		work.MergesUnder = new boolean[this.states-1];
		for (int k=0; k<this.states-1; k++)
		{
			work.Merges[k] = work.GetSize()[0]/2*(k+1);
			work.MergesUnder[k] = true;
		}
		return work;
	}
	
	public Mps(double[] c_in, int states, Mpo op)
	{
		this.states=states;								// Anzahl der moeglichen Zustaende (Anzahl der Matrizen pro Kettenplatz)
		this.N = (int)(Math.log(c_in.length)/Math.log(this.states));

		this.d=2;
		this.d_max=1024;
		this.A = new Matrix[this.N][this.states];
		this.R = new Matrix[this.N][this.d_max];
		this.L = new Matrix[this.N][this.d_max];

		Matrix c = new Matrix(c_in.length, 1);
		for (int i=0; i<c_in.length; i++)
			c.SetContent(i,0,c_in[i]);
 
		for (int j=1; j<this.N; j++)
		{
			Matrix u = new Matrix(), vt = new Matrix();
			Matrix work = ReformMatrix(c);
			Matrix s = work.svd(u, vt);

			c = s.mult(vt);
			for (int k=this.states-1; k>=0; k--)
				this.A[j-1][k] = u.split();

			c.Merges = vt.Merges;
			c.MergesUnder = vt.MergesUnder;
		}
				
		Matrix work = ReformMatrix(c);
		
	
		for (int k=this.states-1; k>=0; k--)
			this.A[N-1][k] = work.split();


		for (int i=0; i<this.A[N-1][0].GetSize()[1]; i++)
		{
			this.R[N-1][i] = new Matrix(1,1);
			this.R[N-1][i].SetContent(0,0,1);
		}
		for (int i=0; i<this.A[N-1][0].GetSize()[1]; i++)
		{
			this.L[N-1][i] = new Matrix(1,1);
			this.L[N-1][i].SetContent(0,0,1);
		}		
		this.DoLeftSweepFull(op);		
	}

	public void FillARandom()
	{
		for (int i=0; i<this.N; i++)
			for (int j=0; j<this.states; j++)
				for (int x=0; x<this.A[i][j].GetSize()[1]; x++)
					for (int y=0; y<this.A[i][j].GetSize()[0]; y++)
						this.A[i][j].SetContent(x,y,(int)(Math.random()*10+1));
	}
	
	public String toString()
	{
		String newString = "";
		for (int i=0; i<this.N; i++)
			for (int j=0; j<this.states; j++)
				newString = newString + this.A[i][j];
		return newString;
	}
	
	public void printreadable()
	{
			for (int i=0; i<this.N; i++)
		{
			for (int k=0; k<this.states; k++)
			{
				System.out.println("A[" + i + "][" + k + "]");
				this.A[i][k].printreadable();
			}
		}	
	}

	public void printLreadable()
	{
		for (int i=0; i<this.N-1; i++)
		{
			for (int k=0; k<this.A[i][0].GetSize()[1]; k++)
			{
				System.out.println("L[" + i + "][" + k + "]");
				this.L[i][k].printreadable();
			}
		}	
	}
	
	public void printRreadable()
	{
		for (int i=0; i<this.N; i++)
		{
			for (int k=0; k<this.A[i][0].GetSize()[1]; k++)
			{
				System.out.println("R[" + i + "][" + k + "]");
				this.R[i][k].printreadable();
			}
		}	
	}
	
	public void DoRightSweepFull(Mpo op)
	{		
		for (int j=0; j<this.N-1; j++)
			this.DoRightSweepStep(j, op);
	}
	
	public void DoLeftSweepFull(Mpo op)
	{		
		for (int j=this.N-1; j>0; j--)
			this.DoLeftSweepStep(j, op);
	}
	
	public Matrix DoRightSweepStep(int j, Mpo op)
	{
		Matrix AAllStates = new Matrix(this.A[j][0]);
		for (int k=1; k<this.states; k++)
			AAllStates.merge(this.A[j][k], true);

		Matrix u=new Matrix();
		Matrix vt=new Matrix();
		Matrix s=AAllStates.svd(u, vt);
		
		for (int i=u.Merges.length;i>=0;i--)
			this.A[j][i]=u.split();
		for (int i=0; i<this.states; i++)
			this.A[j+1][i]=(s.mult(vt)).mult(this.A[j+1][i]);
		this.BuildL(op, j);
		return s;
	}
	
	public Matrix DoLeftSweepStep(int j, Mpo op)
	{
		Matrix AAllStates = new Matrix(this.A[j][0]);
		for (int k=1; k<this.states; k++)
			AAllStates.merge(this.A[j][k], false);
		Matrix u=new Matrix();
		Matrix vt=new Matrix();
		Matrix s=AAllStates.svd(u, vt);
		for (int i=vt.Merges.length;i>=0;i--)
			this.A[j][i]=vt.split();
		for (int i=0; i<this.states; i++)
			this.A[j-1][i]=(this.A[j-1][i].mult(u).mult(s));
		this.BuildR(op, j-1);
		return s;
	}

	public void testOrtho()
	{
		Matrix test = new Matrix();
		String out = "";
		for (int i=0; i<this.N; i++)
		{
			boolean a = true;
			boolean b = true;

			test = this.A[i][0].trans().mult(this.A[i][0]);

			for (int j=1; j<this.states; j++)
				test = test.add(this.A[i][j].trans().mult(this.A[i][j]));
			for (int k=0; k<test.GetSize()[0]; k++)
				for (int l=0; l<test.GetSize()[1]; l++)
					if ((Math.round(test.GetContent(k,l))!=0 && k!=l) || (Math.round(test.GetContent(k,l))!=1 && k==l))
						a = false;

			test = this.A[i][0].mult(this.A[i][0].trans());
			for (int j=1; j<this.states; j++)
				test = test.add(this.A[i][j].mult(this.A[i][j].trans()));
			for (int k=0; k<test.GetSize()[0]; k++)
				for (int l=0; l<test.GetSize()[1]; l++)
					if ((Math.round(test.GetContent(k,l))!=0 && k!=l) || (Math.round(test.GetContent(k,l))!=1 && k==l))
						b = false;

			if (a && b)
				out = out + "X";
			else if (a)
				out = out + "A";
			else if (b)
				out = out + "B";
			else
				out = out + "M";
		}
		System.out.println(out);
	}
	
	public double ScalarProd(Mps b)
	{
		if (this.N != b.N)
			return 0.0;
		Matrix out = new Matrix();
		Matrix pre_out = new Matrix();
		pre_out = this.A[0][0].trans().mult(b.A[0][0]);
		for (int s=1; s<this.states;s++)
			pre_out = pre_out.add(this.A[0][s].trans().mult(b.A[0][s]));
		out = new Matrix(pre_out);
		for (int i=1;i<this.N;i++)
		{
			pre_out = this.A[i][0].trans().mult(out).mult(b.A[i][0]);
			for ( int s=1; s<this.states;s++)
				pre_out = pre_out.add(this.A[i][s].trans().mult(out).mult(b.A[i][s]));
			out = new Matrix(pre_out);
		}
		return out.GetContent(0,0);
	}
	
	public void BuildR(Mpo H)
	{
		System.out.println(this.N);
		for (int i=this.N-2; i>=0; i--)
			this.BuildR(H, i);
	}
	
	public void BuildR(Mpo H, int N)
	{
		Matrix Eins = new Matrix(1,1); Eins.SetContent(0,0,1);
		for (int a=0; a<this.A[N+1][0].GetSize()[0];a++)
		{
			this.R[N][a] = new Matrix(H.W[N+1][0][0].GetSize()[0], this.A[N+1][0].GetSize()[0]);
			for (int j=0; j<this.states; j++)
				for (int k=0; k<this.states; k++)
					for (int l=0; l<this.A[N+1][0].GetSize()[1]; l++)
						if (N>this.N-2)
							for (int m=0; m<this.A[N+1][j].GetSize()[0]; m++)
								for (int o=0; o<H.W[N+1][j][k].GetSize()[0]; o++)
									this.R[N][a].SetContent(o,m, this.R[N][a].GetContent(o,m)+this.A[N+1][j].GetContent(0,m)*H.W[N+1][j][k].GetContent(0,o)*this.A[N+1][k].GetContent(0,a));
						else
							this.R[N][a] = new Matrix(this.R[N][a].add(H.W[N+1][j][k].mult(this.R[N+1][l].trans()).mult(this.A[N+1][k].trans()).mult(this.A[N+1][j].GetContent(l,a)).trans()));
		}
	}
	
	public void BuildL(Mpo H)
	{
		for (int i=0; i<this.N; i++)
			this.BuildL(H, i);
	}
	
	public void BuildL(Mpo H, int N)
	{
		for (int a=0; a<this.A[N][0].GetSize()[1];a++)
		{
			this.L[N][a] = new Matrix(H.W[N][0][0].GetSize()[1], this.A[N][0].GetSize()[1]);
			for (int j=0; j<this.states; j++)
				for (int k=0; k<this.states; k++)
					for (int l=0; l<this.A[N][0].GetSize()[0]; l++)
						if (N==0)
							this.L[N][a] = new Matrix(this.L[N][a].add(H.W[N][j][k].trans().mult(this.A[N][k]).mult(this.A[N][j].GetContent(a,0)).trans()));
						else
							this.L[N][a] = new Matrix(this.L[N][a].add((this.L[N-1][l].mult(H.W[N][j][k]).trans()).mult(this.A[N][k]).mult(this.A[N][j].GetContent(a,l)).trans()));
		}
	}
	
	public void waitKB() 
	{
		try
		{
			while( System.in.read() != (int)'c' ) 
			{
				Thread.sleep(100); //Let the CPU do more worthwhile things for 0.1 seconds ;)
			}
		}
		catch (java.io.IOException e)
		{
			System.out.println("Fehler beim Schreiben der Datei");
		}
		catch (java.lang.InterruptedException e)
		{
			System.out.println("Fehler beim Schreiben der Datei");
		}
	}
	
	public void printState()
	{
		for (int i=0;i<Math.pow(this.states, this.N);i++)
		{
			Matrix temp = new Matrix(this.A[0][i%this.states]);
			System.out.print("State " + i%this.states + " ");
			for (int n=1; n<this.N; n++)
			{
				System.out.print((i/((int)(Math.pow(this.states,n)))%this.states) + " ");
				temp = new Matrix(temp.mult(this.A[n][(i/((int)(Math.pow(this.states,n)))%this.states)]));
			}
			System.out.println("--> " + temp.GetContent(0, 0));
		}
	}
	
}
