Changeset 4649


Ignore:
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:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/benchmark_least_squares.py

    r4648 r4649  
    9191        #m0 = None on windows
    9292        m0 = mem_usage()
    93        
    94         profile_file = "P" + str(num_of_points) + \
     93        if is_fit is True:
     94            op = "Fit_"
     95        else:
     96            op = "Interp_"
     97        profile_file = op + "P" + str(num_of_points) + \
    9598                       "T" + str(len(mesh_dict['triangles'])) + \
    9699                       "PPC" + str(max_points_per_cell) + \
     
    160163                                     mesh_dict['triangles'],
    161164                                 max_vertices_per_cell = max_points_per_cell)
    162                 s = """calc = interp.interpolate(mesh_dict['vertex_attributes']
     165               
     166                if run_profile:
     167                    s="""calc=interp.interpolate(mesh_dict['vertex_attributes']
     168                    ,points_dict['points'],start_blocking_len=blocking_len)"""
     169                    pobject = profile.Profile()
     170                    presult = pobject.runctx(s,
     171                                             vars(sys.modules[__name__]),
     172                                             vars())
     173                    prof_file = tempfile.mktemp(".prof")
     174                    presult.dump_stats(prof_file)
     175                    #
     176                    # Let process these results
     177                    S = pstats.Stats(prof_file)
     178                    saveout = sys.stdout
     179                    pfile = open(profile_file, "w")
     180                    sys.stdout = pfile
     181                    s = S.sort_stats('cumulative').print_stats(30)
     182                    sys.stdout = saveout
     183                    pfile.close()
     184                    os.remove(prof_file)
     185                   
     186                else:
     187                    calc = interp.interpolate(mesh_dict['vertex_attributes']
    163188                                          ,points_dict['points']
    164                                           ,start_blocking_len=blocking_len)"""
    165                
    166                 fileName = tempfile.mktemp(".prof")
    167                 profile.run(s, fileName) #profile_file)
    168                
    169                 S = pstats.Stats(fileName)
    170                 s = S.sort_stats('cumulative').print_stats(30)
    171                 print "***********"
    172                 print s
    173                 print "***********"
    174                 pfile = file.open(profile_file, "w")
    175                 pfile.write(s)
    176                 pfile.close()
     189                                          ,start_blocking_len = 500000)
     190               
    177191        time_taken_sec = (time.time()-t0)
    178192        m1 = mem_usage()
  • anuga_core/source/anuga/fit_interpolate/search_functions.py

    r4599 r4649  
    9797
    9898        # Find triangle that contains x (if any) and interpolate
    99         element_found, sigma0, sigma1, sigma2, k =\
    100                        find_triangle_compute_interpolation(mesh,
    101                                                            triangle_list,
    102                                                            x)
     99       
     100        for k, _ in triangle_list:
     101            element_found, sigma0, sigma1, sigma2, k =\
     102                           find_triangle_compute_interpolation(mesh,
     103                                                               k,
     104                                                               x)
     105
     106            if element_found is True:
     107                # Don't look for any other triangles in the triangle list
     108                break
    103109
    104110        if element_found is True:
    105             # Don't look for any other triangle           
     111            # Don't look for any other triangle_lists from the
     112            # candidate_vertices
    106113            break
    107 
    108114       
    109115    return element_found, sigma0, sigma1, sigma2, k
     
    111117
    112118           
    113 def find_triangle_compute_interpolation(mesh, triangle_list, x):
     119def find_triangle_compute_interpolation(mesh, k, x):
    114120    """Compute linear interpolation of point x and triangle k in mesh.
    115121    It is assumed that x belongs to triangle k.
    116122    """
    117123
    118     element_found = False
    119     for k, _ in triangle_list:
    120         # Get the three vertex_points of candidate triangle k
    121         xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)           
    122        
    123         # Get the three normals
    124         n0 = mesh.get_normal(k, 0)
    125         n1 = mesh.get_normal(k, 1)
    126         n2 = mesh.get_normal(k, 2)           
     124    # Get the three vertex_points of candidate triangle k
     125    xi0, xi1, xi2 = mesh.get_vertex_coordinates(triangle_id=k)
     126
     127    # this is where we can call some fast c code.
     128    xmax = max(xi0[0], xi1[0], xi2[0])
     129    xmin = min(xi0[0], xi1[0], xi2[0])
     130    ymax = max(xi0[1], xi1[1], xi2[1])
     131    ymin = min(xi0[1], xi1[1], xi2[1])
     132
     133    # Integrity check - machine precision is too hard
     134    # so we use hardwired single precision
     135    epsilon = 1.0e-6
     136   
     137    if  x[0] > xmax + epsilon:
     138        return False,0,0,0,0
     139    if  x[0] < xmin - epsilon:
     140        return False,0,0,0,0
     141    if  x[1] > ymax + epsilon:
     142        return False,0,0,0,0
     143    if  x[1] < ymin - epsilon:
     144        return False,0,0,0,0
     145   
     146   
     147    # Get the three normals
     148    n0 = mesh.get_normal(k, 0)
     149    n1 = mesh.get_normal(k, 1)
     150    n2 = mesh.get_normal(k, 2)           
    127151       
    128152       
    129         # Compute interpolation
    130         sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
    131         sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
    132         sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
     153    # Compute interpolation
     154    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
     155    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
     156    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
    133157
    134         # Integrity check - machine precision is too hard
    135         # so we use hardwired single precision
    136         epsilon = 1.0e-6
    137         delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
    138         msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
    139               %(delta, epsilon)
    140         assert delta < epsilon, msg
     158    delta = abs(sigma0+sigma1+sigma2-1.0) # Should be close to zero
     159    msg = 'abs(sigma0+sigma1+sigma2-1) = %.15e, eps = %.15e'\
     160          %(delta, epsilon)
     161    assert delta < epsilon, msg
    141162
    142163
    143         # Check that this triangle contains the data point
    144         # Sigmas are allowed to get negative within
    145         # machine precision on some machines (e.g. nautilus)
    146         epsilon = get_machine_precision() * 2
    147         if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
    148             element_found = True
    149             break
    150            
     164    # Check that this triangle contains the data point
     165    # Sigmas are allowed to get negative within
     166    # machine precision on some machines (e.g. nautilus)
     167    epsilon = get_machine_precision() * 2
     168    if sigma0 >= -epsilon and sigma1 >= -epsilon and sigma2 >= -epsilon:
     169        element_found = True
     170    else:
     171        element_found = False
    151172    return element_found, sigma0, sigma1, sigma2, k
  • anuga_core/source/anuga/fit_interpolate/ticket178_benchmark.py

    r4648 r4649  
    1414ofile = 'lbm_resultsII.csv'
    1515delimiter = ','
    16 
     16run_profile = True
    1717use_least_squares_list = [False]
    18 is_fit_list = [True]
    19 num_of_points_list = [200, 600, 2000, 6000, 10000, 20000]
     18is_fit_list = [True, False]
     19num_of_points_list = [3, 200, 600, 2000, 6000, 10000, 20000]
    2020maxArea_list = [ 0.008, 0.0016, 0.0008]
    2121max_points_per_cell_list = [2,4,8,16,30,64]
    2222use_file_type_list = ['pts']
    23 
    24 
    25 num_of_points_list = [200]
     23#num_of_points_list = [3, 200]
    2624maxArea_list = [ 0.008]
    2725max_points_per_cell_list = [30]
     
    3634         "max_points_per_cell" + delimiter +
    3735         "is_fit" + delimiter +
     36         "is_profiling" + delimiter +
    3837         "mem"  + delimiter +
    3938         "time" + delimiter + "\n")
     
    5352                                                   ,use_file_type=use_file_type
    5453                                                   ,save=True
    55                                                    ,run_profile=True
     54                                                   ,run_profile=run_profile
    5655                                               )
    5756                    print "time",time
     
    6362                             str(max_points_per_cell) + delimiter +
    6463                             str(is_fit) + delimiter +
     64                             str(run_profile) + delimiter +
    6565                             str(mem)  + delimiter +
    6666                             str(time) + delimiter + "\n")
Note: See TracChangeset for help on using the changeset viewer.