source: branches/numpy/obsolete_code/least_squares.py @ 8353

Last change on this file since 8353 was 3564, checked in by ole, 18 years ago

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