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.
401 lines
11 KiB
Python
401 lines
11 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()
|