from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import os
import collections



def main():


	#scales=[4,6,8,16,24,32]
	#numofseeds=np.array([10,10,10,48,100,200])
	#startseed=1

	scales=[2,4,8,12,16,20,26,32]
	numofseeds=np.array([1,2,12,16,20,25,30,50])

	startseed=1
	dim=3

	numofseeds=numofseeds+startseed

	mapa=np.loadtxt(('vecconec.txt')).astype(int)

	if dim==2:
		LL=int(np.sqrt(mapa.shape[0]))
		mapa=mapa.reshape(LL,LL)

	if dim==3:
		LL=int(np.cbrt(mapa.shape[0]))
		mapa=mapa.reshape(LL,LL,LL)
	res, names=doforsubS_computeCmap(mapa,scales,postConec, compCon,dim,[],numofseeds)

	with open('keysCon.txt', 'w') as f:
		for item in names:
			f.write("%s\n" % item)

	np.save('ConResScales.npy',res)


	return

def doforsubS_computeCmap(mapa,scales,funpost, funcompCmap,dim,args,numofseeds):

	L=mapa.shape[0]
	res=dict()
	names=[]


	with open('Kfield.don') as f:
		seed = int(f.readline())

	for iscale in range(len(scales)):

		l=scales[iscale]

		if numofseeds[iscale] > seed:   #guarda aca

			nblocks=L//l #for each dimension

			if dim==2:
				for i in range(nblocks):
					for j in range(nblocks):
						cmapa=funcompCmap(mapa[i*l:(i+1)*l,j*l:(j+1)*l],dim)
						dats,names=funpost(cmapa,dim,args)
						if i== 0 and j==0:
							for icon in range(len(names)):
								res[l,names[icon]]=[]
						for icon in range(len(names)):
							res[l,names[icon]]+=[dats[icon]]	

	
			if dim==3:
				for i in range(nblocks):
					for j in range(nblocks):
						for k in range(nblocks):
							cmapa=funcompCmap(mapa[i*l:(i+1)*l,j*l:(j+1)*l,k*l:(k+1)*l],dim)
							dats, names=funpost(cmapa,dim,args)
							if i== 0 and j==0 and k==0:
								for icon in range(len(names)):
									res[l,names[icon]]=[]
							for icon in range(len(names)):
								res[l,names[icon]]+=[dats[icon]]						



	return res, names


def ConConfig(L,dim):

	params=[]
	if dim==2:
		params=['1','4','imap.txt',str(L)+' '+str(L),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF']
		execCon='conec2d'

	if dim==3:
		params=['1','6','imap.txt',str(L)+' '+str(L)+' ' +str(L),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF']
		execCon='conec3d'
	return params, execCon



def compCon(mapa,dim):

	exeDir='./'
	L=mapa.shape[0]
	params,execCon=ConConfig(L,dim)

	with open(exeDir+'coninput.txt', 'w') as f:
		for item in params:
			f.write("%s\n" % item)
	np.savetxt(exeDir+params[2],mapa.reshape(-1))

	#wiam=os.getcwd()
	#os.chdir(exeDir)
	os.system('cp ../../../tools/conec3d ./')
	os.system(' ./'+execCon +'>/dev/null') #'cd ' +exeDir+
	
	cmapa=np.loadtxt(params[-2]).reshape(mapa.shape).astype(int) #exeDir+
	#os.chdir(wiam)
	return cmapa





def postConec(cmap,dim,args): 

	names=['PPHA','VOLALE','ZNCC','zintcc','spaninning','npz','npy','npx',]


	L=cmap.shape[0]
	results=[]
	names=[]


	y = np.bincount(cmap.reshape(-1))
	ii = np.nonzero(y)[0]
	cf=np.vstack((ii,y[ii])).T #numero de cluster, frecuencia

	if cf[0,0]==0:
		cf=cf[1:,:] #me quedo solo con la distr de tamanos, elimino info cluster cero

	if cf.shape[0]>0:
	#	headers=['N','p','Csize','CLenX','CLenY','CmaxVol','MaxLenX','MaxLenY','NpcX','NpcY']

		nper=np.sum(cf[:,1]) #num de celdas permeables
		nclus=cf.shape[0] #cantidad de clusters


		#ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA
		results+=[nper/np.size(cmap)]  #ppha
		results+=[np.max(cf[:,1])/nper] #volale #corregido va entre [0,p]
		results+=[nclus]   #zncc
		results+=[np.sum(cf[:,1]**2)/np.size(cmap)/nper]  #gamma, recordar zintcc =gamma*p
	
		spanning, pclusZ, pclusY, pclusX =get_perco(cmap,dim)
		results+=[spanning, len(pclusZ), len(pclusY), len(pclusX)]
	

		results+=Plen(spanning,cmap,cf,dim)


	names+=['PPHA']
	names+=['VOLALE']
	names+=['ZNCC']
	names+=['ZINTCC']
	names+=['spanning', 'npz', 'npy', 'npx']
	names+=['Plen','S','P']

	if cf.shape[0]==0:
		for i in range(len(names)):
			results+=[0]
	return  results, names


#ZINTCC,VOLALE,ZGAMMA,ZIPZ,ZNCC,PPHA


def get_pos2D(cmap,cdis):
	
	Ns=cdis.shape[0]
	pos=dict()
	i=0
	for cnum in cdis[:,0]:
		pos[cnum]=np.zeros((cdis[i,1]+1,2)) #+1 porque uso de flag
		i+=1

	for i in range(cmap.shape[0]):
		for j in range(cmap.shape[1]):
			if cmap[i,j] != 0:
				flag=int(pos[cmap[i,j]][0,0])+1
				pos[cmap[i,j]][0,0]=flag
				pos[cmap[i,j]][flag,0]=i
				pos[cmap[i,j]][flag,1]=j



	return pos


def get_pos3D(cmap,cdis):
	
	Ns=cdis.shape[0]
	pos=dict()
	i=0
	for cnum in cdis[:,0]:
		pos[cnum]=np.zeros((cdis[i,1]+1,3))
		i+=1
	for i in range(cmap.shape[0]):
		for j in range(cmap.shape[1]):
			for k in range(cmap.shape[2]):

				if cmap[i,j,k] != 0:
					flag=int(pos[cmap[i,j,k]][0,0])+1
					pos[cmap[i,j,k]][0,0]=flag
					pos[cmap[i,j,k]][flag,0]=i
					pos[cmap[i,j,k]][flag,1]=j
					pos[cmap[i,j,k]][flag,2]=k


	return pos

def Plen(spannng,cmap,cdis,dim):

	if dim==2:
		return P_len2D(spannng,cmap,cdis)
	if dim==3:
		return P_len3D(spannng,cmap,cdis)
	return []

def P_len2D(spanning,cmap,cdis):


	pos = get_pos2D(cmap,cdis)
	#print(summary['NpcY'],summary['NpcX'],summary['PPHA'])
	
	den=0
	num=0

	nperm=np.sum(cdis[:,1])
	if spanning > 0:
		amax=np.argmax(cdis[:,1])
		P=cdis[amax,1]/nperm
		cdis=np.delete(cdis,amax,axis=0)

	else:
		P=0

	i=0
	if cdis.shape[0]> 0:
		S=np.sum(cdis[:,1])/(cdis.shape[0])
		for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k-
			mposx, mposy = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]) #el 1: de sacar el flag
			Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2)  #Rs cuadrado  ecuacion 12.9 libro Harvey Gould, Jan Tobochnik
			num += cdis[i,1]**2 * Rs
			den+=cdis[i,1]**2
			i+=1
		return [np.sqrt(num/den), S, P]
	else:
		return [0,0,P]




def P_len3D(spanning,cmap,cdis):


	pos = get_pos3D(cmap,cdis)
	#print(summary['NpcY'],summary['NpcX'],summary['PPHA'])
	
	den=0
	num=0

	nperm=np.sum(cdis[:,1])
	if spanning > 0:
		amax=np.argmax(cdis[:,1])
		P=cdis[amax,1]/nperm
		cdis=np.delete(cdis,amax,axis=0)

	else:
		P=0

	i=0
	if cdis.shape[0]> 0:
		S=np.sum(cdis[:,1])/(cdis.shape[0])
		for cnum in cdis[:,0]: #los clusters estan numerados a partir de 1, cluster cero es k-
			mposx, mposy, mposz = np.mean(pos[cnum][1:,0]), np.mean(pos[cnum][1:,1]), np.mean(pos[cnum][1:,2]) #el 1: de sacar el flag
			Rs =np.mean((pos[cnum][1:,0]-mposx)**2 +(pos[cnum][1:,1]-mposy)**2+(pos[cnum][1:,2]-mposz)**2)  #Rs cuadrado  ecuacion 12.9 libro Harvey Gould, Jan Tobochnik
			num += cdis[i,1]**2 * Rs
			den+=cdis[i,1]**2
			i+=1
		return [np.sqrt(num/den), S, P]
	else:
		return [0,0,P]




def get_perco(cmap,dim):

	if dim==2:

		pclusY=[] #list of the percolating clusters
		for i in range(cmap.shape[0]):
			if cmap[i,0] != 0:
				if cmap[i,0] not in pclusY:
					if cmap[i,0] in cmap[:,-1]:
						pclusY+=[cmap[i,0]]


		pclusZ=[] #list of the percolating clusters Z direction, this one is the main flow in Ndar.py, the fixed dimension is the direction used to see if pecolates
		for i in range(cmap.shape[1]):
			if cmap[0,i] != 0:
				if cmap[0,i] not in pclusZ:
					if cmap[0,i] in cmap[-1,:]:  #viendo sin en la primer cara esta el mismo cluster que en la ultima
						pclusZ+=[cmap[0,i]]


		pclusX=[]
		spanning=0
		if len(pclusZ)==1 and pclusZ==pclusY:
			spanning=1
		

	if dim==3:

		
		pclusX=[] #list of the percolating clusters
		for i in range(cmap.shape[0]): #  Z
			for j in range(cmap.shape[1]): #X
				if cmap[i,j,0] != 0:
					if cmap[i,j,0] not in pclusX:
						if cmap[i,j,0] in cmap[:,:,-1]:
							pclusX+=[cmap[i,j,0]]

		pclusY=[] #list of the percolating clusters
		for i in range(cmap.shape[0]): #  Z
			for k in range(cmap.shape[2]): #X
				if cmap[i,0,k] != 0:
					if cmap[i,0,k] not in pclusY:
						if cmap[i,0,k] in cmap[:,-1,:]:
							pclusY+=[cmap[i,0,k]]

		pclusZ=[] #list of the percolating clusters
		for k in range(cmap.shape[2]): #x  
			for j in range(cmap.shape[1]): #y
				if cmap[0,j,k] != 0:
					if cmap[0,j,k] not in pclusZ:
						if cmap[0,j,k] in cmap[-1,:,:]:
							pclusZ+=[cmap[0,j,k]]  #this is the one

		spanning=0
		if len(pclusZ)==1 and pclusZ==pclusY and pclusZ==pclusX:
			spanning=1


	return spanning, pclusZ, pclusY, pclusX


main()