source: inundation/pyvolution/least_squares.py @ 2682

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

updates to onslow
add print outs to least_squares and mesh
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        print 'mesh.check_integrity'
444        self.mesh.check_integrity()
445        print 'data_origin'
446        self.data_origin = data_origin
447
448        self.point_indices = None
449
450        #Smoothing parameter
451        if alpha is None:
452            self.alpha = DEFAULT_ALPHA
453        else:   
454            self.alpha = alpha
455
456
457        if point_coordinates is not None:
458            if verbose: print 'Building interpolation matrix'
459            self.build_interpolation_matrix_A(point_coordinates,
460                                              verbose = verbose,
461                                              expand_search = expand_search,
462                                              interp_only = interp_only, 
463                                              max_points_per_cell =\
464                                              max_points_per_cell,
465                                              data_origin = data_origin,
466                                              precrop = precrop)
467        #Build coefficient matrices
468        if interp_only == False:
469            self.build_coefficient_matrix_B(point_coordinates,
470                                        verbose = verbose,
471                                        expand_search = expand_search,
472                                        max_points_per_cell =\
473                                        max_points_per_cell,
474                                        data_origin = data_origin,
475                                        precrop = precrop)
476        print 'finished interpolation'
477
478    def set_point_coordinates(self, point_coordinates,
479                              data_origin = None,
480                              verbose = False,
481                              precrop = True):
482        """
483        A public interface to setting the point co-ordinates.
484        """
485        if point_coordinates is not None:
486            if verbose: print 'Building interpolation matrix'
487            self.build_interpolation_matrix_A(point_coordinates,
488                                              verbose = verbose,
489                                              data_origin = data_origin,
490                                              precrop = precrop)
491        self.build_coefficient_matrix_B(point_coordinates, data_origin)
492
493    def build_coefficient_matrix_B(self, point_coordinates=None,
494                                   verbose = False, expand_search = True,
495                                   max_points_per_cell=30,
496                                   data_origin = None,
497                                   precrop = False):
498        """Build final coefficient matrix"""
499
500
501        if self.alpha <> 0:
502            if verbose: print 'Building smoothing matrix'
503            self.build_smoothing_matrix_D()
504
505        if point_coordinates is not None:
506            if self.alpha <> 0:
507                self.B = self.AtA + self.alpha*self.D
508            else:
509                self.B = self.AtA
510
511            #Convert self.B matrix to CSR format for faster matrix vector
512            self.B = Sparse_CSR(self.B)
513
514    def build_interpolation_matrix_A(self, point_coordinates,
515                                     verbose = False, expand_search = True,
516                                     max_points_per_cell=30,
517                                     data_origin = None,
518                                     precrop = False,
519                                     interp_only = False):
520        """Build n x m interpolation matrix, where
521        n is the number of data points and
522        m is the number of basis functions phi_k (one per vertex)
523
524        This algorithm uses a quad tree data structure for fast binning of data points
525        origin is a 3-tuple consisting of UTM zone, easting and northing.
526        If specified coordinates are assumed to be relative to this origin.
527
528        This one will override any data_origin that may be specified in
529        interpolation instance
530
531        """
532
533
534
535        #FIXME (Ole): Check that this function is memeory efficient.
536        #6 million datapoints and 300000 basis functions
537        #causes out-of-memory situation
538        #First thing to check is whether there is room for self.A and self.AtA
539        #
540        #Maybe we need some sort of blocking
541
542        from pyvolution.quad import build_quadtree
543        from utilities.polygon import inside_polygon
544       
545
546        if data_origin is None:
547            data_origin = self.data_origin #Use the one from
548                                           #interpolation instance
549
550        #Convert input to Numeric arrays just in case.
551        point_coordinates = ensure_numeric(point_coordinates, Float)
552
553        #Keep track of discarded points (if any).
554        #This is only registered if precrop is True
555        self.cropped_points = False
556
557        #Shift data points to same origin as mesh (if specified)
558
559        #FIXME this will shift if there was no geo_ref.
560        #But all this should be removed anyhow.
561        #change coords before this point
562        mesh_origin = self.mesh.geo_reference.get_origin()
563        if point_coordinates is not None:
564            if data_origin is not None:
565                if mesh_origin is not None:
566
567                    #Transformation:
568                    #
569                    #Let x_0 be the reference point of the point coordinates
570                    #and xi_0 the reference point of the mesh.
571                    #
572                    #A point coordinate (x + x_0) is then made relative
573                    #to xi_0 by
574                    #
575                    # x_new = x + x_0 - xi_0
576                    #
577                    #and similarly for eta
578
579                    x_offset = data_origin[1] - mesh_origin[1]
580                    y_offset = data_origin[2] - mesh_origin[2]
581                else: #Shift back to a zero origin
582                    x_offset = data_origin[1]
583                    y_offset = data_origin[2]
584
585                point_coordinates[:,0] += x_offset
586                point_coordinates[:,1] += y_offset
587            else:
588                if mesh_origin is not None:
589                    #Use mesh origin for data points
590                    point_coordinates[:,0] -= mesh_origin[1]
591                    point_coordinates[:,1] -= mesh_origin[2]
592
593
594
595        #Remove points falling outside mesh boundary
596        #This reduced one example from 1356 seconds to 825 seconds
597
598       
599        if precrop is True:
600            from Numeric import take
601
602            if verbose: print 'Getting boundary polygon'
603            P = self.mesh.get_boundary_polygon()
604
605            if verbose: print 'Getting indices inside mesh boundary'
606            indices = inside_polygon(point_coordinates, P, verbose = verbose)
607
608
609            if len(indices) != point_coordinates.shape[0]:
610                self.cropped_points = True
611                if verbose:
612                    print 'Done - %d points outside mesh have been cropped.'\
613                          %(point_coordinates.shape[0] - len(indices))
614
615            point_coordinates = take(point_coordinates, indices)
616            self.point_indices = indices
617
618
619
620
621        #Build n x m interpolation matrix
622        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
623        n = point_coordinates.shape[0]     #Nbr of data points
624
625        if verbose: print 'Number of datapoints: %d' %n
626        if verbose: print 'Number of basis functions: %d' %m
627
628        #FIXME (Ole): We should use CSR here since mat-mat mult is now OK.
629        #However, Sparse_CSR does not have the same methods as Sparse yet
630        #The tests will reveal what needs to be done
631
632        #
633        #self.A = Sparse_CSR(Sparse(n,m))
634        #self.AtA = Sparse_CSR(Sparse(m,m))
635        self.A = Sparse(n,m)
636        self.AtA = Sparse(m,m)
637
638        #Build quad tree of vertices (FIXME: Is this the right spot for that?)
639        root = build_quadtree(self.mesh,
640                              max_points_per_cell = max_points_per_cell)
641        #root.show()
642        self.expanded_quad_searches = []
643        #Compute matrix elements
644        for i in range(n):
645            #For each data_coordinate point
646
647            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
648            x = point_coordinates[i]
649
650            #Find vertices near x
651            candidate_vertices = root.search(x[0], x[1])
652            is_more_elements = True
653
654            element_found, sigma0, sigma1, sigma2, k = \
655                self.search_triangles_of_vertices(candidate_vertices, x)
656            first_expansion = True
657            while not element_found and is_more_elements and expand_search:
658                #if verbose: print 'Expanding search'
659                if first_expansion == True:
660                    self.expanded_quad_searches.append(1)
661                    first_expansion = False
662                else:
663                    end = len(self.expanded_quad_searches) - 1
664                    assert end >= 0
665                    self.expanded_quad_searches[end] += 1
666                candidate_vertices, branch = root.expand_search()
667                if branch == []:
668                    # Searching all the verts from the root cell that haven't
669                    # been searched.  This is the last try
670                    element_found, sigma0, sigma1, sigma2, k = \
671                      self.search_triangles_of_vertices(candidate_vertices, x)
672                    is_more_elements = False
673                else:
674                    element_found, sigma0, sigma1, sigma2, k = \
675                      self.search_triangles_of_vertices(candidate_vertices, x)
676
677               
678            #Update interpolation matrix A if necessary
679            if element_found is True:
680                #Assign values to matrix A
681
682                j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
683                j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
684                j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
685
686                sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
687                js     = [j0,j1,j2]
688
689                for j in js:
690                    self.A[i,j] = sigmas[j]
691                    for k in js:
692                        if interp_only == False:
693                            self.AtA[j,k] += sigmas[j]*sigmas[k]
694            else:
695                pass
696                #Ok if there is no triangle for datapoint
697                #(as in brute force version)
698                #raise 'Could not find triangle for point', x
699
700
701
702    def search_triangles_of_vertices(self, candidate_vertices, x):
703            #Find triangle containing x:
704            element_found = False
705
706            # This will be returned if element_found = False
707            sigma2 = -10.0
708            sigma0 = -10.0
709            sigma1 = -10.0
710            k = -10.0
711            #print "*$* candidate_vertices", candidate_vertices
712            #For all vertices in same cell as point x
713            for v in candidate_vertices:
714                #FIXME (DSG-DSG): this catches verts with no triangle.
715                #Currently pmesh is producing these.
716                #this should be stopped,
717                if self.mesh.vertexlist[v] is None:
718                    continue
719                #for each triangle id (k) which has v as a vertex
720                for k, _ in self.mesh.vertexlist[v]:
721
722                    #Get the three vertex_points of candidate triangle
723                    xi0 = self.mesh.get_vertex_coordinate(k, 0)
724                    xi1 = self.mesh.get_vertex_coordinate(k, 1)
725                    xi2 = self.mesh.get_vertex_coordinate(k, 2)
726
727                    #print "PDSG - k", k
728                    #print "PDSG - xi0", xi0
729                    #print "PDSG - xi1", xi1
730                    #print "PDSG - xi2", xi2
731                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
732                    #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
733
734                    #Get the three normals
735                    n0 = self.mesh.get_normal(k, 0)
736                    n1 = self.mesh.get_normal(k, 1)
737                    n2 = self.mesh.get_normal(k, 2)
738
739
740                    #Compute interpolation
741                    sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
742                    sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
743                    sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
744
745                    #print "PDSG - sigma0", sigma0
746                    #print "PDSG - sigma1", sigma1
747                    #print "PDSG - sigma2", sigma2
748
749                    #FIXME: Maybe move out to test or something
750                    epsilon = 1.0e-6
751                    assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
752
753                    #Check that this triangle contains the data point
754
755                    #Sigmas can get negative within
756                    #machine precision on some machines (e.g nautilus)
757                    #Hence the small eps
758                    eps = 1.0e-15
759                    if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
760                        element_found = True
761                        break
762
763                if element_found is True:
764                    #Don't look for any other triangle
765                    break
766            return element_found, sigma0, sigma1, sigma2, k
767
768
769
770    def build_interpolation_matrix_A_brute(self, point_coordinates):
771        """Build n x m interpolation matrix, where
772        n is the number of data points and
773        m is the number of basis functions phi_k (one per vertex)
774
775        This is the brute force which is too slow for large problems,
776        but could be used for testing
777        """
778
779
780        #Convert input to Numeric arrays
781        point_coordinates = ensure_numeric(point_coordinates, Float)
782
783        #Build n x m interpolation matrix
784        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
785        n = point_coordinates.shape[0]     #Nbr of data points
786
787        self.A = Sparse(n,m)
788        self.AtA = Sparse(m,m)
789
790        #Compute matrix elements
791        for i in range(n):
792            #For each data_coordinate point
793
794            x = point_coordinates[i]
795            element_found = False
796            k = 0
797            while not element_found and k < len(self.mesh):
798                #For each triangle (brute force)
799                #FIXME: Real algorithm should only visit relevant triangles
800
801                #Get the three vertex_points
802                xi0 = self.mesh.get_vertex_coordinate(k, 0)
803                xi1 = self.mesh.get_vertex_coordinate(k, 1)
804                xi2 = self.mesh.get_vertex_coordinate(k, 2)
805
806                #Get the three normals
807                n0 = self.mesh.get_normal(k, 0)
808                n1 = self.mesh.get_normal(k, 1)
809                n2 = self.mesh.get_normal(k, 2)
810
811                #Compute interpolation
812                sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
813                sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
814                sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
815
816                #FIXME: Maybe move out to test or something
817                epsilon = 1.0e-6
818                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
819
820                #Check that this triangle contains data point
821                if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
822                    element_found = True
823                    #Assign values to matrix A
824
825                    j0 = self.mesh.triangles[k,0] #Global vertex id
826                    #self.A[i, j0] = sigma0
827
828                    j1 = self.mesh.triangles[k,1] #Global vertex id
829                    #self.A[i, j1] = sigma1
830
831                    j2 = self.mesh.triangles[k,2] #Global vertex id
832                    #self.A[i, j2] = sigma2
833
834                    sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
835                    js     = [j0,j1,j2]
836
837                    for j in js:
838                        self.A[i,j] = sigmas[j]
839                        for k in js:
840                            self.AtA[j,k] += sigmas[j]*sigmas[k]
841                k = k+1
842
843
844
845    def get_A(self):
846        return self.A.todense()
847
848    def get_B(self):
849        return self.B.todense()
850
851    def get_D(self):
852        return self.D.todense()
853
854        #FIXME: Remember to re-introduce the 1/n factor in the
855        #interpolation term
856
857    def build_smoothing_matrix_D(self):
858        """Build m x m smoothing matrix, where
859        m is the number of basis functions phi_k (one per vertex)
860
861        The smoothing matrix is defined as
862
863        D = D1 + D2
864
865        where
866
867        [D1]_{k,l} = \int_\Omega
868           \frac{\partial \phi_k}{\partial x}
869           \frac{\partial \phi_l}{\partial x}\,
870           dx dy
871
872        [D2]_{k,l} = \int_\Omega
873           \frac{\partial \phi_k}{\partial y}
874           \frac{\partial \phi_l}{\partial y}\,
875           dx dy
876
877
878        The derivatives \frac{\partial \phi_k}{\partial x},
879        \frac{\partial \phi_k}{\partial x} for a particular triangle
880        are obtained by computing the gradient a_k, b_k for basis function k
881        """
882
883        #FIXME: algorithm might be optimised by computing local 9x9
884        #"element stiffness matrices:
885
886        m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
887
888        self.D = Sparse(m,m)
889
890        #For each triangle compute contributions to D = D1+D2
891        for i in range(len(self.mesh)):
892
893            #Get area
894            area = self.mesh.areas[i]
895
896            #Get global vertex indices
897            v0 = self.mesh.triangles[i,0]
898            v1 = self.mesh.triangles[i,1]
899            v2 = self.mesh.triangles[i,2]
900
901            #Get the three vertex_points
902            xi0 = self.mesh.get_vertex_coordinate(i, 0)
903            xi1 = self.mesh.get_vertex_coordinate(i, 1)
904            xi2 = self.mesh.get_vertex_coordinate(i, 2)
905
906            #Compute gradients for each vertex
907            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
908                              1, 0, 0)
909
910            a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
911                              0, 1, 0)
912
913            a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
914                              0, 0, 1)
915
916            #Compute diagonal contributions
917            self.D[v0,v0] += (a0*a0 + b0*b0)*area
918            self.D[v1,v1] += (a1*a1 + b1*b1)*area
919            self.D[v2,v2] += (a2*a2 + b2*b2)*area
920
921            #Compute contributions for basis functions sharing edges
922            e01 = (a0*a1 + b0*b1)*area
923            self.D[v0,v1] += e01
924            self.D[v1,v0] += e01
925
926            e12 = (a1*a2 + b1*b2)*area
927            self.D[v1,v2] += e12
928            self.D[v2,v1] += e12
929
930            e20 = (a2*a0 + b2*b0)*area
931            self.D[v2,v0] += e20
932            self.D[v0,v2] += e20
933
934
935    def fit(self, z):
936        """Fit a smooth surface to given 1d array of data points z.
937
938        The smooth surface is computed at each vertex in the underlying
939        mesh using the formula given in the module doc string.
940
941        Pre Condition:
942          self.A, self.AtA and self.B have been initialised
943
944        Inputs:
945          z: Single 1d vector or array of data at the point_coordinates.
946        """
947
948        #Convert input to Numeric arrays
949        z = ensure_numeric(z, Float)
950
951        if len(z.shape) > 1 :
952            raise VectorShapeError, 'Can only deal with 1d data vector'
953
954        if self.point_indices is not None:
955            #Remove values for any points that were outside mesh
956            z = take(z, self.point_indices)
957
958        #Compute right hand side based on data
959        #FIXME (DSG-DsG): could Sparse_CSR be used here?  Use this format
960        # after a matrix is built, before calcs.
961        Atz = self.A.trans_mult(z)
962
963
964        #Check sanity
965        n, m = self.A.shape
966        if n<m and self.alpha == 0.0:
967            msg = 'ERROR (least_squares): Too few data points\n'
968            msg += 'There are only %d data points and alpha == 0. ' %n
969            msg += 'Need at least %d\n' %m
970            msg += 'Alternatively, set smoothing parameter alpha to a small '
971            msg += 'positive value,\ne.g. 1.0e-3.'
972            raise msg
973
974
975
976        return conjugate_gradient(self.B, Atz, Atz, imax=2*len(Atz) )
977        #FIXME: Should we store the result here for later use? (ON)
978
979
980    def fit_points(self, z, verbose=False):
981        """Like fit, but more robust when each point has two or more attributes
982        FIXME (Ole): The name fit_points doesn't carry any meaning
983        for me. How about something like fit_multiple or fit_columns?
984        """
985
986        try:
987            if verbose: print 'Solving penalised least_squares problem'
988            return self.fit(z)
989        except VectorShapeError, e:
990            # broadcasting is not supported.
991
992            #Convert input to Numeric arrays
993            z = ensure_numeric(z, Float)
994
995            #Build n x m interpolation matrix
996            m = self.mesh.coordinates.shape[0] #Number of vertices
997            n = z.shape[1]                     #Number of data points
998
999            f = zeros((m,n), Float) #Resulting columns
1000
1001            for i in range(z.shape[1]):
1002                f[:,i] = self.fit(z[:,i])
1003
1004            return f
1005
1006
1007    def interpolate(self, f):
1008        """Evaluate smooth surface f at data points implied in self.A.
1009
1010        The mesh values representing a smooth surface are
1011        assumed to be specified in f. This argument could,
1012        for example have been obtained from the method self.fit()
1013
1014        Pre Condition:
1015          self.A has been initialised
1016
1017        Inputs:
1018          f: Vector or array of data at the mesh vertices.
1019          If f is an array, interpolation will be done for each column as
1020          per underlying matrix-matrix multiplication
1021
1022        Output:
1023          Interpolated values at data points implied in self.A
1024
1025        """
1026
1027        return self.A * f
1028
1029    def cull_outsiders(self, f):
1030        pass
1031
1032
1033
1034
1035class Interpolation_function:
1036    """Interpolation_function - creates callable object f(t, id) or f(t,x,y)
1037    which is interpolated from time series defined at vertices of
1038    triangular mesh (such as those stored in sww files)
1039
1040    Let m be the number of vertices, n the number of triangles
1041    and p the number of timesteps.
1042
1043    Mandatory input
1044        time:               px1 array of monotonously increasing times (Float)
1045        quantities:         Dictionary of arrays or 1 array (Float)
1046                            The arrays must either have dimensions pxm or mx1.
1047                            The resulting function will be time dependent in
1048                            the former case while it will be constant with
1049                            respect to time in the latter case.
1050       
1051    Optional input:
1052        quantity_names:     List of keys into the quantities dictionary
1053        vertex_coordinates: mx2 array of coordinates (Float)
1054        triangles:          nx3 array of indices into vertex_coordinates (Int)
1055        interpolation_points: Nx2 array of coordinates to be interpolated to
1056        verbose:            Level of reporting
1057   
1058   
1059    The quantities returned by the callable object are specified by
1060    the list quantities which must contain the names of the
1061    quantities to be returned and also reflect the order, e.g. for
1062    the shallow water wave equation, on would have
1063    quantities = ['stage', 'xmomentum', 'ymomentum']
1064
1065    The parameter interpolation_points decides at which points interpolated
1066    quantities are to be computed whenever object is called.
1067    If None, return average value
1068    """
1069
1070   
1071   
1072    def __init__(self,
1073                 time,
1074                 quantities,
1075                 quantity_names = None, 
1076                 vertex_coordinates = None,
1077                 triangles = None,
1078                 interpolation_points = None,
1079                 verbose = False):
1080        """Initialise object and build spatial interpolation if required
1081        """
1082
1083        from Numeric import array, zeros, Float, alltrue, concatenate,\
1084             reshape, ArrayType
1085
1086
1087        from config import time_format
1088        import types
1089
1090
1091
1092        #Check temporal info
1093        time = ensure_numeric(time)       
1094        msg = 'Time must be a monotonuosly '
1095        msg += 'increasing sequence %s' %time
1096        assert alltrue(time[1:] - time[:-1] >= 0 ), msg
1097
1098
1099        #Check if quantities is a single array only
1100        if type(quantities) != types.DictType:
1101            quantities = ensure_numeric(quantities)
1102            quantity_names = ['Attribute']
1103
1104            #Make it a dictionary
1105            quantities = {quantity_names[0]: quantities}
1106
1107
1108        #Use keys if no names are specified
1109        if quantity_names is None:
1110            quantity_names = quantities.keys()
1111
1112
1113        #Check spatial info
1114        if vertex_coordinates is None:
1115            self.spatial = False
1116        else:   
1117            vertex_coordinates = ensure_numeric(vertex_coordinates)
1118
1119            assert triangles is not None, 'Triangles array must be specified'
1120            triangles = ensure_numeric(triangles)
1121            self.spatial = True           
1122           
1123
1124 
1125        #Save for use with statistics
1126        self.quantity_names = quantity_names       
1127        self.quantities = quantities       
1128        self.vertex_coordinates = vertex_coordinates
1129        self.interpolation_points = interpolation_points
1130        self.T = time[:]  # Time assumed to be relative to starttime
1131        self.index = 0    # Initial time index
1132        self.precomputed_values = {}
1133           
1134
1135
1136        #Precomputed spatial interpolation if requested
1137        if interpolation_points is not None:
1138            if self.spatial is False:
1139                raise 'Triangles and vertex_coordinates must be specified'
1140           
1141            try:
1142                self.interpolation_points = ensure_numeric(interpolation_points)
1143            except:
1144                msg = 'Interpolation points must be an N x 2 Numeric array '+\
1145                      'or a list of points\n'
1146                msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
1147                                      '...')
1148                raise msg
1149
1150
1151            m = len(self.interpolation_points)
1152            p = len(self.T)
1153           
1154            for name in quantity_names:
1155                self.precomputed_values[name] = zeros((p, m), Float)
1156
1157            #Build interpolator
1158            interpol = Interpolation(vertex_coordinates,
1159                                     triangles,
1160                                     point_coordinates = \
1161                                     self.interpolation_points,
1162                                     alpha = 0,
1163                                     precrop = False, 
1164                                     verbose = verbose)
1165
1166            if verbose: print 'Interpolate'
1167            for i, t in enumerate(self.T):
1168                #Interpolate quantities at this timestep
1169                if verbose and i%((p+10)/10)==0:
1170                    print ' time step %d of %d' %(i, p)
1171                   
1172                for name in quantity_names:
1173                    if len(quantities[name].shape) == 2:
1174                        result = interpol.interpolate(quantities[name][i,:])
1175                    else:
1176                       #Assume no time dependency
1177                       result = interpol.interpolate(quantities[name][:])
1178                       
1179                    self.precomputed_values[name][i, :] = result
1180                   
1181                       
1182
1183            #Report
1184            if verbose:
1185                print self.statistics()
1186                #self.print_statistics()
1187           
1188        else:
1189            #Store quantitites as is
1190            for name in quantity_names:
1191                self.precomputed_values[name] = quantities[name]
1192
1193
1194        #else:
1195        #    #Return an average, making this a time series
1196        #    for name in quantity_names:
1197        #        self.values[name] = zeros(len(self.T), Float)
1198        #
1199        #    if verbose: print 'Compute mean values'
1200        #    for i, t in enumerate(self.T):
1201        #        if verbose: print ' time step %d of %d' %(i, len(self.T))
1202        #        for name in quantity_names:
1203        #           self.values[name][i] = mean(quantities[name][i,:])
1204
1205
1206
1207
1208    def __repr__(self):
1209        #return 'Interpolation function (spatio-temporal)'
1210        return self.statistics()
1211   
1212
1213    def __call__(self, t, point_id = None, x = None, y = None):
1214        """Evaluate f(t), f(t, point_id) or f(t, x, y)
1215
1216        Inputs:
1217          t: time - Model time. Must lie within existing timesteps
1218          point_id: index of one of the preprocessed points.
1219          x, y:     Overrides location, point_id ignored
1220         
1221          If spatial info is present and all of x,y,point_id
1222          are None an exception is raised
1223                   
1224          If no spatial info is present, point_id and x,y arguments are ignored
1225          making f a function of time only.
1226
1227         
1228          FIXME: point_id could also be a slice
1229          FIXME: What if x and y are vectors?
1230          FIXME: What about f(x,y) without t?
1231        """
1232
1233        from math import pi, cos, sin, sqrt
1234        from Numeric import zeros, Float
1235        from utilities.numerical_tools import mean       
1236
1237        if self.spatial is True:
1238            if point_id is None:
1239                if x is None or y is None:
1240                    msg = 'Either point_id or x and y must be specified'
1241                    raise msg
1242            else:
1243                if self.interpolation_points is None:
1244                    msg = 'Interpolation_function must be instantiated ' +\
1245                          'with a list of interpolation points before parameter ' +\
1246                          'point_id can be used'
1247                    raise msg
1248
1249
1250        msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
1251        msg += ' does not match model time: %s\n' %t
1252        if t < self.T[0]: raise msg
1253        if t > self.T[-1]: raise msg
1254
1255        oldindex = self.index #Time index
1256
1257        #Find current time slot
1258        while t > self.T[self.index]: self.index += 1
1259        while t < self.T[self.index]: self.index -= 1
1260
1261        if t == self.T[self.index]:
1262            #Protect against case where t == T[-1] (last time)
1263            # - also works in general when t == T[i]
1264            ratio = 0
1265        else:
1266            #t is now between index and index+1
1267            ratio = (t - self.T[self.index])/\
1268                    (self.T[self.index+1] - self.T[self.index])
1269
1270        #Compute interpolated values
1271        q = zeros(len(self.quantity_names), Float)
1272
1273        for i, name in enumerate(self.quantity_names):
1274            Q = self.precomputed_values[name]
1275
1276            if self.spatial is False:
1277                #If there is no spatial info               
1278                assert len(Q.shape) == 1
1279
1280                Q0 = Q[self.index]
1281                if ratio > 0: Q1 = Q[self.index+1]
1282
1283            else:
1284                if x is not None and y is not None:
1285                    #Interpolate to x, y
1286                   
1287                    raise 'x,y interpolation not yet implemented'
1288                else:
1289                    #Use precomputed point
1290                    Q0 = Q[self.index, point_id]
1291                    if ratio > 0: Q1 = Q[self.index+1, point_id]
1292
1293            #Linear temporal interpolation   
1294            if ratio > 0:
1295                q[i] = Q0 + ratio*(Q1 - Q0)
1296            else:
1297                q[i] = Q0
1298
1299
1300        #Return vector of interpolated values
1301        #if len(q) == 1:
1302        #    return q[0]
1303        #else:
1304        #    return q
1305
1306
1307        #Return vector of interpolated values
1308        #FIXME:
1309        if self.spatial is True:
1310            return q
1311        else:
1312            #Replicate q according to x and y
1313            #This is e.g used for Wind_stress
1314            if x is None or y is None: 
1315                return q
1316            else:
1317                try:
1318                    N = len(x)
1319                except:
1320                    return q
1321                else:
1322                    from Numeric import ones, Float
1323                    #x is a vector - Create one constant column for each value
1324                    N = len(x)
1325                    assert len(y) == N, 'x and y must have same length'
1326                    res = []
1327                    for col in q:
1328                        res.append(col*ones(N, Float))
1329                       
1330                return res
1331
1332
1333    def statistics(self):
1334        """Output statistics about interpolation_function
1335        """
1336       
1337        vertex_coordinates = self.vertex_coordinates
1338        interpolation_points = self.interpolation_points               
1339        quantity_names = self.quantity_names
1340        quantities = self.quantities
1341        precomputed_values = self.precomputed_values                 
1342               
1343        x = vertex_coordinates[:,0]
1344        y = vertex_coordinates[:,1]               
1345
1346        str =  '------------------------------------------------\n'
1347        str += 'Interpolation_function (spatio-temporal) statistics:\n'
1348        str += '  Extent:\n'
1349        str += '    x in [%f, %f], len(x) == %d\n'\
1350               %(min(x), max(x), len(x))
1351        str += '    y in [%f, %f], len(y) == %d\n'\
1352               %(min(y), max(y), len(y))
1353        str += '    t in [%f, %f], len(t) == %d\n'\
1354               %(min(self.T), max(self.T), len(self.T))
1355        str += '  Quantities:\n'
1356        for name in quantity_names:
1357            q = quantities[name][:].flat
1358            str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
1359
1360        if interpolation_points is not None:   
1361            str += '  Interpolation points (xi, eta):'\
1362                   ' number of points == %d\n' %interpolation_points.shape[0]
1363            str += '    xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
1364                                            max(interpolation_points[:,0]))
1365            str += '    eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
1366                                             max(interpolation_points[:,1]))
1367            str += '  Interpolated quantities (over all timesteps):\n'
1368       
1369            for name in quantity_names:
1370                q = precomputed_values[name][:].flat
1371                str += '    %s at interpolation points in [%f, %f]\n'\
1372                       %(name, min(q), max(q))
1373        str += '------------------------------------------------\n'
1374
1375        return str
1376
1377        #FIXME: Delete
1378        #print '------------------------------------------------'
1379        #print 'Interpolation_function statistics:'
1380        #print '  Extent:'
1381        #print '    x in [%f, %f], len(x) == %d'\
1382        #      %(min(x), max(x), len(x))
1383        #print '    y in [%f, %f], len(y) == %d'\
1384        #      %(min(y), max(y), len(y))
1385        #print '    t in [%f, %f], len(t) == %d'\
1386        #      %(min(self.T), max(self.T), len(self.T))
1387        #print '  Quantities:'
1388        #for name in quantity_names:
1389        #    q = quantities[name][:].flat
1390        #    print '    %s in [%f, %f]' %(name, min(q), max(q))
1391        #print '  Interpolation points (xi, eta):'\
1392        #      ' number of points == %d ' %interpolation_points.shape[0]
1393        #print '    xi in [%f, %f]' %(min(interpolation_points[:,0]),
1394        #                             max(interpolation_points[:,0]))
1395        #print '    eta in [%f, %f]' %(min(interpolation_points[:,1]),
1396        #                              max(interpolation_points[:,1]))
1397        #print '  Interpolated quantities (over all timesteps):'
1398        #
1399        #for name in quantity_names:
1400        #    q = precomputed_values[name][:].flat
1401        #    print '    %s at interpolation points in [%f, %f]'\
1402        #          %(name, min(q), max(q))
1403        #print '------------------------------------------------'
1404
1405
1406#-------------------------------------------------------------
1407if __name__ == "__main__":
1408    """
1409    Load in a mesh and data points with attributes.
1410    Fit the attributes to the mesh.
1411    Save a new mesh file.
1412    """
1413    import os, sys
1414    usage = "usage: %s mesh_input.tsh point.xya mesh_output.tsh [expand|no_expand][vervose|non_verbose] [alpha] [display_errors|no_display_errors]"\
1415            %os.path.basename(sys.argv[0])
1416
1417    if len(sys.argv) < 4:
1418        print usage
1419    else:
1420        mesh_file = sys.argv[1]
1421        point_file = sys.argv[2]
1422        mesh_output_file = sys.argv[3]
1423
1424        expand_search = False
1425        if len(sys.argv) > 4:
1426            if sys.argv[4][0] == "e" or sys.argv[4][0] == "E":
1427                expand_search = True
1428            else:
1429                expand_search = False
1430
1431        verbose = False
1432        if len(sys.argv) > 5:
1433            if sys.argv[5][0] == "n" or sys.argv[5][0] == "N":
1434                verbose = False
1435            else:
1436                verbose = True
1437
1438        if len(sys.argv) > 6:
1439            alpha = sys.argv[6]
1440        else:
1441            alpha = DEFAULT_ALPHA
1442
1443        # This is used more for testing
1444        if len(sys.argv) > 7:
1445            if sys.argv[7][0] == "n" or sys.argv[5][0] == "N":
1446                display_errors = False
1447            else:
1448                display_errors = True
1449           
1450        t0 = time.time()
1451        try:
1452            fit_to_mesh_file(mesh_file,
1453                         point_file,
1454                         mesh_output_file,
1455                         alpha,
1456                         verbose= verbose,
1457                         expand_search = expand_search,
1458                         display_errors = display_errors)
1459        except IOError,e:
1460            import sys; sys.exit(1)
1461
1462        print 'That took %.2f seconds' %(time.time()-t0)
1463
Note: See TracBrowser for help on using the repository browser.