import numpy as np
import os
import time
from JoinCmaps import *
#k[x,y,z]


def div_veccon(kc,kh,nbl,rundir):

	t0=time.time()
	kc=np.where(kc==kh,1,0).astype(int)

	tcmaps=time.time()
	kc=get_smallCmap(kc,nbl,rundir)
	tcmaps=time.time()-tcmaps
	#if s_scale<kc.shape[0]:
	kc=join(kc,nbl)


	y = np.bincount(kc.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
	nclus=cf.shape[0] #cantidad de clusters
	nper=np.sum(cf[:,1]) #num de celdas permeables

	print(nbl,nclus,float(nper)/(kc.size),	time.time()-t0)

	return np.array([nbl,nclus,float(nper)/(kc.size),time.time()-t0,	tcmaps,tcmaps/(time.time()-t0)])



def get_smallCmap(vec,nbl,rundir):


	Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2]
	sx,sy,sz = Nx//nbl,Ny//nbl,Nz//nbl
	params, execCon = ConConfig(sx,sy,sz,Nz,rundir)
	if Nz==1:
		nblz=1
		sz=1
	else:
		nblz=nbl
	for i in range(nbl):
		for j in range(nbl):
			for k in range(nblz):
				vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz]=connec(vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],execCon,params,rundir)
	return vec


def connec(vec,execCon,params,rundir):
		np.savetxt(rundir+params[2],vec.reshape(-1))
		os.system(rundir+execCon +'>/dev/null') #'cd ' +exeDir++'>/dev/null'
		vec=np.loadtxt(params[-2]).reshape(vec.shape[0],vec.shape[1],vec.shape[2]).astype(int)
		return vec

def ConConfig(sx,sy,sz,Nz,rundir):

	params=[]
	if Nz==1:
		params=['1','4','vecconec.txt',str(sx)+' '+str(sy),'1.0 1.0','pardol.STA','pardol.CCO','pardol.COF']
		execCon='conec2d'

	else:
		params=['1','6','vecconec.txt',str(sx)+' '+str(sy)+' ' +str(sz),'1.0 1.0 1.0','30','pardol.STA','pardol.CCO','pardol.COF']
		execCon='conec3d'


	with open(rundir+'coninput.txt', 'w') as f:
		for item in params:
			f.write("%s\n" % item)

	return params, execCon

def join(vec,nbl):


	Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2]
	sx,sy,sz = Nx//nbl,Ny//nbl,Nz//nbl
	ex,ey,ez=np.log2(Nx),np.log2(Ny),np.log2(Nz)

	if Nz==1:
		sz=1
		nbz=1
		ez=1
		esz=1
	else:
		esz=np.log2(sz)
		
		
	esx,esy=np.log2(sx),np.log2(sy)

		

	for bs in range(0,int(ex-esx)):

		nbx,nby = int(2**(ex-esx-bs-1)),int(2**(ey-esy-bs-1))
		if Nz==1:
			sz=1
			nbz=1
		else:
			nbz=int(2**(ez-esz-bs-1))
			sz=Nz//nbz
		sx,sy=Nx//nbx,Ny//nby

		for i in range(nbx):
			for j in range(nby):
				for k in range(nbz):
					a=2
					vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz]=joinBox(vec[i*sx:(i+1)*sx,j*sy:(j+1)*sy,k*sz:(k+1)*sz],True,False)

	return vec

'''
job=0
k=np.load('../../data/'+str(job)+'/k.npy')
div_veccon(k,100,1,'./')	
div_veccon(k,100,2,'./')
div_veccon(k,100,4,'./')
'''

for job in range(6):
	k=np.load('../../data/'+str(job)+'/k.npy')
	print(job)
	res=div_veccon(k,100,4,'./')	
	np.savetxt('../../data/'+str(job)+'/Cmap_res.txt',res)
	res=div_veccon(k,100,1,'./')	
	#div_veccon(k,100,64,'./')
	#div_veccon(k,100,128,'./')