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