You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
369 lines
8.0 KiB
Python
369 lines
8.0 KiB
Python
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()
|
|
|
|
|
|
|
|
|