import datetime, pathlib
from bisect import bisect_left, bisect_right
from math import hypot, fabs
from p_ggen import iterby2
from p_gmath import thanSegSeguw
from .dtmvar import ThanDTMDEM, _uniqint


class ThanDTMlines(ThanDTMDEM):
    "A set of lines which behaves as a Digital Terrain Model."

    def __init__(self, dxmax=None, dext=None):
        "Initialize DEM."
        if dxmax is None and dext is None:   #See manual_dtmlines.txt
            dxmax = dext = 20.0
        elif dxmax is None:                  #See manual_dtmlines.txt
            dxmax = dext/5.0
        elif dext is None:                   #See manual_dtmlines.txt
            dext = dxmax
        self.thanDxmax = dxmax #Max X distance of the end points of a line segment
        self.thanDext = dext   #X or Y distance that we extend from a point in order to find an intersection
        self.thanCena = (0.0, 0.0, 0.0)  #Centroid of the area of the dtm
        self.thanNori = 0      #Number of original line segments (just for information)
        self.thanLines = []


    def thanExportToPython(self, dtmname="dtm1", dir="."):
        "Create a python module which recreates this DTM."
        dir = pathlib.Path(dir)
        filnam = dir/(dtmname + ".py")
        try:
            fw = open(str(filnam), "w")
        except Exception as why:
            return 11, "Could not write Python code to file {}".format(filnam)
        fw.write("#DTMlines DTM autogenerated by ThanDTMlines.thanExportToPython()\n")
        fw.write("#{}\n".format(datetime.datetime.now().ctime()))
        fw.write("import p_gtri\n")

        fw.write("\ndef makedtm():\n")
        fw.write("    clines = (\n")
        for line1 in self.thanLines:
            assert len(line1) == 2, "All lines should have exactly 2 points: otherwise expthc() and impthc() are also wrong."
            fw.write("    ({:23.15e},{:23.15e},{:23.15e}),\n".format(*line1[0]))
            fw.write("    ({:23.15e},{:23.15e},{:23.15e}),\n".format(*line1[1]))
        fw.write("    )\n")
        fw.write("    d = p_gtri.ThanDTMlines(dxmax={:.15e}, dext={:.15e})\n".format(self.thanDxmax, self.thanDext))
        fw.write("    it = iter(clines)\n")
        fw.write("    for ca in it:\n")
        fw.write("        cb = next(it)\n")
        fw.write("        d.thanAddLine1((ca, cb))\n")
        fw.write("    d.thanRecreate()\n")
        fw.write("    return d\n")
        fw.write("\n{} = makedtm()\n".format(dtmname))
        fw.write("del makedtm\n")
        fw.close()
        return 0, ""


    def thanAddLine1(self, cp):
        "Add a (3d) line to the DEM."
        for ca, cb in iterby2(cp):
            ca = tuple(ca[:3])
            cb = tuple(cb[:3])
            self.thanNori += 1
            if ca[0] > cb[0]: ca, cb = cb, ca
            dx = cb[0]-ca[0]
            if dx <= self.thanDxmax:
                self.thanLines.append((ca, cb))
                continue
            n = int(dx/self.thanDxmax) + 1
            #ddx = self.thanDxmax/n   #Thanasis2022_06_11 commented out: 
            ddx = dx/n                #Thanasis2022_06_11
            #print("thanDxmax=", self.thanDxmax, "dx=", dx, "ddx=", ddx)
            x = 0.0
            cc = ca
            for i in range(n-1):
                x += ddx
                cd = [ca1+(cb1-ca1)/dx*x for (ca1,cb1) in zip(ca, cb)]
                self.thanLines.append((tuple(cc), tuple(cd)))
                cc = cd
            self.thanLines.append((tuple(cc), cb))

    def thanRecreate(self):
        "Recreate DEM after lines additions."
        if len(self.thanLines) < 1: return False, "DTM is empty!"
        self.thanLines.sort()
        self.thanCentroidCompute()
        return True, None

    def thanCentroidCompute(self):
        "Compute the centroid of all lines."
        cc = [0.0]*3
        sdis = 0.0
        xr = range(3)
        for ca, cb in self.thanLines:
            dis = hypot(cb[0]-ca[0], cb[1]-ca[1])
            for i in xr:
                cc[i] += (ca[i]+cb[i])*0.5 * dis
            sdis += dis
        self.thanCena = tuple(cc1/dis for cc1 in cc)

    def thanMinxy(self):
        "Find the minimun x and y coordinates (the DTM is assumed ordered)."
        cmin = list(self.thanLines[0][0])
        cmin[1] = min(c[1] for lin1 in self.thanLines for c in lin1)
        return cmin

    def thanIntersegZ(self, ca, cb, native=False):
        "Compute intersections of segment with DEM lines; don't sort intersections from ca to cb."
        ca = tuple(ca)
        cb = tuple(cb)
#        print "ca=", ca, "cb=", cb
        rev = False
        if ca[0] > cb[0]: ca, cb = cb, ca; rev=True
#        print "ca=", ca, "cb=", cb, "rev=", rev
        i = bisect_left(self.thanLines, ((ca[0]-self.thanDxmax,), (-1.0e30,)))
        j = bisect_right(self.thanLines, ((cb[0],),(+1.0e30,)))
        #print("i=", i, "j=", j, "len=", len(self.thanLines))
        #for ii in range(i, j):
        #    print("thanlines", ii,":", self.thanLines[ii])
        cint = []
        for k in range(i, j):
            c1, c2 = self.thanLines[k]
            uw = thanSegSeguw(ca, cb, c1, c2)
            if uw is None: continue
            cd = [ca1+(cb1-ca1)*uw[1] for (ca1,cb1) in zip(c1, c2)]
            if rev: cint.append((1-uw[0], cd))
            else:   cint.append((uw[0], cd))
        return cint


    def thanPointZ(self, cp, native=False):
        "Calculate the z coordinate of a point."
        for p in (1, 2, 4):
            cp = list(cp)     # bisect (called in __z1) needs cp as a list
            ca = list(cp)
            cb = list(cp)
            ca[0] -= self.thanDext*p           # |cb-ca| (1)
            cb[0] += self.thanDext*p
            d1, z1 = self.__z1(cp, ca, cb)     # Try to find intersections along x direction
#            if z1 is not None and fabs(z1) < 1.0: stop()
            ca = list(cp)
            cb = list(cp)
            ca[1] -= self.thanDext*p           # |cb-ca| must be the same as in (1), or the d1, d2..
            cb[1] += self.thanDext*p           # ..returned by __z1 will not be comparable
            d2, z2 = self.__z1(cp, ca, cb)     # Try to find intersections along y direction
#            if z2 is not None and fabs(z2) < 1.0: stop()

            #print("thanPointZ: d1, z1=", d1, z1)
            #print("thanPointZ: d2, z2=", d2, z2)
            if   z1 is None: z = z2            # Note that z2 may be None
            elif z2 is None: z = z1
            elif d1 < d2:    z = z1
            else:            z = z2
#            print "thanPointZ: z     =", z
            if z is not None: break
#        if z is not None and fabs(z) < 1.0: stop()
        return z


    def __z1(self, cp, ca, cb):
        "Calculate z with cp between ca and cb."
        cint = _uniqint(self.thanIntersegZ(ca, cb, native=True))
        #print("__z1   uniq: cint=", cint, "len=", len(cint))
        if len(cint) < 1: return None, None
        if len(cint) < 2:
            d1 = fabs(cint[0][0]-0.5)
            if d1 < 0.001: return d1, cint[0][1][2]  # 1 intersection, but exactly on the isoline, with 0.1% tolerance
            return None, None
        if cint[0][0] > 0.5:
            if cint[0][0] < 0.501: return cint[0][0]-0.5, cint[0][1][2]     # Allow for 0.1% error
            return None, None
        if cint[-1][0] < 0.5:
            if cint[-1][0] > 0.499: return 0.5-cint[-1][0], cint[-1][1][2]  # Allow for 0.1% error
            return None, None
        j = bisect_left(cint, (0.5, cp))
#        print "cint[j], j=", j
        i = j-1
        di = cint[i][0]
        zi = cint[i][1][2]
        dj = cint[j][0]
        zj = cint[j][1][2]
        if dj-di == 0.0:
#           The problem is as follows. Assume the very common case that two line segments
#           share an endnode. And now assume the rare case (but in certain applications
#           very common indeed) where we seek the z of exactly
#           this endnode. The program will find two intersections (one intersection for
#           each segment) that will have exactly the same coordinates, and the di and dj
#           (generalized distance) will be exactly 0.5 (for both intersections). This can
#           also happen if we have duplicate line segments.
#           Because we seek distance 0.5 in the bisect_left() function, the position
#           returned should be the first 0.5 (counting from left), so that if we have
#           only 2 points in cint, the i=j-1 will return -1 and we have the problem that
#           di==dj.
#           However bisect_left() does not work, because if distances are the same, then
#           it compares the lists (coordinates) that follow the distance. THUS PARALLEL
#           LISTS SHOULD BE USED. But then how do we do the sort?
#           Dirty solution follows:
#            for i,c in enumerate(cint): print i, ":", c
#            print "i=", i, "j=", j
            if dj == 0.5: return dj-0.5, zj
            for i in range(j-2, -1, -1):
                if cint[i][0] < dj: break
            else:
                assert 0, "Dirty solution unsuccessful; see comments above this statement!"
            di = cint[i][0]
            zi = cint[i][1][2]
        z1 = zi+(zj-zi)/(dj-di)*(0.5-di)
        d1 = min((0.5-di, dj-0.5))
#        print "z1=", z1, "d1=", d1
        return d1, z1


    def thanExpThc1(self, fw):
        "Saves the lines of the DTM to a .thcx file."
        f = fw.formFloat
        fw.writeAtt("maxdistancex", f % (self.thanDxmax,))
        fw.writeAtt("extensionx", f % (self.thanDext,))
        fw.writeSnode("centroid", 3, self.thanCena)
        fw.writeAtt("Original_segments", "%d" % (self.thanNori,))
        fw.writeNodes(c for lin in self.thanLines for c in lin)

    def thanImpThc1(self, fr, ver, than):
        "Reads the lines of the DTM from a .thc file."
        self.thanDxmax  = float(fr.readAtt("maxdistancex")[0])
        self.thanDext   = float(fr.readAtt("extensionx")[0])
        self.thanCena   = tuple(fr.readSnode("centroid", 3))
        self.thanNori   = int(fr.readAtt("Original_segments")[0])
        it = fr.iterNodes()
        for ca in it:
            cb = next(it)
            self.thanLines.append((tuple(ca), tuple(cb)))
