Changeset 1221

Timestamp:

Apr 12, 2005, 3:43:56 PM (20 years ago)

Author:

prow

Message:

fixed. point_on_line needed an allclose.

Location:

inundation/ga/storm_surge/pmesh

Files:

• mesh.py (modified) (16 diffs)
• pmesh.py (modified) (1 diff)

inundation/ga/storm_surge/pmesh/mesh.py

@@ -2179 +2179 @@
     def triangles_to_polySet(self,setName):
         seg2line = self.seg2line
-        ver2point = self.ver2point
+        ver2point= self.ver2point
         line2seg = self.line2seg
-        point2ver = self.point2ver
+        point2ver= self.point2ver
 
         from Numeric import array,allclose
@@ -2189 +2189 @@
         for triangle in Triangles:
             Triangles_dict[triangle]=None
+ 
 
         point_keys = {}
@@ -2198 +2199 @@
                 userVertices[vertex]=True
 
+        #######
+        for vert in userVertices.keys():
+            assert point_keys[(vert.x,vert.y)]==vert
+        assert len(point_keys.keys())==len(userVertices.keys())
+        #######
+
+
         userSegments = self.segs_to_dict(self.userSegments)
         affine_lines = Affine_Linespace()
@@ -2214 +2222 @@
                     a = triangle.vertices[i-1]
                     b = triangle.vertices[i-2]
-                    for vert in (a,b):
-                        point = ver2point(vert)
-                        if not point_keys.has_key(point):
-                            userVertices[vert]=True
-                            point_keys[point]=vert
-
                     #Get possible matches:
                     point_a = ver2point(a)
@@ -2230 +2232 @@
                     for user_line in possible_lines:
                         if self.point_on_line(midpoint,user_line):
-                            if user_line in userSegments:
+                            if userSegments.has_key(user_line):
                                 parent_segment = userSegments.pop(user_line)
                             else:
@@ -2257 +2259 @@
                         for point in line:
                             if point_keys.has_key(point):
-                                vert = point_keys[point]
+                                vert = point_keys[point]
                             else:
                                 vert = Vertex(point[0],point[1])
@@ -2267 +2269 @@
                         affine_lines.append(line)
         
+        #######
+        for vert in userVertices.keys():
+            assert point_keys[(vert.x,vert.y)]==vert
+        assert len(point_keys.keys())==len(userVertices.keys())
+        #######
         self.userVerticies = []
         self.userSegments = []
@@ -2277 +2284 @@
         #FIXME change alpha to user and sync with pmesh
         #ie self.alphaSegments.append(alphaSegments[key])
-
+        #print point_keys.keys()
+        #print userVertices.keys()
+        #print userVertices.keys()
         return userVertices,userSegments,alphaSegments
 
@@ -2343 +2352 @@
         y0=line[0][1]
         y1=line[1][1]
-        return point_on_line(x,y,x0,y0,x1,y1)
-        #a wrapper for the c func.
+        from Numeric import array, dot, allclose
+        from math import sqrt
+
+        a = array([x - x0, y - y0]) 
+        a_normal = array([a[1], -a[0]])
+
+        b = array([x1 - x0, y1 - y0])
+    
+        if allclose(dot(a_normal, b),0):
+            #Point is somewhere on the infinite extension of the line
+
+            len_a = sqrt(sum(a**2))            
+            len_b = sqrt(sum(b**2))                             
+            if dot(a, b) >= 0 and len_a <= len_b:
+               return True
+            else:   
+               return False
+        else:
+          return False  
 
     def line_length(self,line):
@@ -3122 +3148 @@
     #NOT IMPLIMENTED
     at a.precision = 2:
-        a.round_up[0.0]=
+        a.round_up_rel[0.0]=
         a.round_flat[0.0]=
-        a.round_downp0.0]=
+        a.round_down_rel[0.0]=
 
         a.up((0.1,2.04))=
 
-    If tolerance = 0, nothing gets rounded into
-    two bins. If tolerance = 0.5, everything does.
+    If t_rel = 0, nothing gets rounded into
+    two bins. If t_rel = 0.5, everything does.
 
     Ideally, precision can be set high, so that multiple
-    entries are rarely in the same bin. And tolerance should
+    entries are rarely in the same bin. And t_rel should
     be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
     so that it is rare to put items in mutiple bins.
 
     Ex bins per entry = product(a,b,c...,n)
-    a = 1 or 2 s.t. Ex(a) = 1+2*tolerance
+    a = 1 or 2 s.t. Ex(a) = 1+2*t_rel
     b = 1 or 2 ... 
 
     BUT!!! to avoid missing the right bin:
-    (-10)**(precision+1)*tolerance must be greater than the 
+    (-10)**(precision+1)*t_rel must be greater than the 
     greatest possible variation that an identical element
     can display.
-    """
-    def __init__(self,precision = 6,tolerance = 0.01):
-        self.__precision__ = precision
-        self.__tolerance__ = tolerance
-        assert tolerance <= 0.5
+
+
+    Note that if tol = 0.5 (the max allowed) 0.6 will round to .7 and .5
+    but not .6 - this looks wrong, but note that *everything* will round,
+    so .6 wont be missed as everything close to it will check in .7 and .5.
+    """
+    def __init__(self,p_rel = 6,t_rel = 0.01):
+        self.__p_rel__ = p_rel
+        self.__t_rel__ = t_rel
+
+        self.__p_abs__ = p_rel+1
+        self.__t_abs__ = t_rel
+
+        assert t_rel <= 0.5
         self.__items__ = {}
         from math import frexp
@@ -3164 +3199 @@
     def __rounded_keys__(self,key):
         keys = []
-        up_key=list(key)
-        for i in range(len(up_key)):
-            up_key[i]=self.round_up(key[i])
-        keys.append(tuple(up_key))
-
-        #round down the value
-        for i in range(len(up_key)):
-            if not keys[0][i]==self.round_down(key[i]):
-            #unless round_up = round_down 
-            #so the keys stay unique
-                for j in range(len(keys)):
-                    down_key = list(keys[j])
-                    down_key[i]=self.round_down(key[i])
-                    down_key=tuple(down_key)
-                    keys.append(down_key)
+        rounded_key=list(key)
+        rounded_values=list(key)
+
+        roundings = [self.round_up_rel,\
+        self.round_down_rel,self.round_flat_rel,\
+        self.round_down_abs,self.round_up_abs,\
+        self.round_flat_abs]
+
+
+        roundings = \
+        [self.round_up_rel, self.round_down_rel, \
+         self.round_up_rel2,self.round_down_rel2,\
+         self.round_up_abs, self.round_down_abs]
+
+        #initialise rounded_values
+        round = roundings.pop(0)
+        for i in range(len(rounded_values)):
+            rounded_key[i]=round(key[i])
+            rounded_values[i]={}
+            rounded_values[i][rounded_key[i]]=None
+        keys.append(tuple(rounded_key))
+        
+        for round in roundings:
+            for i in range(len(rounded_key)):
+                rounded_value=round(key[i])
+                if not rounded_values[i].has_key(rounded_value):
+                    #ie unless round_up_rel = round_down_rel
+                    #so the keys stay unique
+                    for j in range(len(keys)):
+                        rounded_key = list(keys[j])
+                        rounded_key[i]=rounded_value
+                        keys.append(tuple(rounded_key))
         return keys
 
@@ -3230 +3282 @@
         return self.__contains__(item)
 
-    def round_up(self,value):
-         tolerance=self.__tolerance__
+    def round_up_rel2(self,value):
+         t_rel=self.__t_rel__
+         #Rounding up the value
+         m,e = self.frexp(value)
+         m = m/2
+         e = e + 1
+         #m is the mantissa, e the exponent
+         # 0.5 < |m| < 1.0
+         m = m+t_rel*(10**-(self.__p_rel__))
+         #bump m up
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_down_rel2(self,value):
+         t_rel=self.__t_rel__
+         #Rounding down the value
+         m,e = self.frexp(value)
+         m = m/2
+         e = e + 1
+         #m is the mantissa, e the exponent
+         # 0.5 < m < 1.0
+         m = m-t_rel*(10**-(self.__p_rel__))
+         #bump the |m| down, by 5% or whatever
+         #self.p_rel dictates
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_flat_rel2(self,value):
+    #redundant
+         m,e = self.frexp(value)
+         m = m/2
+         e = e + 1
+         m = round(m,self.__p_rel__)
+         return m*(2.**e)
+
+    def round_up_rel(self,value):
+         t_rel=self.__t_rel__
          #Rounding up the value
          m,e = self.frexp(value)
          #m is the mantissa, e the exponent
          # 0.5 < |m| < 1.0
-         m = m*(1+tolerance*(10**-(self.__precision__+1))/((m**2)**0.5))
-         #bump m up, by 5% or whatever
-         #self.precision dictates
-         m = round(m,self.__precision__)
-     #round m off so that the bump up 
-     #will be discernable if and only if
+         m = m+t_rel*(10**-(self.__p_rel__))
+         #bump m up
+         m = round(m,self.__p_rel__)
          return m*(2.**e)
 
-    def round_down(self,value):
-         tolerance=self.__tolerance__
+    def round_down_rel(self,value):
+         t_rel=self.__t_rel__
          #Rounding down the value
          m,e = self.frexp(value)
          #m is the mantissa, e the exponent
          # 0.5 < m < 1.0
-         m = m*(1-tolerance*(10**-(self.__precision__+1))/((m**2)**0.5))
+         m = m-t_rel*(10**-(self.__p_rel__))
          #bump the |m| down, by 5% or whatever
-         #self.precision dictates
-         m = round(m,self.__precision__)
-     #round m off so that the bump down 
-     #will be discernable if and only if
+         #self.p_rel dictates
+         m = round(m,self.__p_rel__)
          return m*(2.**e)
 
-    def round_flat(self,value):
+    def round_flat_rel(self,value):
     #redundant
          m,e = self.frexp(value)
-         m = round(m,self.__precision__)
+         m = round(m,self.__p_rel__)
          return m*(2.**e)
+
+    def round_up_abs(self,value):
+         t_abs=self.__t_abs__
+         #Rounding up the value
+         m = value+t_abs*(10**-(self.__p_abs__))
+         #bump m up
+         m = round(m,self.__p_abs__)
+         return m
+
+    def round_down_abs(self,value):
+         t_abs=self.__t_abs__
+         #Rounding down the value
+         m = value-t_abs*(10**-(self.__p_abs__))
+         #bump the |m| down, by 5% or whatever
+         #self.p_rel dictates
+         m = round(m,self.__p_abs__)
+         return m
+
+    def round_flat_abs(self,value):
+    #redundant?
+         m = round(value,self.__p_abs__)
+         return m
 
     def keys(self):
@@ -3289 +3393 @@
     a[elephant] -> []
     """
-    def __init__(self,mapping,precision = 6, tolerance=0.01):
+    def __init__(self,mapping,p_rel = 6, t_rel=0.01):
         Discretised_Tuple_Set.__init__\
-         (self,precision,tolerance = tolerance)
+         (self,p_rel,t_rel = t_rel)
         self.mapping = mapping
 
@@ -3308 +3412 @@
     Nearly.
     """
-    def __init__(self,precision = 4,tolerance=0.04):
+    def __init__(self,p_rel = 3,t_rel=0.1):
         Mapped_Discretised_Tuple_Set.__init__\
             (self,self.affine_line,\
-            precision=precision,tolerance=tolerance)
+            p_rel=p_rel,t_rel=t_rel)
 
     def affine_line(self,line):
@@ -3357 +3461 @@
         else:
             sign_a = a/((a**2)**0.5)
+        #This does not change the mathematical
+        #properties, but it makes comparism
         if b == 0:
             sign_b = 1.
@@ -3370 +3476 @@
         return a,b,c
 
+
     
 if __name__ == "__main__":

inundation/ga/storm_surge/pmesh/pmesh.py

@@ -354 +354 @@
         print 'len(alphaSegments.keys())'
         print len(alphaSegments.keys())
+
+
+        #######
+        point_keys = {}
+        for vert in userVertices.keys():
+            assert not point_keys.has_key((vert.x,vert.y))
+            point_keys[(vert.x,vert.y)]=vert
+        assert len(point_keys.keys())==len(userVertices.keys())
+        #######
+
         for v in userVertices.keys():
             if userVertices[v] is True: