package dmrgvisu2;

import javax.swing.JTextArea;

public class Dmrg
{
	private static final long serialVersionUID = 1L;
	Mps state;
	Mpo operator;
	
	int CurrentSite;
	boolean Direction;	// True == left; False == right
	double CurrentEnergy;
	double MinEnergy;
	double[] CurrentDm;
	double CurrentErw;
	double CurrentEe;
	int CurrentSweep;
	JTextArea EnergyText;
	JTextArea SweepText;
	JTextArea SiteText;
	double[] CurrentSxExV;
	double[] CurrentSzExV;
	
	Mpo Sz;
	Mpo Sx;
	Matrix s_z;
	Matrix s_x;
	double sx_full;
	double sz_full;
	double sx_full_without_abs;
	double sz_full_without_abs;	
	int UseVisuRef;
	
	visu visuRef;
		
	public Dmrg(Mps state, Mpo operator, visu visuRef)
	{
		this.state = state;
		this.operator = operator;
		this.CurrentSite = 0;
		this.Direction = false;
		this.CurrentEnergy = 0;
		this.MinEnergy = 0;
		this.CurrentDm = new double[0];
		this.CurrentErw = 0;
		this.CurrentEe = 0;
		this.CurrentSweep = 0;
		this.CurrentSxExV = new double[this.state.N];
		this.CurrentSzExV = new double[this.state.N];
		for (int i=0; i<this.state.N; i++)
		{
			this.CurrentSzExV[i]=0;
			this.CurrentSxExV[i]=0;
		}
		this.visuRef = visuRef;
		this.UseVisuRef = 1;
	}
	
	public Dmrg(Mps state, Mpo operator, Mpo Sz, Mpo Sx, Matrix s_z, Matrix s_x)
	{
		this.state = state;
		this.operator = operator;
		this.CurrentSite = 0;
		this.Direction = false;
		this.CurrentEnergy = 0;
		this.MinEnergy = 0;
		this.CurrentDm = new double[0];
		this.CurrentErw = 0;
		this.CurrentEe = 0;
		this.CurrentSweep = 0;
		this.CurrentSxExV = new double[this.state.N];
		this.CurrentSzExV = new double[this.state.N];
		for (int i=0; i<this.state.N; i++)
		{
			this.CurrentSzExV[i]=0;
			this.CurrentSxExV[i]=0;
		}
		this.visuRef = visuRef;
		this.Sz=Sz;
		this.Sx=Sx;
		this.s_z = s_z;
		this.s_x = s_x;
		this.UseVisuRef = 0;
	}
	
	public int DoNextStep()
	{
		if (this.CurrentSite == this.state.N-1)
			this.Direction = true;
		if (this.CurrentSite == 0)
		{
			this.Direction = false;
			this.CurrentSweep++;
		}
		if (this.Direction == false)
		{
			this.DoRightStep(this.CurrentSite);
			this.CurrentSite++;
			return this.CurrentSite-1;
		}
		else
		{
			this.DoLeftStep(this.CurrentSite);
			this.CurrentSite--;
			return this.CurrentSite+1;
		}
		
	}
	
	public void DoRightSweep()
	{
		for (int j=0; j<this.state.N-1; j++)
			this.DoRightStep(j);
	}
	
	public void DoLeftSweep()
	{
		for (int j=this.state.N-1; j>0; j--)
			this.DoLeftStep(j);
	}
	
	public void DoRightStep(int N)
	{
		Matrix H = this.BuildHeff(N);
		double[] EigenValue = H.Eigen();
		Matrix[] M = new Matrix[this.state.states];
		for (int i=0; i<this.state.states; i++)
		{
			M[i] = new Matrix(this.state.A[N][0].GetSize()[1], this.state.A[N][0].GetSize()[0]);
			for (int m=0; m<this.state.A[N][0].GetSize()[0]; m++)
				for (int o=0; o<this.state.A[N][0].GetSize()[1]; o++)
					M[i].SetContent(o, m, H.GetContent(0, (i*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1])+(m*this.state.A[N][0].GetSize()[1])+o));
			this.state.A[N][i] = new Matrix(M[i]);
		}
		this.CalcSxSz();
		Matrix s = this.state.DoRightSweepStep(N, this.operator);
		s = s.mult(s);
		this.CurrentDm = new double[s.GetSize()[0]];
		for (int i = 0; i<this.CurrentDm.length; i++)
			this.CurrentDm[i] = s.GetContent(i,i);
		this.CurrentEnergy = EigenValue[0];
	}
	
	public void CalcSxSz()
	{
		this.sx_full = 0;
		this.sz_full = 0;
		this.sx_full_without_abs = 0;
		this.sz_full_without_abs = 0;
		if (this.UseVisuRef == 1)
		{
			for (int i=0; i<this.state.N; i++)
			{
				this.visuRef.mpo_sz.SetLocal(i, this.visuRef.s_z);
				this.visuRef.mpo_sx.SetLocal(i, this.visuRef.s_x);
				this.CurrentSzExV[i] = this.state.ScalarProd(this.visuRef.mpo_sz.mult(this.state));
				sz_full += this.CurrentSzExV[i];
				this.CurrentSxExV[i] = this.state.ScalarProd(this.visuRef.mpo_sx.mult(this.state));
				this.sx_full += this.CurrentSxExV[i];
				this.visuRef.mpo_sz.DelLocal(i);
				this.visuRef.mpo_sx.DelLocal(i);
			}
		}
		else
		{			
			for (int i=0; i<this.state.N; i++)
			{
				this.Sz.SetLocal(i, this.s_z);
				this.Sx.SetLocal(i, this.s_x);
				this.CurrentSzExV[i] = this.state.ScalarProd(this.Sz.mult(this.state));
				this.sz_full += Math.abs(this.CurrentSzExV[i]);
				this.sz_full_without_abs += this.CurrentSzExV[i];
				this.CurrentSxExV[i] = this.state.ScalarProd(this.Sx.mult(this.state));
				this.sx_full += Math.abs(this.CurrentSxExV[i]);
				this.sx_full_without_abs += this.CurrentSxExV[i];
				this.Sz.DelLocal(i);
				this.Sx.DelLocal(i);
			}
		}
		//System.out.println("Sx: " + this.sx_full + "    Sz: " + sz_full);
		return;
	}
	public void DoLeftStep(int N)
	{
		Matrix H = this.BuildHeff(N);
		double[] EigenValue = H.Eigen();
		Matrix[] M = new Matrix[this.state.states];
		for (int i=0; i<this.state.states; i++)
		{
			M[i] = new Matrix(this.state.A[N][0].GetSize()[1], this.state.A[N][0].GetSize()[0]);
			for (int n=0; n<this.state.A[N][0].GetSize()[1]; n++)
				for (int p=0; p<this.state.A[N][0].GetSize()[0]; p++)
					M[i].SetContent(n, p, H.GetContent(0, (i*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1])+(p*this.state.A[N][0].GetSize()[1])+n));
			this.state.A[N][i] = new Matrix(M[i]);
		}
		this.CalcSxSz();

		Matrix s = this.state.DoLeftSweepStep(N, this.operator);
		s = s.mult(s);
		this.CurrentDm = new double[s.GetSize()[0]];
		for (int i = 0; i<this.CurrentDm.length; i++)
			this.CurrentDm[i] = s.GetContent(i,i);
		this.CurrentEnergy = EigenValue[0];
	}

	public double ScalarProd()
	{
		return this.state.ScalarProd(this.operator.mult(state));
	}
	
	public double ErwW(int N)
	{
		double out = 0.0;
		Matrix H = this.BuildHeff(N);
		
		for (int i=0; i<this.state.states; i++)
			for (int j=0; j<this.state.states; j++)
				for (int m=0; m<this.state.A[N][0].GetSize()[0]; m++)
					for (int n=0; n<this.state.A[N][0].GetSize()[0]; n++)
						for (int o=0; o<this.state.A[N][0].GetSize()[1]; o++)
							for (int p=0; p<this.state.A[N][0].GetSize()[1]; p++)
								out = out+H.GetContent((j*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1])+(n*this.state.A[N][0].GetSize()[1])+p, (i*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1])+(m*this.state.A[N][0].GetSize()[1])+o)*this.state.A[N][i].GetContent(o,m)*this.state.A[N][j].GetContent(p,n);
		return out;
	}	
	
	public Matrix BuildHeff(int N)
	{
		int Size = this.state.states*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1];
		Matrix H = new Matrix(Size, Size);
		Matrix[][][][] LWR = new Matrix[this.state.states][this.state.states][this.state.A[N][0].GetSize()[0]][this.state.A[N][0].GetSize()[1]];
		for (int i=0; i<this.state.states; i++)
			for (int j=0; j<this.state.states; j++)
				for (int m=0; m<this.state.A[N][0].GetSize()[0]; m++)
					for (int o=0; o<this.state.A[N][0].GetSize()[1]; o++)
						if (N-1 < 0)
							LWR[i][j][m][o] = new Matrix((this.operator.W[N][i][j].mult(this.state.R[N][o].trans())));
						else
							LWR[i][j][m][o] = new Matrix(this.state.L[N-1][m].mult(this.operator.W[N][i][j].mult(this.state.R[N][o].trans())));

		for (int i=0; i<this.state.states; i++)
			for (int j=0; j<this.state.states; j++)
				for (int m=0; m<this.state.A[N][0].GetSize()[0]; m++)
					for (int n=0; n<this.state.A[N][0].GetSize()[0]; n++)
						for (int o=0; o<this.state.A[N][0].GetSize()[1]; o++)
							for (int p=0; p<this.state.A[N][0].GetSize()[1]; p++)
								H.SetContent((j*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1])+(n*this.state.A[N][0].GetSize()[1])+p, (i*this.state.A[N][0].GetSize()[0]*this.state.A[N][0].GetSize()[1])+(m*this.state.A[N][0].GetSize()[1])+o, LWR[i][j][m][o].GetContent(p,n));
		return H;
	}
}