source: inundation/pyvolution/least_squares.py @ 2683

Last change on this file since 2683 was 2683, checked in by nick, 18 years ago

removed print statements from least_squares.py and mesh.py

File size: 51.8 KB
Line 
1"""Least squares smooting and interpolation.
2
3   Implements a penalised least-squares fit and associated interpolations.
4
5   The penalty term (or smoothing term) is controlled by the smoothing
6   parameter alpha.
7   With a value of alpha=0, the fit function will attempt
8   to interpolate as closely as possible in the least-squares sense.
9   With values alpha > 0, a certain amount of smoothing will be applied.
10   A positive alpha is essential in cases where there are too few
11   data points.
12   A negative alpha is not allowed.
13   A typical value of alpha is 1.0e-6
14
15
16   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
17   Geoscience Australia, 2004.
18"""
19
20import exceptions
21class ShapeError(exceptions.Exception): pass
22class FittingError(exceptions.Exception): pass
23
24
25#from general_mesh import General_mesh
26from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, ArrayType, NewAxis
27from pyvolution.mesh import Mesh
28
29from Numeric import zeros, take, compress, array, Float, Int, dot, transpose, concatenate, ArrayType
30from utilities.sparse import Sparse, Sparse_CSR
31from utilities.cg_solve import conjugate_gradient, VectorShapeError
32from utilities.numerical_tools import ensure_numeric, mean, gradient
33
34
35from coordinate_transforms.geo_reference import Geo_reference
36
37import time
38
39
40
41DEFAULT_ALPHA = 0.001
42
43def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
44                     alpha=DEFAULT_ALPHA, verbose= False,
45                     expand_search = False,
46                     data_origin = None,
47                     mesh_origin = None,
48                     precrop = False,
49                     display_errors = True):
50    """
51    Given a mesh file (tsh) and a point attribute file (xya), fit
52    point attributes to the mesh and write a mesh file with the
53    results.
54
55
56    If data_origin is not None it is assumed to be
57    a 3-tuple with geo referenced
58    UTM coordinates (zone, easting, northing)
59
60    NOTE: Throws IOErrors, for a variety of file problems.
61   
62    mesh_origin is the same but refers to the input tsh file.
63    FIXME: When the tsh format contains it own origin, these parameters can go.
64    FIXME: And both origins should be obtained from the specified files.
65    """
66
67    from load_mesh.loadASCII import import_mesh_file, \
68                 import_points_file, export_mesh_file, \
69                 concatinate_attributelist
70
71
72    try:
73        mesh_dict = import_mesh_file(mesh_file)
74    except IOError,e:
75        if display_errors:
76            print "Could not load bad file. ", e
77        raise IOError  #Re-raise exception
78       
79    vertex_coordinates = mesh_dict['vertices']
80    triangles = mesh_dict['triangles']
81    if type(mesh_dict['vertex_attributes']) == ArrayType:
82        old_point_attributes = mesh_dict['vertex_attributes'].tolist()
83    else:
84        old_point_attributes = mesh_dict['vertex_attributes']
85
86    if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
87        old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
88    else:
89        old_title_list = mesh_dict['vertex_attribute_titles']
90
91    if verbose: print 'tsh file %s loaded' %mesh_file
92
93    # load in the .pts file
94    try:
95        point_dict = import_points_file(point_file, verbose=verbose)
96    except IOError,e:
97        if display_errors:
98            print "Could not load bad file. ", e
99        raise IOError  #Re-raise exception 
100
101    point_coordinates = point_dict['pointlist']
102    title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
103
104    if point_dict.has_key('geo_reference') and not point_dict['geo_reference'] is None:
105        data_origin = point_dict['geo_reference'].get_origin()
106    else:
107        data_origin = (56, 0, 0) #FIXME(DSG-DSG)
108
109    if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
110        mesh_origin = mesh_dict['geo_reference'].get_origin()
111    else:
112        mesh_origin = (56, 0, 0) #FIXME(DSG-DSG)
113
114    if verbose: print "points file loaded"
115    if verbose: print "fitting to mesh"
116    f = fit_to_mesh(vertex_coordinates,
117                    triangles,
118                    point_coordinates,
119                    point_attributes,
120                    alpha = alpha,
121                    verbose = verbose,
122                    expand_search = expand_search,
123                    data_origin = data_origin,
124                    mesh_origin = mesh_origin,
125                    precrop = precrop)
126    if verbose: print "finished fitting to mesh"
127
128    # convert array to list of lists
129    new_point_attributes = f.tolist()
130    #FIXME have this overwrite attributes with the same title - DSG
131    #Put the newer attributes last
132    if old_title_list <> []:
133        old_title_list.extend(title_list)
134        #FIXME can this be done a faster way? - DSG
135        for i in range(len(old_point_attributes)):
136            old_point_attributes[i].extend(new_point_attributes[i])
137        mesh_dict['vertex_attributes'] = old_point_attributes
138        mesh_dict['vertex_attribute_titles'] = old_title_list
139    else:
140        mesh_dict['vertex_attributes'] = new_point_attributes
141        mesh_dict['vertex_attribute_titles'] = title_list
142
143    #FIXME (Ole): Remember to output mesh_origin as well
144    if verbose: print "exporting to file ", mesh_output_file
145
146    try:
147        export_mesh_file(mesh_output_file, mesh_dict)
148    except IOError,e:
149        if display_errors:
150            print "Could not write file. ", e
151        raise IOError
152
153def fit_to_mesh(vertex_coordinates,
154                triangles,
155                point_coordinates,
156                point_attributes,
157                alpha = DEFAULT_ALPHA,
158                verbose = False,
159                acceptable_overshoot = 1.01,
160                expand_search = False,
161                data_origin = None,
162                mesh_origin = None,
163                precrop = False,
164                use_cache = False):
165    """
166    Fit a smooth surface to a triangulation,
167    given data points with attributes.
168
169
170        Inputs:
171
172          vertex_coordinates: List of coordinate pairs [xi, eta] of points
173          constituting mesh (or a an m x 2 Numeric array)
174
175          triangles: List of 3-tuples (or a Numeric array) of
176          integers representing indices of all vertices in the mesh.
177
178          point_coordinates: List of coordinate pairs [x, y] of data points
179          (or an nx2 Numeric array)
180
181          alpha: Smoothing parameter.
182
183          acceptable overshoot: controls the allowed factor by which fitted values
184          may exceed the value of input data. The lower limit is defined
185          as min(z) - acceptable_overshoot*delta z and upper limit
186          as max(z) + acceptable_overshoot*delta z
187         
188
189          point_attributes: Vector or array of data at the point_coordinates.
190
191          data_origin and mesh_origin are 3-tuples consisting of
192          UTM zone, easting and northing. If specified
193          point coordinates and vertex coordinates are assumed to be
194          relative to their respective origins.
195
196    """
197
198    if use_cache is True:
199        from caching.caching import cache
200        interp = cache(_interpolation,
201                       (vertex_coordinates,
202                        triangles,
203                        point_coordinates),
204                       {'alpha': alpha,
205                        'verbose': verbose,
206                        'expand_search': expand_search,
207                        'data_origin': data_origin,
208                        'mesh_origin': mesh_origin,
209                        'precrop': precrop},
210                       verbose = verbose)       
211       
212    else:
213        interp = Interpolation(vertex_coordinates,
214                               triangles,
215                               point_coordinates,
216                               alpha = alpha,
217                               verbose = verbose,
218                               expand_search = expand_search,
219                               data_origin = data_origin,
220                               mesh_origin = mesh_origin,
221                               precrop = precrop)
222
223    vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
224
225
226    #Sanity check
227    point_coordinates = ensure_numeric(point_coordinates)
228    vertex_coordinates = ensure_numeric(vertex_coordinates)
229
230    #Data points
231    X = point_coordinates[:,0]
232    Y = point_coordinates[:,1] 
233    Z = ensure_numeric(point_attributes)
234    if len(Z.shape) == 1:
235        Z = Z[:, NewAxis]
236       
237
238    #Data points inside mesh boundary
239    indices = interp.point_indices
240    if indices is not None:   
241        Xc = take(X, indices)
242        Yc = take(Y, indices)   
243        Zc = take(Z, indices)
244    else:
245        Xc = X
246        Yc = Y 
247        Zc = Z       
248   
249    #Vertex coordinates
250    Xi = vertex_coordinates[:,0]
251    Eta = vertex_coordinates[:,1]       
252    Zeta = ensure_numeric(vertex_attributes)
253    if len(Zeta.shape) == 1:
254        Zeta = Zeta[:, NewAxis]   
255
256    for i in range(Zeta.shape[1]): #For each attribute
257        zeta = Zeta[:,i]
258        z = Z[:,i]               
259        zc = Zc[:,i]
260
261        max_zc = max(zc)
262        min_zc = min(zc)
263        delta_zc = max_zc-min_zc
264        upper_limit = max_zc + delta_zc*acceptable_overshoot
265        lower_limit = min_zc - delta_zc*acceptable_overshoot       
266       
267
268        if max(zeta) > upper_limit or min(zeta) < lower_limit:
269            msg = 'Least sqares produced values outside the allowed '
270            msg += 'range [%f, %f].\n' %(lower_limit, upper_limit)
271            msg += 'z in [%f, %f], zeta in [%f, %f].\n' %(min_zc, max_zc,
272                                                          min(zeta), max(zeta))
273            msg += 'If greater range is needed, increase the value of '
274            msg += 'acceptable_fit_overshoot (currently %.2f).\n' %(acceptable_overshoot)
275
276
277            offending_vertices = (zeta > upper_limit or zeta < lower_limit)
278            Xi_c = compress(offending_vertices, Xi)
279            Eta_c = compress(offending_vertices, Eta)
280            offending_coordinates = concatenate((Xi_c[:, NewAxis],
281                                                 Eta_c[:, NewAxis]),
282                                                axis=1)
283
284            msg += 'Offending locations:\n %s' %(offending_coordinates)
285           
286            raise FittingError, msg
287
288
289   
290        if verbose:
291            print '+------------------------------------------------'
292            print 'Least squares statistics'
293            print '+------------------------------------------------'   
294            print 'points: %d points' %(len(z))
295            print '    x in [%f, %f]'%(min(X), max(X))
296            print '    y in [%f, %f]'%(min(Y), max(Y))
297            print '    z in [%f, %f]'%(min(z), max(z))
298            print
299
300            if indices is not None:
301                print 'Cropped points: %d points' %(len(zc))
302                print '    x in [%f, %f]'%(min(Xc), max(Xc))
303                print '    y in [%f, %f]'%(min(Yc), max(Yc))
304                print '    z in [%f, %f]'%(min(zc), max(zc))
305                print
306           
307
308            print 'Mesh: %d vertices' %(len(zeta))
309            print '    xi in [%f, %f]'%(min(Xi), max(Xi))
310            print '    eta in [%f, %f]'%(min(Eta), max(Eta))
311            print '    zeta in [%f, %f]'%(min(zeta), max(zeta))
312            print '+------------------------------------------------'
313
314    return vertex_attributes
315
316
317
318def pts2rectangular(pts_name, M, N, alpha = DEFAULT_ALPHA,
319                    verbose = False, reduction = 1):
320    """Fits attributes from pts file to MxN rectangular mesh
321
322    Read pts file and create rectangular mesh of resolution MxN such that
323    it covers all points specified in pts file.
324
325    FIXME: This may be a temporary function until we decide on
326    netcdf formats etc
327
328    FIXME: Uses elevation hardwired
329    """
330
331    import  mesh_factory
332    from load_mesh.loadASCII import import_points_file
333   
334    if verbose: print 'Read pts'
335    points_dict = import_points_file(pts_name)
336    #points, attributes = util.read_xya(pts_name)
337
338    #Reduce number of points a bit
339    points = points_dict['pointlist'][::reduction]
340    elevation = points_dict['attributelist']['elevation']  #Must be elevation
341    elevation = elevation[::reduction]
342
343    if verbose: print 'Got %d data points' %len(points)
344
345    if verbose: print 'Create mesh'
346    #Find extent
347    max_x = min_x = points[0][0]
348    max_y = min_y = points[0][1]
349    for point in points[1:]:
350        x = point[0]
351        if x > max_x: max_x = x
352        if x < min_x: min_x = x
353        y = point[1]
354        if y > max_y: max_y = y
355        if y < min_y: min_y = y
356
357    #Create appropriate mesh
358    vertex_coordinates, triangles, boundary =\
359         mesh_factory.rectangular(M, N, max_x-min_x, max_y-min_y,
360                                (min_x, min_y))
361
362    #Fit attributes to mesh
363    vertex_attributes = fit_to_mesh(vertex_coordinates,
364                        triangles,
365                        points,
366                        elevation, alpha=alpha, verbose=verbose)
367
368
369
370    return vertex_coordinates, triangles, boundary, vertex_attributes
371
372
373def _interpolation(*args, **kwargs):
374    """Private function for use with caching. Reason is that classes
375    may change their byte code between runs which is annoying.
376    """
377   
378    return Interpolation(*args, **kwargs)
379
380
381class Interpolation:
382
383    def __init__(self,
384                 vertex_coordinates,
385                 triangles,
386                 point_coordinates = None,
387                 alpha = None,
388                 verbose = False,
389                 expand_search = True,
390                 interp_only = False,
391                 max_points_per_cell = 30,
392                 mesh_origin = None,
393                 data_origin = None,
394                 precrop = False):
395
396
397        """ Build interpolation matrix mapping from
398        function values at vertices to function values at data points
399
400        Inputs:
401
402          vertex_coordinates: List of coordinate pairs [xi, eta] of
403          points constituting mesh (or a an m x 2 Numeric array)
404          Points may appear multiple times
405          (e.g. if vertices have discontinuities)
406
407          triangles: List of 3-tuples (or a Numeric array) of
408          integers representing indices of all vertices in the mesh.
409
410          point_coordinates: List of coordinate pairs [x, y] of
411          data points (or an nx2 Numeric array)
412          If point_coordinates is absent, only smoothing matrix will
413          be built
414
415          alpha: Smoothing parameter
416
417          data_origin and mesh_origin are 3-tuples consisting of
418          UTM zone, easting and northing. If specified
419          point coordinates and vertex coordinates are assumed to be
420          relative to their respective origins.
421
422        """
423
424        #Convert input to Numeric arrays
425        triangles = ensure_numeric(triangles, Int)
426        vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
427
428        #Build underlying mesh
429        if verbose: print 'Building mesh'
430        #self.mesh = General_mesh(vertex_coordinates, triangles,
431        #FIXME: Trying the normal mesh while testing precrop,
432        #       The functionality of boundary_polygon is needed for that
433
434        #FIXME - geo ref does not have to go into mesh.
435        # Change the point co-ords to conform to the
436        # mesh co-ords early in the code
437        if mesh_origin is None:
438            geo = None
439        else:
440            geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
441        self.mesh = Mesh(vertex_coordinates, triangles,
442                         geo_reference = geo)
443
444        if verbose: print 'Checking mesh integrity'
445        self.mesh.check_integrity()
446
447        if verbose: print 'Mesh integrity checked'
448       
449        self.data_origin = data_origin
450
451        self.point_indices = None
452
453        #Smoothing parameter
454        if alpha is None:
455            self.alpha = DEFAULT_ALPHA
456        else:   
457            self.alpha = alpha
458
459
460        if point_coordinates is not None:
461            if verbose: print 'Building interpolation matrix'
462            self.build_interpolation_matrix_A(point_coordinates,
463                                              verbose = verbose,
464                                              expand_search = expand_search,
465                                              interp_only = interp_only, 
466                                              max_points_per_cell =\
467                                              max_points_per_cell,
468                                              data_origin = data_origin,
469                                              precrop = precrop)
470        #Build coefficient matrices
471        if interp_only == False:
472            self.build_coefficient_matrix_B(point_coordinates,
473                                        verbose = verbose,
474                                        expand_search = expand_search,
475                                        max_points_per_cell =\
476                                        max_points_per_cell,
477                                        data_origin = data_origin,
478                                        precrop = precrop)
479        if verbose: print 'Finished interpolation'
480
481    def set_point_coordinates(self, point_coordinates,
482                              data_origin = None,
483                              verbose = False,
484                              precrop = True):
485        """
486        A public interface to setting the point co-ordinates.
487        """
488        if point_coordinates is not None:
489            if verbose: print 'Building interpolation matrix'
490            self.build_interpolation_matrix_A(point_coordinates,
491                                              verbose = verbose,
492                                              data_origin = data_origin,
493                                              precrop = precrop)
494        self.build_coefficient_matrix_B(point_coordinates, data_origin)
495
496    def build_coefficient_matrix_B(self, point_coordinates=None,
497                                   verbose = False, expand_search = True,
498                                   max_points_per_cell=30,
499                                   data_origin = None,
500                                   precrop = False):
501        """Build final coefficient matrix"""
502
503
504        if self.alpha <> 0:
505            if verbose: print 'Building smoothing matrix'
506            self.build_smoothing_matrix_D()
507
508        if point_coordinates is not None:
509            if self.alpha <> 0:
510                self.B = self.AtA + self.alpha*self.D
511            else:
512                self.B = self.AtA
513
514            #Convert self.B matrix to CSR format for faster matrix vector
515            self.B = Sparse_CSR(self.B)
516
517    def build_interpolation_matrix_A(self, point_coordinates,
518                                     verbose = False, expand_search = True,
519                                     max_points_per_cell=30,
520                                     data_origin = None,
521                                     precrop = False,
522                                     interp_only = False):
523        """Build n x m interpolation matrix, where
524        n is the number of data points and
525        m is the number of basis functions phi_k (one per vertex)
526
527        This algorithm uses a quad tree data structure for fast binning of data points
528        origin is a 3-tuple consisting of UTM zone, easting and northing.
529        If specified coordinates are assumed to be relative to this origin.
530
531        This one will override any data_origin that may be specified in
532        interpolation instance
533
534        """
535
536
537
538        #FIXME (Ole): Check that this function is memeory efficient.
539        #6 million datapoints and 300000 basis functions
540        #causes out-of-memory situation
541        #First thing to check is whether there is room for self.A and self.AtA
542        #
543        #Maybe we need some sort of blocking
544
545        from pyvolution.quad import build_quadtree
546        from utilities.polygon import inside_polygon
547       
548
549        if data_origin is None:
550            data_origin = self.data_origin #Use the one from
551                                           #interpolation instance
552
553        #Convert input to Numeric arrays just in case.
554        point_coordinates = ensure_numeric(point_coordinates, Float)
555
556        #Keep track of discarded points (if any).
557        #This is only registered if precrop is True
558        self.cropped_points = False
559
560        #Shift data points to same origin as mesh (if specified)
561
562        #FIXME this will shift if there was no geo_ref.
563        #But all this should be removed anyhow.
564        #change coords before this point
565        mesh_origin = self.mesh.geo_reference.get_origin()
566        if point_coordinates is not None:
567            if data_origin is not None:
568                if mesh_origin is not None:
569
570                    #Transformation:
571                    #
572                    #Let x_0 be the reference point of the point coordinates
573                    #and xi_0 the reference point of the mesh.
574                    #
575                    #A point coordinate (x + x_0) is then made relative
576                    #to xi_0 by
577                    #
578                    # x_new = x + x_0 - xi_0
579                    #
580                    #and similarly for eta
581
582                    x_offset = data_origin[1] - mesh_origin[1]
583                    y_offset = data_origin[2] - mesh_origin[2]
584                else: #Shift back to a zero origin
585                    x_offset = data_origin[1]
586                    y_offset = data_origin[2]
587
588                point_coordinates[:,0] += x_offset
589                point_coordinates[:,1] += y_offset
590            else:
591                if mesh_origin is not None:
592                    #Use mesh origin for data points
593                    point_coordinates[:,0] -= mesh_origin[1]
594                    point_coordinates[:,1] -= mesh_origin[2]
595
596
597
598        #Remove points falling outside mesh boundary
599        #This reduced one example from 1356 seconds to 825 seconds
600
601       
602        if precrop is True:
603            from Numeric import take
604
605            if verbose: print 'Getting boundary polygon'
606            P = self.mesh.get_boundary_polygon()
607
608            if verbose: print 'Getting indices inside mesh boundary'
609            indices = inside_polygon(point_coordinates, P, verbose = verbose)
610
611
612            if len(indices) != point_coordinates.shape[0]:
613                self.cropped_points = True
614                if verbose:
615                    print 'Done - %d points outside mesh have been cropped.'\
616                          %(point_coordinates.shape[0] - len(indices))
617
618            point_coordinates = take(point_coordinates, indices)
619            self.point_indices = indices
620
621
622
623
624        #Build n x m interpolation matrix
625        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
626        n = point_coordinates.shape[0]     #Nbr of data points
627
628        if verbose: print 'Number of datapoints: %d' %n
629        if verbose: print 'Number of basis functions: %d' %m
630
631        #FIXME (Ole): We should use CSR here since mat-mat mult is now OK.
632        #However, Sparse_CSR does not have the same methods as Sparse yet
633        #The tests will reveal what needs to be done
634
635        #
636        #self.A = Sparse_CSR(Sparse(n,m))
637        #self.AtA = Sparse_CSR(Sparse(m,m))
638        self.A = Sparse(n,m)
639        self.AtA = Sparse(m,m)
640
641        #Build quad tree of vertices (FIXME: Is this the right spot for that?)
642        root = build_quadtree(self.mesh,
643                              max_points_per_cell = max_points_per_cell)
644        #root.show()
645        self.expanded_quad_searches = []
646        #Compute matrix elements
647        for i in range(n):
648            #For each data_coordinate point
649
650            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
651            x = point_coordinates[i]
652
653            #Find vertices near x
654            candidate_vertices = root.search(x[0], x[1])
655            is_more_elements = True
656
657            element_found, sigma0, sigma1, sigma2, k = \
658                self.search_triangles_of_vertices(candidate_vertices, x)
659            first_expansion = True
660            while not element_found and is_more_elements and expand_search:
661                #if verbose: print 'Expanding search'
662                if first_expansion == True:
663                    self.expanded_quad_searches.append(1)
664                    first_expansion = False
665                else:
666                    end = len(self.expanded_quad_searches) - 1
667                    assert end >= 0
668                    self.expanded_quad_searches[end] += 1
669                candidate_vertices, branch = root.expand_search()
670                if branch == []:
671                    # Searching all the verts from the root cell that haven't
672                    # been searched.  This is the last try
673                    element_found, sigma0, sigma1, sigma2, k = \
674                      self.search_triangles_of_vertices(candidate_vertices, x)
675                    is_more_elements = False
676                else:
677                    element_found, sigma0, sigma1, sigma2, k = \
678                      self.search_triangles_of_vertices(candidate_vertices, x)
679
680               
681            #Update interpolation matrix A if necessary
682            if element_found is True:
683                #Assign values to matrix A
684
685                j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
686                j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
687                j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
688
689                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
690                js     = [j0,j1,j2]
691
692                for j in js:
693                    self.A[i,j] = sigmas[j]
694                    for k in js:
695                        if interp_only == False:
696                            self.AtA[j,k] += sigmas[j]*sigmas[k]
697            else:
698                pass
699                #Ok if there is no triangle for datapoint
700                #(as in brute force version)
701                #raise 'Could not find triangle for point', x
702
703
704
705    def search_triangles_of_vertices(self, candidate_vertices, x):
706            #Find triangle containing x:
707            element_found = False
708
709            # This will be returned if element_found = False
710            sigma2 = -10.0
711            sigma0 = -10.0
712            sigma1 = -10.0
713            k = -10.0
714            #print "*$* candidate_vertices", candidate_vertices
715            #For all vertices in same cell as point x
716            for v in candidate_vertices:
717                #FIXME (DSG-DSG): this catches verts with no triangle.
718                #Currently pmesh is producing these.
719                #this should be stopped,
720                if self.mesh.vertexlist[v] is None:
721                    continue
722                #for each triangle id (k) which has v as a vertex
723                for k, _ in self.mesh.vertexlist[v]:
724
725                    #Get the three vertex_points of candidate triangle
726                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
727                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
728                    xi2 = self.mesh.get_vertex_coordinate(k, 2)
729
730                    #print "PDSG - k", k
731                    #print "PDSG - xi0", xi0
732                    #print "PDSG - xi1", xi1
733                    #print "PDSG - xi2", xi2
734                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
735                    #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
736
737                    #Get the three normals
738                    n0 = self.mesh.get_normal(k, 0)
739                    n1 = self.mesh.get_normal(k, 1)
740                    n2 = self.mesh.get_normal(k, 2)
741
742
743                    #Compute interpolation
744                    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
745                    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
746                    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
747
748                    #print "PDSG - sigma0", sigma0
749                    #print "PDSG - sigma1", sigma1
750                    #print "PDSG - sigma2", sigma2
751
752                    #FIXME: Maybe move out to test or something
753                    epsilon = 1.0e-6
754                    assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
755
756                    #Check that this triangle contains the data point
757
758                    #Sigmas can get negative within
759                    #machine precision on some machines (e.g nautilus)
760                    #Hence the small eps
761                    eps = 1.0e-15
762                    if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
763                        element_found = True
764                        break
765
766                if element_found is True:
767                    #Don't look for any other triangle
768                    break
769            return element_found, sigma0, sigma1, sigma2, k
770
771
772
773    def build_interpolation_matrix_A_brute(self, point_coordinates):
774        """Build n x m interpolation matrix, where
775        n is the number of data points and
776        m is the number of basis functions phi_k (one per vertex)
777
778        This is the brute force which is too slow for large problems,
779        but could be used for testing
780        """
781
782
783        #Convert input to Numeric arrays
784        point_coordinates = ensure_numeric(point_coordinates, Float)
785
786        #Build n x m interpolation matrix
787        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
788        n = point_coordinates.shape[0]     #Nbr of data points
789
790        self.A = Sparse(n,m)
791        self.AtA = Sparse(m,m)
792
793        #Compute matrix elements
794        for i in range(n):
795            #For each data_coordinate point
796
797            x = point_coordinates[i]
798            element_found = False
799            k = 0
800            while not element_found and k < len(self.mesh):
801                #For each triangle (brute force)
802                #FIXME: Real algorithm should only visit relevant triangles
803
804                #Get the three vertex_points
805                xi0 = self.mesh.get_vertex_coordinate(k, 0)
806                xi1 = self.mesh.get_vertex_coordinate(k, 1)
807                xi2 = self.mesh.get_vertex_coordinate(k, 2)
808
809                #Get the three normals
810                n0 = self.mesh.get_normal(k, 0)
811                n1 = self.mesh.get_normal(k, 1)
812                n2 = self.mesh.get_normal(k, 2)
813
814                #Compute interpolation
815                sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
816                sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
817                sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
818
819                #FIXME: Maybe move out to test or something
820                epsilon = 1.0e-6
821                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
822
823                #Check that this triangle contains data point
824                if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
825                    element_found = True
826                    #Assign values to matrix A
827
828                    j0 = self.mesh.triangles[k,0] #Global vertex id
829                    #self.A[i, j0] = sigma0
830
831                    j1 = self.mesh.triangles[k,1] #Global vertex id
832                    #self.A[i, j1] = sigma1
833
834                    j2 = self.mesh.triangles[k,2] #Global vertex id
835                    #self.A[i, j2] = sigma2
836
837                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
838                    js     = [j0,j1,j2]
839
840                    for j in js:
841                        self.A[i,j] = sigmas[j]
842                        for k in js:
843                            self.AtA[j,k] += sigmas[j]*sigmas[k]
844                k = k+1
845
846
847
848    def get_A(self):
849        return self.A.todense()
850
851    def get_B(self):
852        return self.B.todense()
853
854    def get_D(self):
855        return self.D.todense()
856
857        #FIXME: Remember to re-introduce the 1/n factor in the
858        #interpolation term
859
860    def build_smoothing_matrix_D(self):
861        """Build m x m smoothing matrix, where
862        m is the number of basis functions phi_k (one per vertex)
863
864        The smoothing matrix is defined as
865
866        D = D1 + D2
867
868        where
869
870        [D1]_{k,l} = \int_\Omega
871           \frac{\partial \phi_k}{\partial x}
872           \frac{\partial \phi_l}{\partial x}\,
873           dx dy
874
875        [D2]_{k,l} = \int_\Omega
876           \frac{\partial \phi_k}{\partial y}
877           \frac{\partial \phi_l}{\partial y}\,
878           dx dy
879
880
881        The derivatives \frac{\partial \phi_k}{\partial x},
882        \frac{\partial \phi_k}{\partial x} for a particular triangle
883        are obtained by computing the gradient a_k, b_k for basis function k
884        """
885
886        #FIXME: algorithm might be optimised by computing local 9x9
887        #"element stiffness matrices:
888
889        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
890
891        self.D = Sparse(m,m)
892
893        #For each triangle compute contributions to D = D1+D2
894        for i in range(len(self.mesh)):
895
896            #Get area
897            area = self.mesh.areas[i]
898
899            #Get global vertex indices
900            v0 = self.mesh.triangles[i,0]
901            v1 = self.mesh.triangles[i,1]
902            v2 = self.mesh.triangles[i,2]
903
904            #Get the three vertex_points
905            xi0 = self.mesh.get_vertex_coordinate(i, 0)
906            xi1 = self.mesh.get_vertex_coordinate(i, 1)
907            xi2 = self.mesh.get_vertex_coordinate(i, 2)
908
909            #Compute gradients for each vertex
910            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
911                              1, 0, 0)
912
913            a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
914                              0, 1, 0)
915
916            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
917                              0, 0, 1)
918
919            #Compute diagonal contributions
920            self.D[v0,v0] += (a0*a0 + b0*b0)*area
921            self.D[v1,v1] += (a1*a1 + b1*b1)*area
922            self.D[v2,v2] += (a2*a2 + b2*b2)*area
923
924            #Compute contributions for basis functions sharing edges
925            e01 = (a0*a1 + b0*b1)*area
926            self.D[v0,v1] += e01
927            self.D[v1,v0] += e01
928
929            e12 = (a1*a2 + b1*b2)*area
930            self.D[v1,v2] += e12
931            self.D[v2,v1] += e12
932
933            e20 = (a2*a0 + b2*b0)*area
934            self.D[v2,v0] += e20
935            self.D[v0,v2] += e20
936
937
938    def fit(self, z):
939        """Fit a smooth surface to given 1d array of data points z.
940
941        The smooth surface is computed at each vertex in the underlying
942        mesh using the formula given in the module doc string.
943
944        Pre Condition:
945          self.A, self.AtA and self.B have been initialised
946
947        Inputs:
948          z: Single 1d vector or array of data at the point_coordinates.
949        """
950
951        #Convert input to Numeric arrays
952        z = ensure_numeric(z, Float)
953
954        if len(z.shape) > 1 :
955            raise VectorShapeError, 'Can only deal with 1d data vector'
956
957        if self.point_indices is not None:
958            #Remove values for any points that were outside mesh
959            z = take(z, self.point_indices)
960
961        #Compute right hand side based on data
962        #FIXME (DSG-DsG): could Sparse_CSR be used here?  Use this format
963        # after a matrix is built, before calcs.
964        Atz = self.A.trans_mult(z)
965
966
967        #Check sanity
968        n, m = self.A.shape
969        if n<m and self.alpha == 0.0:
970            msg = 'ERROR (least_squares): Too few data points\n'
971            msg += 'There are only %d data points and alpha == 0. ' %n
972            msg += 'Need at least %d\n' %m
973            msg += 'Alternatively, set smoothing parameter alpha to a small '
974            msg += 'positive value,\ne.g. 1.0e-3.'
975            raise msg
976
977
978
979        return conjugate_gradient(self.B, Atz, Atz, imax=2*len(Atz) )
980        #FIXME: Should we store the result here for later use? (ON)
981
982
983    def fit_points(self, z, verbose=False):
984        """Like fit, but more robust when each point has two or more attributes
985        FIXME (Ole): The name fit_points doesn't carry any meaning
986        for me. How about something like fit_multiple or fit_columns?
987        """
988
989        try:
990            if verbose: print 'Solving penalised least_squares problem'
991            return self.fit(z)
992        except VectorShapeError, e:
993            # broadcasting is not supported.
994
995            #Convert input to Numeric arrays
996            z = ensure_numeric(z, Float)
997
998            #Build n x m interpolation matrix
999            m = self.mesh.coordinates.shape[0] #Number of vertices
1000            n = z.shape[1]                     #Number of data points
1001
1002            f = zeros((m,n), Float) #Resulting columns
1003
1004            for i in range(z.shape[1]):
1005                f[:,i] = self.fit(z[:,i])
1006
1007            return f
1008
1009
1010    def interpolate(self, f):
1011        """Evaluate smooth surface f at data points implied in self.A.
1012
1013        The mesh values representing a smooth surface are
1014        assumed to be specified in f. This argument could,
1015        for example have been obtained from the method self.fit()
1016
1017        Pre Condition:
1018          self.A has been initialised
1019
1020        Inputs:
1021          f: Vector or array of data at the mesh vertices.
1022          If f is an array, interpolation will be done for each column as
1023          per underlying matrix-matrix multiplication
1024
1025        Output:
1026          Interpolated values at data points implied in self.A
1027
1028        """
1029
1030        return self.A * f
1031
1032    def cull_outsiders(self, f):
1033        pass
1034
1035
1036
1037
1038class Interpolation_function:
1039    """Interpolation_function - creates callable object f(t, id) or f(t,x,y)
1040    which is interpolated from time series defined at vertices of
1041    triangular mesh (such as those stored in sww files)
1042
1043    Let m be the number of vertices, n the number of triangles
1044    and p the number of timesteps.
1045
1046    Mandatory input
1047        time:               px1 array of monotonously increasing times (Float)
1048        quantities:         Dictionary of arrays or 1 array (Float)
1049                            The arrays must either have dimensions pxm or mx1.
1050                            The resulting function will be time dependent in
1051                            the former case while it will be constant with
1052                            respect to time in the latter case.
1053       
1054    Optional input:
1055        quantity_names:     List of keys into the quantities dictionary
1056        vertex_coordinates: mx2 array of coordinates (Float)
1057        triangles:          nx3 array of indices into vertex_coordinates (Int)
1058        interpolation_points: Nx2 array of coordinates to be interpolated to
1059        verbose:            Level of reporting
1060   
1061   
1062    The quantities returned by the callable object are specified by
1063    the list quantities which must contain the names of the
1064    quantities to be returned and also reflect the order, e.g. for
1065    the shallow water wave equation, on would have
1066    quantities = ['stage', 'xmomentum', 'ymomentum']
1067
1068    The parameter interpolation_points decides at which points interpolated
1069    quantities are to be computed whenever object is called.
1070    If None, return average value
1071    """
1072
1073   
1074   
1075    def __init__(self,
1076                 time,
1077                 quantities,
1078                 quantity_names = None, 
1079                 vertex_coordinates = None,
1080                 triangles = None,
1081                 interpolation_points = None,
1082                 verbose = False):
1083        """Initialise object and build spatial interpolation if required
1084        """
1085
1086        from Numeric import array, zeros, Float, alltrue, concatenate,\
1087             reshape, ArrayType
1088
1089
1090        from config import time_format
1091        import types
1092
1093
1094
1095        #Check temporal info
1096        time = ensure_numeric(time)       
1097        msg = 'Time must be a monotonuosly '
1098        msg += 'increasing sequence %s' %time
1099        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
1100
1101
1102        #Check if quantities is a single array only
1103        if type(quantities) != types.DictType:
1104            quantities = ensure_numeric(quantities)
1105            quantity_names = ['Attribute']
1106
1107            #Make it a dictionary
1108            quantities = {quantity_names[0]: quantities}
1109
1110
1111        #Use keys if no names are specified
1112        if quantity_names is None:
1113            quantity_names = quantities.keys()
1114
1115
1116        #Check spatial info
1117        if vertex_coordinates is None:
1118            self.spatial = False
1119        else:   
1120            vertex_coordinates = ensure_numeric(vertex_coordinates)
1121
1122            assert triangles is not None, 'Triangles array must be specified'
1123            triangles = ensure_numeric(triangles)
1124            self.spatial = True           
1125           
1126
1127 
1128        #Save for use with statistics
1129        self.quantity_names = quantity_names       
1130        self.quantities = quantities       
1131        self.vertex_coordinates = vertex_coordinates
1132        self.interpolation_points = interpolation_points
1133        self.T = time[:]  # Time assumed to be relative to starttime
1134        self.index = 0    # Initial time index
1135        self.precomputed_values = {}
1136           
1137
1138
1139        #Precomputed spatial interpolation if requested
1140        if interpolation_points is not None:
1141            if self.spatial is False:
1142                raise 'Triangles and vertex_coordinates must be specified'
1143           
1144            try:
1145                self.interpolation_points = ensure_numeric(interpolation_points)
1146            except:
1147                msg = 'Interpolation points must be an N x 2 Numeric array '+\
1148                      'or a list of points\n'
1149                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
1150                                      '...')
1151                raise msg
1152
1153
1154            m = len(self.interpolation_points)
1155            p = len(self.T)
1156           
1157            for name in quantity_names:
1158                self.precomputed_values[name] = zeros((p, m), Float)
1159
1160            #Build interpolator
1161            interpol = Interpolation(vertex_coordinates,
1162                                     triangles,
1163                                     point_coordinates = \
1164                                     self.interpolation_points,
1165                                     alpha = 0,
1166                                     precrop = False, 
1167                                     verbose = verbose)
1168
1169            if verbose: print 'Interpolate'
1170            for i, t in enumerate(self.T):
1171                #Interpolate quantities at this timestep
1172                if verbose and i%((p+10)/10)==0:
1173                    print ' time step %d of %d' %(i, p)
1174                   
1175                for name in quantity_names:
1176                    if len(quantities[name].shape) == 2:
1177                        result = interpol.interpolate(quantities[name][i,:])
1178                    else:
1179                       #Assume no time dependency
1180                       result = interpol.interpolate(quantities[name][:])
1181                       
1182                    self.precomputed_values[name][i, :] = result
1183                   
1184                       
1185
1186            #Report
1187            if verbose:
1188                print self.statistics()
1189                #self.print_statistics()
1190           
1191        else:
1192            #Store quantitites as is
1193            for name in quantity_names:
1194                self.precomputed_values[name] = quantities[name]
1195
1196
1197        #else:
1198        #    #Return an average, making this a time series
1199        #    for name in quantity_names:
1200        #        self.values[name] = zeros(len(self.T), Float)
1201        #
1202        #    if verbose: print 'Compute mean values'
1203        #    for i, t in enumerate(self.T):
1204        #        if verbose: print ' time step %d of %d' %(i, len(self.T))
1205        #        for name in quantity_names:
1206        #           self.values[name][i] = mean(quantities[name][i,:])
1207
1208
1209
1210
1211    def __repr__(self):
1212        #return 'Interpolation function (spatio-temporal)'
1213        return self.statistics()
1214   
1215
1216    def __call__(self, t, point_id = None, x = None, y = None):
1217        """Evaluate f(t), f(t, point_id) or f(t, x, y)
1218
1219        Inputs:
1220          t: time - Model time. Must lie within existing timesteps
1221          point_id: index of one of the preprocessed points.
1222          x, y:     Overrides location, point_id ignored
1223         
1224          If spatial info is present and all of x,y,point_id
1225          are None an exception is raised
1226                   
1227          If no spatial info is present, point_id and x,y arguments are ignored
1228          making f a function of time only.
1229
1230         
1231          FIXME: point_id could also be a slice
1232          FIXME: What if x and y are vectors?
1233          FIXME: What about f(x,y) without t?
1234        """
1235
1236        from math import pi, cos, sin, sqrt
1237        from Numeric import zeros, Float
1238        from utilities.numerical_tools import mean       
1239
1240        if self.spatial is True:
1241            if point_id is None:
1242                if x is None or y is None:
1243                    msg = 'Either point_id or x and y must be specified'
1244                    raise msg
1245            else:
1246                if self.interpolation_points is None:
1247                    msg = 'Interpolation_function must be instantiated ' +\
1248                          'with a list of interpolation points before parameter ' +\
1249                          'point_id can be used'
1250                    raise msg
1251
1252
1253        msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
1254        msg += ' does not match model time: %s\n' %t
1255        if t < self.T[0]: raise msg
1256        if t > self.T[-1]: raise msg
1257
1258        oldindex = self.index #Time index
1259
1260        #Find current time slot
1261        while t > self.T[self.index]: self.index += 1
1262        while t < self.T[self.index]: self.index -= 1
1263
1264        if t == self.T[self.index]:
1265            #Protect against case where t == T[-1] (last time)
1266            # - also works in general when t == T[i]
1267            ratio = 0
1268        else:
1269            #t is now between index and index+1
1270            ratio = (t - self.T[self.index])/\
1271                    (self.T[self.index+1] - self.T[self.index])
1272
1273        #Compute interpolated values
1274        q = zeros(len(self.quantity_names), Float)
1275
1276        for i, name in enumerate(self.quantity_names):
1277            Q = self.precomputed_values[name]
1278
1279            if self.spatial is False:
1280                #If there is no spatial info               
1281                assert len(Q.shape) == 1
1282
1283                Q0 = Q[self.index]
1284                if ratio > 0: Q1 = Q[self.index+1]
1285
1286            else:
1287                if x is not None and y is not None:
1288                    #Interpolate to x, y
1289                   
1290                    raise 'x,y interpolation not yet implemented'
1291                else:
1292                    #Use precomputed point
1293                    Q0 = Q[self.index, point_id]
1294                    if ratio > 0: Q1 = Q[self.index+1, point_id]
1295
1296            #Linear temporal interpolation   
1297            if ratio > 0:
1298                q[i] = Q0 + ratio*(Q1 - Q0)
1299            else:
1300                q[i] = Q0
1301
1302
1303        #Return vector of interpolated values
1304        #if len(q) == 1:
1305        #    return q[0]
1306        #else:
1307        #    return q
1308
1309
1310        #Return vector of interpolated values
1311        #FIXME:
1312        if self.spatial is True:
1313            return q
1314        else:
1315            #Replicate q according to x and y
1316            #This is e.g used for Wind_stress
1317            if x is None or y is None: 
1318                return q
1319            else:
1320                try:
1321                    N = len(x)
1322                except:
1323                    return q
1324                else:
1325                    from Numeric import ones, Float
1326                    #x is a vector - Create one constant column for each value
1327                    N = len(x)
1328                    assert len(y) == N, 'x and y must have same length'
1329                    res = []
1330                    for col in q:
1331                        res.append(col*ones(N, Float))
1332                       
1333                return res
1334
1335
1336    def statistics(self):
1337        """Output statistics about interpolation_function
1338        """
1339       
1340        vertex_coordinates = self.vertex_coordinates
1341        interpolation_points = self.interpolation_points               
1342        quantity_names = self.quantity_names
1343        quantities = self.quantities
1344        precomputed_values = self.precomputed_values                 
1345               
1346        x = vertex_coordinates[:,0]
1347        y = vertex_coordinates[:,1]               
1348
1349        str =  '------------------------------------------------\n'
1350        str += 'Interpolation_function (spatio-temporal) statistics:\n'
1351        str += '  Extent:\n'
1352        str += '    x in [%f, %f], len(x) == %d\n'\
1353               %(min(x), max(x), len(x))
1354        str += '    y in [%f, %f], len(y) == %d\n'\
1355               %(min(y), max(y), len(y))
1356        str += '    t in [%f, %f], len(t) == %d\n'\
1357               %(min(self.T), max(self.T), len(self.T))
1358        str += '  Quantities:\n'
1359        for name in quantity_names:
1360            q = quantities[name][:].flat
1361            str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
1362
1363        if interpolation_points is not None:   
1364            str += '  Interpolation points (xi, eta):'\
1365                   ' number of points == %d\n' %interpolation_points.shape[0]
1366            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1367                                            max(interpolation_points[:,0]))
1368            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1369                                             max(interpolation_points[:,1]))
1370            str += '  Interpolated quantities (over all timesteps):\n'
1371       
1372            for name in quantity_names:
1373                q = precomputed_values[name][:].flat
1374                str += '    %s at interpolation points in [%f, %f]\n'\
1375                       %(name, min(q), max(q))
1376        str += '------------------------------------------------\n'
1377
1378        return str
1379
1380        #FIXME: Delete
1381        #print '------------------------------------------------'
1382        #print 'Interpolation_function statistics:'
1383        #print '  Extent:'
1384        #print '    x in [%f, %f], len(x) == %d'\
1385        #      %(min(x), max(x), len(x))
1386        #print '    y in [%f, %f], len(y) == %d'\
1387        #      %(min(y), max(y), len(y))
1388        #print '    t in [%f, %f], len(t) == %d'\
1389        #      %(min(self.T), max(self.T), len(self.T))
1390        #print '  Quantities:'
1391        #for name in quantity_names:
1392        #    q = quantities[name][:].flat
1393        #    print '    %s in [%f, %f]' %(name, min(q), max(q))
1394        #print '  Interpolation points (xi, eta):'\
1395        #      ' number of points == %d ' %interpolation_points.shape[0]
1396        #print '    xi in [%f, %f]' %(min(interpolation_points[:,0]),
1397        #                             max(interpolation_points[:,0]))
1398        #print '    eta in [%f, %f]' %(min(interpolation_points[:,1]),
1399        #                              max(interpolation_points[:,1]))
1400        #print '  Interpolated quantities (over all timesteps):'
1401        #
1402        #for name in quantity_names:
1403        #    q = precomputed_values[name][:].flat
1404        #    print '    %s at interpolation points in [%f, %f]'\
1405        #          %(name, min(q), max(q))
1406        #print '------------------------------------------------'
1407
1408
1409#-------------------------------------------------------------
1410if __name__ == "__main__":
1411    """
1412    Load in a mesh and data points with attributes.
1413    Fit the attributes to the mesh.
1414    Save a new mesh file.
1415    """
1416    import os, sys
1417    usage = "usage: %s mesh_input.tsh point.xya mesh_output.tsh [expand|no_expand][vervose|non_verbose] [alpha] [display_errors|no_display_errors]"\
1418            %os.path.basename(sys.argv[0])
1419
1420    if len(sys.argv) < 4:
1421        print usage
1422    else:
1423        mesh_file = sys.argv[1]
1424        point_file = sys.argv[2]
1425        mesh_output_file = sys.argv[3]
1426
1427        expand_search = False
1428        if len(sys.argv) > 4:
1429            if sys.argv[4][0] == "e" or sys.argv[4][0] == "E":
1430                expand_search = True
1431            else:
1432                expand_search = False
1433
1434        verbose = False
1435        if len(sys.argv) > 5:
1436            if sys.argv[5][0] == "n" or sys.argv[5][0] == "N":
1437                verbose = False
1438            else:
1439                verbose = True
1440
1441        if len(sys.argv) > 6:
1442            alpha = sys.argv[6]
1443        else:
1444            alpha = DEFAULT_ALPHA
1445
1446        # This is used more for testing
1447        if len(sys.argv) > 7:
1448            if sys.argv[7][0] == "n" or sys.argv[5][0] == "N":
1449                display_errors = False
1450            else:
1451                display_errors = True
1452           
1453        t0 = time.time()
1454        try:
1455            fit_to_mesh_file(mesh_file,
1456                         point_file,
1457                         mesh_output_file,
1458                         alpha,
1459                         verbose= verbose,
1460                         expand_search = expand_search,
1461                         display_errors = display_errors)
1462        except IOError,e:
1463            import sys; sys.exit(1)
1464
1465        print 'That took %.2f seconds' %(time.time()-t0)
1466
Note: See TracBrowser for help on using the repository browser.