"""Least squares fitting. Implements a penalised least-squares fit. The penalty term (or smoothing term) is controlled by the smoothing parameter alpha. With a value of alpha=0, the fit function will attempt to interpolate as closely as possible in the least-squares sense. With values alpha > 0, a certain amount of smoothing will be applied. A positive alpha is essential in cases where there are too few data points. A negative alpha is not allowed. A typical value of alpha is 1.0e-6 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou Geoscience Australia, 2004. """ from geospatial_data.geospatial_data import Geospatial_data from fit_interpolate.general_fit_interpolate import FitInterpolate DEFAULT_ALPHA = 0.001 class Fit(FitInterpolate): def __init__(self, vertex_coordinates, triangles, mesh_origin=None, alpha = None, verbose=False, max_vertices_per_cell=30): """ Fit data at points to the vertices of a mesh. Inputs: vertex_coordinates: List of coordinate pairs [xi, eta] of points constituting a mesh (or an m x 2 Numeric array or a geospatial object) Points may appear multiple times (e.g. if vertices have discontinuities) triangles: List of 3-tuples (or a Numeric array) of integers representing indices of all vertices in the mesh. mesh_origin: A geo_reference object or 3-tuples consisting of UTM zone, easting and northing. If specified vertex coordinates are assumed to be relative to their respective origins. max_vertices_per_cell: Number of vertices in a quad tree cell at which the cell is split into 4. Note: Don't supply a vertex coords as a geospatial object and a mesh origin, since geospatial has its own mesh origin. """ # Initialise variabels #self._A_can_be_reused = False #self._point_coordinates = None if alpha is None: self.alpha = DEFAULT_ALPHA else: self.alpha = alpha FitInterpolate.__init__(self, vertex_coordinates, triangles, mesh_origin, verbose, max_vertices_per_cell) m = self.mesh.coordinates.shape[0] #Nbr of basis functions (vertices) #Build Atz and AtA matrix self.AtA = Sparse(m,m) self.Atz = zeros((m,att_num), Float) def _build_coefficient_matrix_B(self, verbose = False): """Build final coefficient matrix""" def _build_smoothing_matrix_D(self): """Build m x m smoothing matrix, where m is the number of basis functions phi_k (one per vertex) The smoothing matrix is defined as D = D1 + D2 where [D1]_{k,l} = \int_\Omega \frac{\partial \phi_k}{\partial x} \frac{\partial \phi_l}{\partial x}\, dx dy [D2]_{k,l} = \int_\Omega \frac{\partial \phi_k}{\partial y} \frac{\partial \phi_l}{\partial y}\, dx dy The derivatives \frac{\partial \phi_k}{\partial x}, \frac{\partial \phi_k}{\partial x} for a particular triangle are obtained by computing the gradient a_k, b_k for basis function k """ def _build_matrix_AtA_Atz(self, point_coordinates, z, verbose = False): """Build: AtA m x m interpolation matrix, and, Atz m x a interpolation matrix where, m is the number of basis functions phi_k (one per vertex) a is the number of data attributes This algorithm uses a quad tree data structure for fast binning of data points. Preconditions z and points are numeric Point_coordindates and mesh vertices have the same origin. """ if verbose: print 'Getting indices inside mesh boundary' #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon() self.inside_poly_indices, self.outside_poly_indices = \ in_and_outside_polygon(point_coordinates, self.mesh.get_boundary_polygon(), closed = True, verbose = verbose) #print "self.inside_poly_indices",self.inside_poly_indices #print "self.outside_poly_indices",self.outside_poly_indices #Compute matrix elements for points inside the mesh for i in self.inside_poly_indices: #For each data_coordinate point if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n) x = point_coordinates[i] element_found, sigma0, sigma1, sigma2, k = \ search_tree_of_vertices(self.root, self.mesh, x) #Update interpolation matrix A if necessary if element_found is True: #Assign values to matrix A j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0 j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1 j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2 sigmas = {j0:sigma0, j1:sigma1, j2:sigma2} js = [j0,j1,j2] for j in js: self.Atz[j] += sigmas[j]*z[i] for k in js: if interp_only == False: self.AtA[j,k] += sigmas[j]*sigmas[k] else: msg = 'Could not find triangle for point', x raise Exception(msg) def fit(self, point_coordinates=point_coordinates, z=z): """Fit a smooth surface to given 1d array of data points z. The smooth surface is computed at each vertex in the underlying mesh using the formula given in the module doc string. Inputs: point_coordinates: The co-ordinates of the data points. List of coordinate pairs [x, y] of data points or an nx2 Numeric array or a Geospatial_data object z: Single 1d vector or array of data at the point_coordinates. """ # build ata and atz # solve fit def build_fit_subset(self, point_coordinates, z): """Fit a smooth surface to given 1d array of data points z. The smooth surface is computed at each vertex in the underlying mesh using the formula given in the module doc string. Inputs: point_coordinates: The co-ordinates of the data points. List of coordinate pairs [x, y] of data points or an nx2 Numeric array or a Geospatial_data object z: Single 1d vector or array of data at the point_coordinates. """ #Note: Don't get the z info from Geospatial_data.attributes yet. # That means fit has to handle attribute title info. #FIXME(DSG-DSG): Check that the vert and point coords #have the same zone. if isinstance(point_coordinates,Geospatial_data): point_coordinates = vertex_coordinates.get_data_points( \ absolute = True) #Convert input to Numeric arrays z = ensure_numeric(z, Float) point_coordinates = ensure_numeric(point_coordinates, Float) ############################################################################ def fit_to_mesh(vertex_coordinates, triangles, point_coordinates, point_attributes, alpha = DEFAULT_ALPHA, verbose = False, acceptable_overshoot = 1.01, expand_search = False, data_origin = None, mesh_origin = None, precrop = False, use_cache = False): """ Fit a smooth surface to a triangulation, given data points with attributes. Inputs: vertex_coordinates: List of coordinate pairs [xi, eta] of points constituting mesh (or a an m x 2 Numeric array) triangles: List of 3-tuples (or a Numeric array) of integers representing indices of all vertices in the mesh. point_coordinates: List of coordinate pairs [x, y] of data points (or an nx2 Numeric array) alpha: Smoothing parameter. acceptable overshoot: controls the allowed factor by which fitted values may exceed the value of input data. The lower limit is defined as min(z) - acceptable_overshoot*delta z and upper limit as max(z) + acceptable_overshoot*delta z point_attributes: Vector or array of data at the point_coordinates. data_origin and mesh_origin are 3-tuples consisting of UTM zone, easting and northing. If specified point coordinates and vertex coordinates are assumed to be relative to their respective origins. """ if use_cache is True: interp = cache(_fit, (vertex_coordinates, triangles), {'verbose': verbose, 'mesh_origin': mesh_origin}, verbose = verbose) else: interp = Interpolation(vertex_coordinates, triangles, verbose = verbose, mesh_origin = mesh_origin) vertex_attributes = interp.fit_points(point_attributes, verbose = verbose) # point_coordinates, # data_origin = data_origin, #Sanity check point_coordinates = ensure_numeric(point_coordinates) vertex_coordinates = ensure_numeric(vertex_coordinates) #Data points X = point_coordinates[:,0] Y = point_coordinates[:,1] Z = ensure_numeric(point_attributes) if len(Z.shape) == 1: Z = Z[:, NewAxis] #Data points inside mesh boundary indices = interp.point_indices if indices is not None: Xc = take(X, indices) Yc = take(Y, indices) Zc = take(Z, indices) else: Xc = X Yc = Y Zc = Z #Vertex coordinates Xi = vertex_coordinates[:,0] Eta = vertex_coordinates[:,1] Zeta = ensure_numeric(vertex_attributes) if len(Zeta.shape) == 1: Zeta = Zeta[:, NewAxis] for i in range(Zeta.shape[1]): #For each attribute zeta = Zeta[:,i] z = Z[:,i] zc = Zc[:,i] max_zc = max(zc) min_zc = min(zc) delta_zc = max_zc-min_zc upper_limit = max_zc + delta_zc*acceptable_overshoot lower_limit = min_zc - delta_zc*acceptable_overshoot if max(zeta) > upper_limit or min(zeta) < lower_limit: msg = 'Least sqares produced values outside the allowed ' msg += 'range [%f, %f].\n' %(lower_limit, upper_limit) msg += 'z in [%f, %f], zeta in [%f, %f].\n' %(min_zc, max_zc, min(zeta), max(zeta)) msg += 'If greater range is needed, increase the value of ' msg += 'acceptable_fit_overshoot (currently %.2f).\n' %(acceptable_overshoot) offending_vertices = (zeta > upper_limit or zeta < lower_limit) Xi_c = compress(offending_vertices, Xi) Eta_c = compress(offending_vertices, Eta) offending_coordinates = concatenate((Xi_c[:, NewAxis], Eta_c[:, NewAxis]), axis=1) msg += 'Offending locations:\n %s' %(offending_coordinates) raise FittingError, msg if verbose: print '+------------------------------------------------' print 'Least squares statistics' print '+------------------------------------------------' print 'points: %d points' %(len(z)) print ' x in [%f, %f]'%(min(X), max(X)) print ' y in [%f, %f]'%(min(Y), max(Y)) print ' z in [%f, %f]'%(min(z), max(z)) print if indices is not None: print 'Cropped points: %d points' %(len(zc)) print ' x in [%f, %f]'%(min(Xc), max(Xc)) print ' y in [%f, %f]'%(min(Yc), max(Yc)) print ' z in [%f, %f]'%(min(zc), max(zc)) print print 'Mesh: %d vertices' %(len(zeta)) print ' xi in [%f, %f]'%(min(Xi), max(Xi)) print ' eta in [%f, %f]'%(min(Eta), max(Eta)) print ' zeta in [%f, %f]'%(min(zeta), max(zeta)) print '+------------------------------------------------' return vertex_attributes def _fit(*args, **kwargs): """Private function for use with caching. Reason is that classes may change their byte code between runs which is annoying. """ return Fit(*args, **kwargs)