import numpy as np
import scipy
from scipy.optimize import minimize
import random
import math
#This code has a couple of loops to decompose T with SVD, and then perform an approximation by adding
#the biggest Dc singular values terms. approxcheck1/approxcheck2 gives the substraction from the exact tensor and
#the approximation. check1/check2 just makes sure the SVD is working appropriately.
D=2             #Bond dimensionality
Dc=15           #Defines the cutoff Dc for the summation of singular values. 
d=7             #Physical index dimension
T=np.zeros((D,D,D,D,d),dtype=complex)        #Original tensor 
TT=np.zeros((D,D,D,D,D,D,D,D), dtype=complex)  #\mathbb{T} in the paper
Tstar=np.conjugate(T)


def main():
    for a in range(D):
        for b in range(D):
            for c in range(D):
                for d in range(D):
                    for e in range(D):
                        for f in range(D):
                            for g in range(D):
                                for h in range(D):
                                    for m in range(D):
                                        T[e][f][g][h][m] =random.uniform(0.4,0.43)+0*np.complex(0,1)
                                        Tstar[e][f][g][h][m]=np.conjugate(T[e][f][g][h][m])
	                                TT[a][b][c][d][e][f][g][h] += Tstar[a][b][c][d][m]*T[e][f][g][h][m]
        

    M=TT.reshape((16,16))

    U, s, V=np.linalg.svd(M)
    X1=np.zeros((D**4,D**4),dtype=complex)
    X2=np.zeros((D**4,D**4),dtype=complex)
    for a in range(D**4):
        for b in range(D**4):
                X1[a][b]=U[a][b]*math.sqrt(s[b])
                X2[a][b]=V[a][b]*math.sqrt(s[a])
            
    S1=X1.reshape((D**2,D**2,D**4))
    S2=X2.reshape((D**4,D**2,D**2))

    reshapedTT=np.zeros((D**2,D**2,D**2,D**2),dtype=complex)
    for a in range(D**2):
        for b in range(D**2):
            for c in range(D**2):
                for d in range(D**2):
                    for n in range(D**4):
                        reshapedTT[a][b][c][d]+=S1[a][b][n]*S2[n][c][d]
                    
    check1 = np.round(reshapedTT.reshape((D,D,D,D,D,D,D,D))-TT,13)      #This two are to check that our decomposition is indeed working. (Mainly the reshape functions)
    check2 = np.round(reshapedTT-TT.reshape(D**2,D**2,D**2,D**2),13)

    approxreshapedTT=np.zeros((D**2,D**2,D**2,D**2),dtype=complex)
    for a in range(D**2):
        for b in range(D**2):
            for c in range(D**2):
                for d in range(D**2):
                    for n in range(Dc):
                        approxreshapedTT[a][b][c][d]+=S1[a][b][n]*S2[n][c][d]

    approxcheck1 = np.round(approxreshapedTT.reshape((D,D,D,D,D,D,D,D))-TT,13)      #This two are to check that our approximation converges.
    approxcheck2 = np.round(approxreshapedTT-TT.reshape(D**2,D**2,D**2,D**2),13)

    sum=0
    for i in range(D**2):
        for j in range(D**2):
            for k in range(D**2):
                for l in range(D**2):
                    sum+=abs(approxcheck2[i,j,k,l]/TT.reshape(D**2,D**2,D**2,D**2)[i,j,k,l])
    return sum/(D**8)
a=0 
for i in range(100):   
    a+=main()
a=a/100.