simulacion-permeabilidad/tools/connec/comp_cmap.py

import numpy as np
import os
import time
from tools.connec.JoinCmaps import *
import subprocess
from tools.connec.PostConec import ConnecInd
#k[x,y,z]
import json

def comp_connec(parser,rundir,nr):

	kc=np.load(rundir+'k.npy')
	keep_aspect = parser.get('Connectivity','keep_aspect')
	kh,sx = float(parser.get('Generation','kh')),int(parser.get('Connectivity','block_size'))
	S_min_post = int(parser.get('Connectivity','indicators_MinBlockSize'))
	nimax =2** int(parser.get('Connectivity','Max_sample_size'))

	gcon =bool(parser.get('Connectivity','compGconec'))

	if S_min_post ==-1 or S_min_post > kc.shape[0]:
		S_min_post=kc.shape[0]  #solo calcula indicadores para mayo escala
	if S_min_post ==0:
		S_min_post=sx  #solo calcula indicadores para escalas a partir del optimo
	if sx > S_min_post:
		sx = get_min_nbl(kc,nimax,nr,S_min_post)  #corta en mas artes para tener mediads de conec

	nbl=kc.shape[0]//sx


	if keep_aspect=='yes':
		keep_aspect=True
	else:
		keep_aspect=False

	t0=time.time()
	kc=np.where(kc==kh,1,0).astype(int)
	tcmaps=time.time()
	kc=get_smallCmap(kc,nbl,rundir,keep_aspect)
	tcmaps=time.time()-tcmaps

	kc,PostConTime=join(kc,nbl,keep_aspect,rundir,S_min_post,gcon)

	ttotal=time.time()-t0

	summary = np.array([nbl,ttotal,tcmaps/ttotal,PostConTime/ttotal])
	np.savetxt(rundir + 'ConnSummary.txt',summary,header='nbl,ttotal,tcmaps/ttotal,PostConTime/ttotal')
	np.save(rundir+'Cmap.npy',kc)
	return 


def get_min_nbl(kc,nimax,nr,smin):

	if kc.shape[2]==1:
		dim=2.0
	else:
		dim=3.0

	if nr>0:
		y=(1/dim)*np.log2(nr*kc.size/(nimax*(smin**dim)))
	else:	
		y=0
	y=int(y)
	s=int((2**y) * smin)
	if s<smin:
		s=smin

	return s
def get_smallCmap(vec,nbl,rundir,keep_aspect):


	Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2]
	sx = Nx//nbl
	if keep_aspect:
		sy,sz = Ny//nbl,Nz//nbl
		nblx, nbly,nblz = nbl, nbl, nbl
	else:
		sy,sz = sx,sx
		nblx=nbl
		nbly, nblz = Ny//sy, Nz//sz

	params, execCon = ConConfig(sx,sy,sz,Nz,rundir)
	if Nz==1:
		nblz=1
		sz=1
	os.system('cp ./tools/connec/'+execCon +' '+rundir)
	for i in range(nblx):
		for j in range(nbly):
			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)
	
	try: 
		temps=['pardol*','conec*d' ,'coninput.txt' ,'vecconec.txt']
		for temp in temps:
			os.system('rm '+rundir+temp)
	except:
		print('No connectivity temps to delete')
	
	return vec


def connec(vec,execCon,params,rundir):
		np.savetxt(rundir+params[2],vec.reshape(-1), fmt='%i')
		wd = os.getcwd()
		os.chdir(rundir)
		os.system('nohup ./'+execCon +' > connec.out 2>&1') #subprocess.call(['./tools/connec/'+execCon],cwd=rundir) #, '>/dev/null' , cwd=rundir
		os.chdir(wd)
		vec=np.loadtxt(rundir+params[-1]).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.CCO']
		execCon='conec2d'

	else:
		params=['1','6','vecconec.txt',str(sx)+' '+str(sy)+' ' +str(sz),'1.0 1.0 1.0','pardol.CCO']
		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,keep_aspect,datadir,S_min_post,gcon):

	Nx, Ny,Nz=vec.shape[0],vec.shape[1],vec.shape[2]
	sx = Nx//nbl

	if keep_aspect:
		sy,sz = Ny//nbl,Nz//nbl
		nblx, nbly,nblz = nbl, nbl, nbl
	else:
		sy,sz = sx,sx
		nblx=nbl
		nbly, nblz = Ny//sy, Nz//sz

	ex=np.log2(Nx)
	esx=np.log2(sx)
	join_z=True
	join_y=True
	if Nz==1:
		sz=1
		nblz=1


	post_time=0
	sxL=[sx]
	for bs in range(0,int(ex-esx)):


		if  vec.shape[0]==vec.shape[1] and sx>=S_min_post:
			t0=time.time()
			ConnecInd(vec,[sx],datadir)
			post_time=time.time()-t0
		sx,sy,sz = 2*sx,2*sy,2*sz
		sxL+=[sx]

		if sz > Nz:
			sz=Nz
			nblz=1
			join_z=False
		if sy > Ny:
			sy=Ny
			nbly=1
			join_y=False

		nblx,nbly,nblz = Nx//sx, Ny//sy, Nz//sz

		for i in range(nblx):
			for j in range(nbly):
				for k in range(nblz):
					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],join_y,join_z)

		if  vec.shape[0]==vec.shape[1] and sx>=S_min_post: #
			t0=time.time()
			ConnecInd(vec,[sx],datadir)
			post_time=post_time+(time.time()-t0)
	if gcon:
		ConnecInd(vec,sxL,datadir+'Global')
	return vec, post_time