Changeset 4649

Timestamp:

Aug 1, 2007, 10:22:11 AM (17 years ago)

Author:

duncan

Message:

Optimising search_functions.py

Location:

anuga_core/source/anuga/fit_interpolate

Files:

• benchmark_least_squares.py (modified) (2 diffs)
• search_functions.py (modified) (2 diffs)
• ticket178_benchmark.py (modified) (4 diffs)

anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py

@@ -91 +91 @@
         #m0 = None on windows
         m0 = mem_usage()
-        
-        profile_file = "P" + str(num_of_points) + \
+        if is_fit is True:
+            op = "Fit_"
+        else:
+            op = "Interp_"
+        profile_file = op + "P" + str(num_of_points) + \
                        "T" + str(len(mesh_dict['triangles'])) + \
                        "PPC" + str(max_points_per_cell) + \
@@ -160 +163 @@
                                      mesh_dict['triangles'], 
                                  max_vertices_per_cell = max_points_per_cell)
-                s = """calc = interp.interpolate(mesh_dict['vertex_attributes']
+                
+                if run_profile:
+                    s="""calc=interp.interpolate(mesh_dict['vertex_attributes']
+                    ,points_dict['points'],start_blocking_len=blocking_len)"""
+                    pobject = profile.Profile()
+                    presult = pobject.runctx(s,
+                                             vars(sys.modules[__name__]),
+                                             vars())
+                    prof_file = tempfile.mktemp(".prof")
+                    presult.dump_stats(prof_file)
+                    #
+                    # Let process these results
+                    S = pstats.Stats(prof_file)
+                    saveout = sys.stdout 
+                    pfile = open(profile_file, "w")
+                    sys.stdout = pfile
+                    s = S.sort_stats('cumulative').print_stats(30)
+                    sys.stdout = saveout 
+                    pfile.close()
+                    os.remove(prof_file)
+                    
+                else:
+                    calc = interp.interpolate(mesh_dict['vertex_attributes']
                                           ,points_dict['points']
-                                          ,start_blocking_len=blocking_len)"""
-                
-                fileName = tempfile.mktemp(".prof")
-                profile.run(s, fileName) #profile_file)
-                
-                S = pstats.Stats(fileName)
-                s = S.sort_stats('cumulative').print_stats(30)
-                print "***********"
-                print s
-                print "***********"
-                pfile = file.open(profile_file, "w")
-                pfile.write(s)
-                pfile.close()
+                                          ,start_blocking_len = 500000)
+                
         time_taken_sec = (time.time()-t0)
         m1 = mem_usage()

anuga_core/source/anuga/fit_interpolate/search_functions.py

@@ -97 +97 @@
 
         # Find triangle that contains x (if any) and interpolate
-        element_found, sigma0, sigma1, sigma2, k =\
-                       find_triangle_compute_interpolation(mesh,
-                                                           triangle_list,
-                                                           x)
+        
+        for k, _ in triangle_list:
+            element_found, sigma0, sigma1, sigma2, k =\
+                           find_triangle_compute_interpolation(mesh,
+                                                               k,
+                                                               x)
+
+            if element_found is True:
+                # Don't look for any other triangles in the triangle list
+                break
 
         if element_found is True:
-            # Don't look for any other triangle           
+            # Don't look for any other triangle_lists from the
+            # candidate_vertices
             break
-
         
     return element_found, sigma0, sigma1, sigma2, k
@@ -111 +117 @@
 
             
-def find_triangle_compute_interpolation(mesh, triangle_list, x):
+def find_triangle_compute_interpolation(mesh, k, x):
     """Compute linear interpolation of point x and triangle k in mesh.
     It is assumed that x belongs to triangle k.
     """
 
-    element_found = False
-    for k, _ in triangle_list:
-        # Get the three vertex_points of candidate triangle k
-        xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)            
-        
-        # Get the three normals 
-        n0 = mesh.get_normal(k, 0)
-        n1 = mesh.get_normal(k, 1)
-        n2 = mesh.get_normal(k, 2)            
+    # Get the three vertex_points of candidate triangle k
+    xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)
+
+    # this is where we can call some fast c code.
+    xmax = max(xi0[0], xi1[0], xi2[0])
+    xmin = min(xi0[0], xi1[0], xi2[0])
+    ymax = max(xi0[1], xi1[1], xi2[1])
+    ymin = min(xi0[1], xi1[1], xi2[1])
+
+    # Integrity check - machine precision is too hard
+    # so we use hardwired single precision 
+    epsilon = 1.0e-6
+    
+    if  x[0] > xmax + epsilon:
+        return False,0,0,0,0
+    if  x[0] < xmin - epsilon:
+        return False,0,0,0,0
+    if  x[1] > ymax + epsilon:
+        return False,0,0,0,0
+    if  x[1] < ymin - epsilon:
+        return False,0,0,0,0
+    
+    
+    # Get the three normals 
+    n0 = mesh.get_normal(k, 0)
+    n1 = mesh.get_normal(k, 1)
+    n2 = mesh.get_normal(k, 2)            
         
         
-        # Compute interpolation
-        sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
-        sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
-        sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
+    # Compute interpolation
+    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
+    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
+    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
 
-        # Integrity check - machine precision is too hard
-        # so we use hardwired single precision 
-        epsilon = 1.0e-6
-        delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
-        msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
-              %(delta, epsilon)
-        assert delta < epsilon, msg
+    delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
+    msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
+          %(delta, epsilon)
+    assert delta < epsilon, msg
 
 
-        # Check that this triangle contains the data point
-        # Sigmas are allowed to get negative within
-        # machine precision on some machines (e.g. nautilus)
-        epsilon = get_machine_precision() * 2
-        if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
-            element_found = True
-            break
-            
+    # Check that this triangle contains the data point
+    # Sigmas are allowed to get negative within
+    # machine precision on some machines (e.g. nautilus)
+    epsilon = get_machine_precision() * 2
+    if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
+        element_found = True
+    else:
+        element_found = False 
     return element_found, sigma0, sigma1, sigma2, k

anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py

@@ -14 +14 @@
 ofile = 'lbm_resultsII.csv'
 delimiter = ','
-
+run_profile = True
 use_least_squares_list = [False]
-is_fit_list = [True]
-num_of_points_list = [200, 600, 2000, 6000, 10000, 20000] 
+is_fit_list = [True, False]
+num_of_points_list = [3, 200, 600, 2000, 6000, 10000, 20000] 
 maxArea_list = [ 0.008, 0.0016, 0.0008]
 max_points_per_cell_list = [2,4,8,16,30,64]
 use_file_type_list = ['pts']
-
-
-num_of_points_list = [200] 
+#num_of_points_list = [3, 200]
 maxArea_list = [ 0.008]
 max_points_per_cell_list = [30]
@@ -36 +34 @@
          "max_points_per_cell" + delimiter +
          "is_fit" + delimiter +
+         "is_profiling" + delimiter +
          "mem"  + delimiter +
          "time" + delimiter + "\n")
@@ -53 +52 @@
                                                    ,use_file_type=use_file_type
                                                    ,save=True
-                                                   ,run_profile=True
+                                                   ,run_profile=run_profile
                                                )
                     print "time",time
@@ -63 +62 @@
                              str(max_points_per_cell) + delimiter +
                              str(is_fit) + delimiter +
+                             str(run_profile) + delimiter +
                              str(mem)  + delimiter +
                              str(time) + delimiter + "\n")