Context Navigation

source: inundation/fit_interpolate/fit.py @ 3085

Last change on this file since 3085 was 3030, checked in by duncan, 19 years ago
improve warning message
File size: 20.5 KB

Rev	Line
[2802]	1	"""Least squares fitting.
	2
	3	Implements a penalised least-squares fit.
	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.
[2897]	18
	19	TO DO
	20	* test geo_ref, geo_spatial
[2802]	21	"""
	22
[3012]	23	from Numeric import zeros, Float, ArrayType,take
[2897]	24
	25	from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
[2802]	26	from fit_interpolate.general_fit_interpolate import FitInterpolate
[2897]	27	from utilities.sparse import Sparse, Sparse_CSR
	28	from utilities.polygon import in_and_outside_polygon
	29	from fit_interpolate.search_functions import search_tree_of_vertices
	30	from utilities.cg_solve import conjugate_gradient
	31	from utilities.numerical_tools import ensure_numeric, gradient
[2802]	32
[2897]	33	import exceptions
	34	class ToFewPointsError(exceptions.Exception): pass
[3014]	35	class VertsWithNoTrianglesError(exceptions.Exception): pass
[2897]	36
[2802]	37	DEFAULT_ALPHA = 0.001
	38
	39
	40	class Fit(FitInterpolate):
	41
	42	def __init__(self,
	43	vertex_coordinates,
	44	triangles,
	45	mesh_origin=None,
	46	alpha = None,
	47	verbose=False,
	48	max_vertices_per_cell=30):
	49
	50
	51	"""
	52	Fit data at points to the vertices of a mesh.
	53
	54	Inputs:
	55
	56	vertex_coordinates: List of coordinate pairs [xi, eta] of
	57	points constituting a mesh (or an m x 2 Numeric array or
	58	a geospatial object)
	59	Points may appear multiple times
	60	(e.g. if vertices have discontinuities)
	61
	62	triangles: List of 3-tuples (or a Numeric array) of
	63	integers representing indices of all vertices in the mesh.
	64
	65	mesh_origin: A geo_reference object or 3-tuples consisting of
	66	UTM zone, easting and northing.
	67	If specified vertex coordinates are assumed to be
	68	relative to their respective origins.
	69
	70	max_vertices_per_cell: Number of vertices in a quad tree cell
	71	at which the cell is split into 4.
	72
	73	Note: Don't supply a vertex coords as a geospatial object and
	74	a mesh origin, since geospatial has its own mesh origin.
	75	"""
	76	# Initialise variabels
	77
	78	if alpha is None:
[2897]	79
[2802]	80	self.alpha = DEFAULT_ALPHA
	81	else:
	82	self.alpha = alpha
	83	FitInterpolate.__init__(self,
	84	vertex_coordinates,
	85	triangles,
	86	mesh_origin,
	87	verbose,
	88	max_vertices_per_cell)
	89
	90	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (vertices)
	91
[2897]	92	self.AtA = None
	93	self.Atz = None
[2802]	94
[2897]	95	self.point_count = 0
	96	if self.alpha <> 0:
	97	if verbose: print 'Building smoothing matrix'
	98	self._build_smoothing_matrix_D()
	99
[2802]	100	def _build_coefficient_matrix_B(self,
	101	verbose = False):
[2897]	102	"""
	103	Build final coefficient matrix
[2802]	104
[2897]	105	Precon
	106	If alpha is not zero, matrix D has been built
	107	Matrix Ata has been built
	108	"""
[2802]	109
[2897]	110	if self.alpha <> 0:
	111	#if verbose: print 'Building smoothing matrix'
	112	#self._build_smoothing_matrix_D()
	113	self.B = self.AtA + self.alpha*self.D
	114	else:
	115	self.B = self.AtA
	116
	117	#Convert self.B matrix to CSR format for faster matrix vector
	118	self.B = Sparse_CSR(self.B)
	119
[2802]	120	def _build_smoothing_matrix_D(self):
	121	"""Build m x m smoothing matrix, where
	122	m is the number of basis functions phi_k (one per vertex)
	123
	124	The smoothing matrix is defined as
	125
	126	D = D1 + D2
	127
	128	where
	129
	130	[D1]_{k,l} = \int_\Omega
	131	\frac{\partial \phi_k}{\partial x}
	132	\frac{\partial \phi_l}{\partial x}\,
	133	dx dy
	134
	135	[D2]_{k,l} = \int_\Omega
	136	\frac{\partial \phi_k}{\partial y}
	137	\frac{\partial \phi_l}{\partial y}\,
	138	dx dy
	139
	140
	141	The derivatives \frac{\partial \phi_k}{\partial x},
	142	\frac{\partial \phi_k}{\partial x} for a particular triangle
	143	are obtained by computing the gradient a_k, b_k for basis function k
	144	"""
[2897]	145
	146	#FIXME: algorithm might be optimised by computing local 9x9
	147	#"element stiffness matrices:
[2802]	148
[2897]	149	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
	150
	151	self.D = Sparse(m,m)
	152
	153	#For each triangle compute contributions to D = D1+D2
	154	for i in range(len(self.mesh)):
	155
	156	#Get area
	157	area = self.mesh.areas[i]
	158
	159	#Get global vertex indices
	160	v0 = self.mesh.triangles[i,0]
	161	v1 = self.mesh.triangles[i,1]
	162	v2 = self.mesh.triangles[i,2]
	163
	164	#Get the three vertex_points
	165	xi0 = self.mesh.get_vertex_coordinate(i, 0)
	166	xi1 = self.mesh.get_vertex_coordinate(i, 1)
	167	xi2 = self.mesh.get_vertex_coordinate(i, 2)
	168
	169	#Compute gradients for each vertex
	170	a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
	171	1, 0, 0)
	172
	173	a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
	174	0, 1, 0)
	175
	176	a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
	177	0, 0, 1)
	178
	179	#Compute diagonal contributions
	180	self.D[v0,v0] += (a0a0 + b0b0)*area
	181	self.D[v1,v1] += (a1a1 + b1b1)*area
	182	self.D[v2,v2] += (a2a2 + b2b2)*area
	183
	184	#Compute contributions for basis functions sharing edges
	185	e01 = (a0a1 + b0b1)*area
	186	self.D[v0,v1] += e01
	187	self.D[v1,v0] += e01
	188
	189	e12 = (a1a2 + b1b2)*area
	190	self.D[v1,v2] += e12
	191	self.D[v2,v1] += e12
	192
	193	e20 = (a2a0 + b2b0)*area
	194	self.D[v2,v0] += e20
	195	self.D[v0,v2] += e20
	196
	197
	198	def get_D(self):
	199	return self.D.todense()
	200
	201
[2802]	202	def _build_matrix_AtA_Atz(self,
	203	point_coordinates,
	204	z,
	205	verbose = False):
	206	"""Build:
	207	AtA m x m interpolation matrix, and,
	208	Atz m x a interpolation matrix where,
	209	m is the number of basis functions phi_k (one per vertex)
	210	a is the number of data attributes
	211
	212	This algorithm uses a quad tree data structure for fast binning of
	213	data points.
	214
[2939]	215	If Ata is None, the matrices AtA and Atz are created.
	216
	217	This function can be called again and again, with sub-sets of
	218	the point coordinates. Call fit to get the results.
	219
[2802]	220	Preconditions
	221	z and points are numeric
	222	Point_coordindates and mesh vertices have the same origin.
[2897]	223
	224	The number of attributes of the data points does not change
[2802]	225	"""
[2897]	226	#Build n x m interpolation matrix
	227
	228	if self.AtA == None:
	229	# AtA and Atz need ot be initialised.
	230	m = self.mesh.coordinates.shape[0] #Nbr of vertices
	231	if len(z.shape) > 1:
	232	att_num = z.shape[1]
	233	self.Atz = zeros((m,att_num), Float)
	234	else:
	235	att_num = 1
	236	self.Atz = zeros((m,), Float)
	237	assert z.shape[0] == point_coordinates.shape[0]
	238
	239	self.AtA = Sparse(m,m)
	240	self.point_count += point_coordinates.shape[0]
	241	#print "_build_matrix_AtA_Atz - self.point_count", self.point_count
[2802]	242	if verbose: print 'Getting indices inside mesh boundary'
[2897]	243	#print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon()
[2802]	244	self.inside_poly_indices, self.outside_poly_indices = \
	245	in_and_outside_polygon(point_coordinates,
	246	self.mesh.get_boundary_polygon(),
	247	closed = True, verbose = verbose)
	248	#print "self.inside_poly_indices",self.inside_poly_indices
	249	#print "self.outside_poly_indices",self.outside_poly_indices
	250
	251
[2897]	252	n = len(self.inside_poly_indices)
[2802]	253	#Compute matrix elements for points inside the mesh
	254	for i in self.inside_poly_indices:
	255	#For each data_coordinate point
	256	if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
	257	x = point_coordinates[i]
	258	element_found, sigma0, sigma1, sigma2, k = \
	259	search_tree_of_vertices(self.root, self.mesh, x)
[2939]	260
[2802]	261	if element_found is True:
	262	j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
	263	j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
	264	j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
	265
	266	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
	267	js = [j0,j1,j2]
	268
	269	for j in js:
[2897]	270	self.Atz[j] += sigmas[j]*z[i]
	271	#print "self.Atz building", self.Atz
	272	#print "self.Atz[j]", self.Atz[j]
	273	#print " sigmas[j]", sigmas[j]
	274	#print "z[i]",z[i]
	275	#print "result", sigmas[j]*z[i]
	276
[2802]	277	for k in js:
[2897]	278	self.AtA[j,k] += sigmas[j]*sigmas[k]
[2802]	279	else:
	280	msg = 'Could not find triangle for point', x
	281	raise Exception(msg)
	282
	283
[2897]	284	def fit(self, point_coordinates=None, z=None,
	285	verbose = False,
	286	point_origin = None):
[2802]	287	"""Fit a smooth surface to given 1d array of data points z.
	288
	289	The smooth surface is computed at each vertex in the underlying
	290	mesh using the formula given in the module doc string.
	291
	292	Inputs:
	293	point_coordinates: The co-ordinates of the data points.
	294	List of coordinate pairs [x, y] of
	295	data points or an nx2 Numeric array or a Geospatial_data object
	296	z: Single 1d vector or array of data at the point_coordinates.
	297
	298	"""
[2897]	299	if point_coordinates is None:
	300	assert self.AtA <> None
	301	assert self.Atz <> None
	302	#FIXME (DSG) - do a message
	303	else:
	304	point_coordinates = ensure_absolute(point_coordinates,
	305	geo_reference=point_origin)
	306	self.build_fit_subset(point_coordinates, z, verbose)
[3012]	307
[2897]	308	#Check sanity
	309	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
	310	n = self.point_count
	311	if n<m and self.alpha == 0.0:
	312	msg = 'ERROR (least_squares): Too few data points\n'
	313	msg += 'There are only %d data points and alpha == 0. ' %n
	314	msg += 'Need at least %d\n' %m
	315	msg += 'Alternatively, set smoothing parameter alpha to a small '
	316	msg += 'positive value,\ne.g. 1.0e-3.'
	317	raise ToFewPointsError(msg)
[2802]	318
[2897]	319	self._build_coefficient_matrix_B(verbose)
[3012]	320	loners = self.mesh.get_lone_vertices()
[3014]	321	# FIXME - make this as error message.
	322	# test with
	323	# Not_yet_test_smooth_att_to_mesh_with_excess_verts.
	324	if len(loners)>0:
	325	msg = 'WARNING: (least_squares): \nVertices with no triangles\n'
	326	msg += 'All vertices should be part of a triangle.\n'
	327	msg += 'In the future this will be inforced.\n'
[3030]	328	msg += 'The following vertices are not part of a triangle;\n'
[3014]	329	msg += str(loners)
	330	print msg
	331	#raise VertsWithNoTrianglesError(msg)
[3012]	332
[3013]	333
[2897]	334	return conjugate_gradient(self.B, self.Atz, self.Atz,
	335	imax=2*len(self.Atz) )
	336
[2802]	337
[2897]	338	def build_fit_subset(self, point_coordinates, z,
	339	verbose = False):
[2802]	340	"""Fit a smooth surface to given 1d array of data points z.
	341
	342	The smooth surface is computed at each vertex in the underlying
	343	mesh using the formula given in the module doc string.
	344
	345	Inputs:
	346	point_coordinates: The co-ordinates of the data points.
	347	List of coordinate pairs [x, y] of
	348	data points or an nx2 Numeric array or a Geospatial_data object
	349	z: Single 1d vector or array of data at the point_coordinates.
	350
	351	"""
	352	#Note: Don't get the z info from Geospatial_data.attributes yet.
[2897]	353	# If we did fit would have to handle attribute title info.
[2802]	354
	355	#FIXME(DSG-DSG): Check that the vert and point coords
	356	#have the same zone.
	357	if isinstance(point_coordinates,Geospatial_data):
[2897]	358	point_coordinates = point_coordinates.get_data_points( \
[2802]	359	absolute = True)
	360
	361	#Convert input to Numeric arrays
	362	z = ensure_numeric(z, Float)
	363	point_coordinates = ensure_numeric(point_coordinates, Float)
	364
[2897]	365	self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
[2802]	366
	367
	368	############################################################################
	369
	370	def fit_to_mesh(vertex_coordinates,
	371	triangles,
	372	point_coordinates,
	373	point_attributes,
	374	alpha = DEFAULT_ALPHA,
	375	verbose = False,
	376	acceptable_overshoot = 1.01,
[2897]	377	mesh_origin = None,
[2802]	378	data_origin = None,
	379	use_cache = False):
	380	"""
	381	Fit a smooth surface to a triangulation,
	382	given data points with attributes.
	383
	384
	385	Inputs:
[2897]	386	vertex_coordinates: List of coordinate pairs [xi, eta] of
	387	points constituting a mesh (or an m x 2 Numeric array or
	388	a geospatial object)
	389	Points may appear multiple times
	390	(e.g. if vertices have discontinuities)
[2802]	391
	392	triangles: List of 3-tuples (or a Numeric array) of
	393	integers representing indices of all vertices in the mesh.
	394
	395	point_coordinates: List of coordinate pairs [x, y] of data points
	396	(or an nx2 Numeric array)
	397
	398	alpha: Smoothing parameter.
	399
	400	acceptable overshoot: controls the allowed factor by which fitted values
	401	may exceed the value of input data. The lower limit is defined
	402	as min(z) - acceptable_overshoot*delta z and upper limit
	403	as max(z) + acceptable_overshoot*delta z
[2897]	404
	405	mesh_origin: A geo_reference object or 3-tuples consisting of
	406	UTM zone, easting and northing.
	407	If specified vertex coordinates are assumed to be
	408	relative to their respective origins.
[2802]	409
	410
[2897]	411	point_attributes: Vector or array of data at the
	412	point_coordinates.
[2802]	413
	414	"""
[2897]	415	#Since this is a wrapper for fit, lets handle the geo_spatial att's
[2802]	416	if use_cache is True:
	417	interp = cache(_fit,
	418	(vertex_coordinates,
	419	triangles),
	420	{'verbose': verbose,
[2897]	421	'mesh_origin': mesh_origin,
	422	'alpha':alpha},
[2802]	423	verbose = verbose)
	424
	425	else:
[2897]	426	interp = Fit(vertex_coordinates,
	427	triangles,
	428	verbose = verbose,
	429	mesh_origin = mesh_origin,
	430	alpha=alpha)
[2802]	431
[2897]	432	vertex_attributes = interp.fit(point_coordinates,
	433	point_attributes,
	434	point_origin = data_origin,
	435	verbose = verbose)
[2802]	436
	437
[2897]	438	# Add the value checking stuff that's in least squares.
	439	# Maybe this stuff should get pushed down into Fit.
	440	# at least be a method of Fit.
	441	# Or intigrate it into the fit method, saving teh max and min's
	442	# as att's.
[2802]	443
	444	return vertex_attributes
	445
	446
[2939]	447	def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
	448	alpha=DEFAULT_ALPHA, verbose= False,
	449	display_errors = True):
	450	"""
	451	Given a mesh file (tsh) and a point attribute file (xya), fit
	452	point attributes to the mesh and write a mesh file with the
	453	results.
	454
	455
	456	If data_origin is not None it is assumed to be
	457	a 3-tuple with geo referenced
	458	UTM coordinates (zone, easting, northing)
	459
	460	NOTE: Throws IOErrors, for a variety of file problems.
	461
	462	"""
	463
	464	# Question
	465	# should data_origin and mesh_origin be passed in?
	466	# No they should be in the data structure
	467	#
	468	#Should the origin of the mesh be changed using this function?
	469	# That is overloading this function. Have it as a seperate
	470	# method, at least initially.
	471
	472	from load_mesh.loadASCII import import_mesh_file, \
	473	import_points_file, export_mesh_file, \
	474	concatinate_attributelist
	475
	476	# FIXME: Use geospatial instead of import_points_file
	477	try:
	478	mesh_dict = import_mesh_file(mesh_file)
	479	except IOError,e:
	480	if display_errors:
	481	print "Could not load bad file. ", e
	482	raise IOError #Re-raise exception
	483
	484	vertex_coordinates = mesh_dict['vertices']
	485	triangles = mesh_dict['triangles']
	486	if type(mesh_dict['vertex_attributes']) == ArrayType:
	487	old_point_attributes = mesh_dict['vertex_attributes'].tolist()
	488	else:
	489	old_point_attributes = mesh_dict['vertex_attributes']
	490
	491	if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
	492	old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
	493	else:
	494	old_title_list = mesh_dict['vertex_attribute_titles']
	495
	496	if verbose: print 'tsh file %s loaded' %mesh_file
	497
	498	# load in the .pts file
	499	try:
	500	point_dict = import_points_file(point_file, verbose=verbose)
	501	except IOError,e:
	502	if display_errors:
	503	print "Could not load bad file. ", e
	504	raise IOError #Re-raise exception
	505
	506	point_coordinates = point_dict['pointlist']
	507	title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
	508
	509	if point_dict.has_key('geo_reference') and not point_dict['geo_reference'] is None:
	510	data_origin = point_dict['geo_reference'].get_origin()
	511	else:
	512	data_origin = (56, 0, 0) #FIXME(DSG-DSG)
	513
	514	if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
	515	mesh_origin = mesh_dict['geo_reference'].get_origin()
	516	else:
	517	mesh_origin = (56, 0, 0) #FIXME(DSG-DSG)
	518
	519	if verbose: print "points file loaded"
	520	if verbose: print "fitting to mesh"
	521	f = fit_to_mesh(vertex_coordinates,
	522	triangles,
	523	point_coordinates,
	524	point_attributes,
	525	alpha = alpha,
	526	verbose = verbose,
	527	data_origin = data_origin,
	528	mesh_origin = mesh_origin)
	529	if verbose: print "finished fitting to mesh"
	530
	531	# convert array to list of lists
	532	new_point_attributes = f.tolist()
	533	#FIXME have this overwrite attributes with the same title - DSG
	534	#Put the newer attributes last
	535	if old_title_list <> []:
	536	old_title_list.extend(title_list)
	537	#FIXME can this be done a faster way? - DSG
	538	for i in range(len(old_point_attributes)):
	539	old_point_attributes[i].extend(new_point_attributes[i])
	540	mesh_dict['vertex_attributes'] = old_point_attributes
	541	mesh_dict['vertex_attribute_titles'] = old_title_list
	542	else:
	543	mesh_dict['vertex_attributes'] = new_point_attributes
	544	mesh_dict['vertex_attribute_titles'] = title_list
	545
	546	if verbose: print "exporting to file ", mesh_output_file
	547
	548	try:
	549	export_mesh_file(mesh_output_file, mesh_dict)
	550	except IOError,e:
	551	if display_errors:
	552	print "Could not write file. ", e
	553	raise IOError
	554
	555
[2802]	556	def _fit(args, *kwargs):
	557	"""Private function for use with caching. Reason is that classes
	558	may change their byte code between runs which is annoying.
	559	"""
	560
	561	return Fit(args, *kwargs)
	562

Note: See TracBrowser for help on using the repository browser.

Download in other formats: