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()