Context Navigation

source: branches/numpy/anuga/fit_interpolate/fit.py @ 6553

Last change on this file since 6553 was 6553, checked in by rwilson, 16 years ago
Merged trunk into numpy, all tests and validations work.
File size: 25.8 KB

Line
1	"""Least squares fitting.
2
3	Implements a penalised least-squares fit.
4	putting point data onto the mesh.
5
6	The penalty term (or smoothing term) is controlled by the smoothing
7	parameter alpha.
8	With a value of alpha=0, the fit function will attempt
9	to interpolate as closely as possible in the least-squares sense.
10	With values alpha > 0, a certain amount of smoothing will be applied.
11	A positive alpha is essential in cases where there are too few
12	data points.
13	A negative alpha is not allowed.
14	A typical value of alpha is 1.0e-6
15
16
17	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
18	Geoscience Australia, 2004.
19
20	TO DO
21	* test geo_ref, geo_spatial
22
23	IDEAS
24	* (DSG-) Change the interface of fit, so a domain object can
25	be passed in. (I don't know if this is feasible). If could
26	save time/memory.
27	"""
28	import types
29
30	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
31	from anuga.caching import cache
32	from anuga.geospatial_data.geospatial_data import Geospatial_data, \
33	ensure_absolute
34	from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
35	from anuga.utilities.sparse import Sparse, Sparse_CSR
36	from anuga.utilities.polygon import inside_polygon, is_inside_polygon
37	from anuga.fit_interpolate.search_functions import search_tree_of_vertices
38
39	from anuga.utilities.cg_solve import conjugate_gradient
40	from anuga.utilities.numerical_tools import ensure_numeric, gradient
41	from anuga.config import default_smoothing_parameter as DEFAULT_ALPHA
42
43	import exceptions
44	class TooFewPointsError(exceptions.Exception): pass
45	class VertsWithNoTrianglesError(exceptions.Exception): pass
46
47	import numpy as num
48
49
50	class Fit(FitInterpolate):
51
52	def __init__(self,
53	vertex_coordinates=None,
54	triangles=None,
55	mesh=None,
56	mesh_origin=None,
57	alpha = None,
58	verbose=False,
59	max_vertices_per_cell=None):
60
61
62	"""
63	Fit data at points to the vertices of a mesh.
64
65	Inputs:
66
67	vertex_coordinates: List of coordinate pairs [xi, eta] of
68	points constituting a mesh (or an m x 2 numeric array or
69	a geospatial object)
70	Points may appear multiple times
71	(e.g. if vertices have discontinuities)
72
73	triangles: List of 3-tuples (or a numeric array) of
74	integers representing indices of all vertices in the mesh.
75
76	mesh_origin: A geo_reference object or 3-tuples consisting of
77	UTM zone, easting and northing.
78	If specified vertex coordinates are assumed to be
79	relative to their respective origins.
80
81	max_vertices_per_cell: Number of vertices in a quad tree cell
82	at which the cell is split into 4.
83
84	Note: Don't supply a vertex coords as a geospatial object and
85	a mesh origin, since geospatial has its own mesh origin.
86
87
88	Usage,
89	To use this in a blocking way, call build_fit_subset, with z info,
90	and then fit, with no point coord, z info.
91
92	"""
93	# Initialise variabels
94	if alpha is None:
95	self.alpha = DEFAULT_ALPHA
96	else:
97	self.alpha = alpha
98
99	FitInterpolate.__init__(self,
100	vertex_coordinates,
101	triangles,
102	mesh,
103	mesh_origin,
104	verbose,
105	max_vertices_per_cell)
106
107	m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
108
109	self.AtA = None
110	self.Atz = None
111
112	self.point_count = 0
113	if self.alpha <> 0:
114	if verbose: print 'Building smoothing matrix'
115	self._build_smoothing_matrix_D()
116
117	bd_poly = self.mesh.get_boundary_polygon()
118	self.mesh_boundary_polygon = ensure_numeric(bd_poly)
119
120	def _build_coefficient_matrix_B(self,
121	verbose = False):
122	"""
123	Build final coefficient matrix
124
125	Precon
126	If alpha is not zero, matrix D has been built
127	Matrix Ata has been built
128	"""
129
130	if self.alpha <> 0:
131	#if verbose: print 'Building smoothing matrix'
132	#self._build_smoothing_matrix_D()
133	self.B = self.AtA + self.alpha*self.D
134	else:
135	self.B = self.AtA
136
137	# Convert self.B matrix to CSR format for faster matrix vector
138	self.B = Sparse_CSR(self.B)
139
140	def _build_smoothing_matrix_D(self):
141	"""Build m x m smoothing matrix, where
142	m is the number of basis functions phi_k (one per vertex)
143
144	The smoothing matrix is defined as
145
146	D = D1 + D2
147
148	where
149
150	[D1]_{k,l} = \int_\Omega
151	\frac{\partial \phi_k}{\partial x}
152	\frac{\partial \phi_l}{\partial x}\,
153	dx dy
154
155	[D2]_{k,l} = \int_\Omega
156	\frac{\partial \phi_k}{\partial y}
157	\frac{\partial \phi_l}{\partial y}\,
158	dx dy
159
160
161	The derivatives \frac{\partial \phi_k}{\partial x},
162	\frac{\partial \phi_k}{\partial x} for a particular triangle
163	are obtained by computing the gradient a_k, b_k for basis function k
164	"""
165
166	# FIXME: algorithm might be optimised by computing local 9x9
167	# "element stiffness matrices:
168
169	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
170
171	self.D = Sparse(m,m)
172
173	# For each triangle compute contributions to D = D1+D2
174	for i in range(len(self.mesh)):
175
176	# Get area
177	area = self.mesh.areas[i]
178
179	# Get global vertex indices
180	v0 = self.mesh.triangles[i,0]
181	v1 = self.mesh.triangles[i,1]
182	v2 = self.mesh.triangles[i,2]
183
184	# Get the three vertex_points
185	xi0 = self.mesh.get_vertex_coordinate(i, 0)
186	xi1 = self.mesh.get_vertex_coordinate(i, 1)
187	xi2 = self.mesh.get_vertex_coordinate(i, 2)
188
189	# Compute gradients for each vertex
190	a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
191	1, 0, 0)
192
193	a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
194	0, 1, 0)
195
196	a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
197	0, 0, 1)
198
199	# Compute diagonal contributions
200	self.D[v0,v0] += (a0a0 + b0b0)*area
201	self.D[v1,v1] += (a1a1 + b1b1)*area
202	self.D[v2,v2] += (a2a2 + b2b2)*area
203
204	# Compute contributions for basis functions sharing edges
205	e01 = (a0a1 + b0b1)*area
206	self.D[v0,v1] += e01
207	self.D[v1,v0] += e01
208
209	e12 = (a1a2 + b1b2)*area
210	self.D[v1,v2] += e12
211	self.D[v2,v1] += e12
212
213	e20 = (a2a0 + b2b0)*area
214	self.D[v2,v0] += e20
215	self.D[v0,v2] += e20
216
217	def get_D(self):
218	return self.D.todense()
219
220
221
222	def _build_matrix_AtA_Atz(self,
223	point_coordinates,
224	z,
225	verbose = False):
226	"""Build:
227	AtA m x m interpolation matrix, and,
228	Atz m x a interpolation matrix where,
229	m is the number of basis functions phi_k (one per vertex)
230	a is the number of data attributes
231
232	This algorithm uses a quad tree data structure for fast binning of
233	data points.
234
235	If Ata is None, the matrices AtA and Atz are created.
236
237	This function can be called again and again, with sub-sets of
238	the point coordinates. Call fit to get the results.
239
240	Preconditions
241	z and points are numeric
242	Point_coordindates and mesh vertices have the same origin.
243
244	The number of attributes of the data points does not change
245	"""
246
247	# Build n x m interpolation matrix
248	if self.AtA == None:
249	# AtA and Atz need to be initialised.
250	m = self.mesh.number_of_nodes
251	if len(z.shape) > 1:
252	att_num = z.shape[1]
253	self.Atz = num.zeros((m,att_num), num.float)
254	else:
255	att_num = 1
256	self.Atz = num.zeros((m,), num.float)
257	assert z.shape[0] == point_coordinates.shape[0]
258
259	AtA = Sparse(m,m)
260	# The memory damage has been done by now.
261	else:
262	AtA = self.AtA #Did this for speed, did ~nothing
263	self.point_count += point_coordinates.shape[0]
264
265
266	inside_indices = inside_polygon(point_coordinates,
267	self.mesh_boundary_polygon,
268	closed=True,
269	verbose=False) # Too much output if True
270
271
272	n = len(inside_indices)
273
274	# Compute matrix elements for points inside the mesh
275	triangles = self.mesh.triangles # Shorthand
276	for d, i in enumerate(inside_indices):
277	# For each data_coordinate point
278	# if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
279	x = point_coordinates[i]
280
281	element_found, sigma0, sigma1, sigma2, k = \
282	search_tree_of_vertices(self.root, x)
283
284	if element_found is True:
285	j0 = triangles[k,0] #Global vertex id for sigma0
286	j1 = triangles[k,1] #Global vertex id for sigma1
287	j2 = triangles[k,2] #Global vertex id for sigma2
288
289	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
290	js = [j0,j1,j2]
291
292	for j in js:
293	self.Atz[j] += sigmas[j]*z[i]
294	#print "self.Atz building", self.Atz
295	#print "self.Atz[j]", self.Atz[j]
296	#print " sigmas[j]", sigmas[j]
297	#print "z[i]",z[i]
298	#print "result", sigmas[j]*z[i]
299
300	for k in js:
301	AtA[j,k] += sigmas[j]*sigmas[k]
302	else:
303	flag = is_inside_polygon(x,
304	self.mesh_boundary_polygon,
305	closed=True,
306	verbose=False) # Too much output if True
307	msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
308	assert flag is True, msg
309
310	# FIXME(Ole): This is the message referred to in ticket:314
311	minx = min(self.mesh_boundary_polygon[:,0])
312	maxx = max(self.mesh_boundary_polygon[:,0])
313	miny = min(self.mesh_boundary_polygon[:,1])
314	maxy = max(self.mesh_boundary_polygon[:,1])
315	msg = 'Could not find triangle for point %s. ' % str(x)
316	msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
317	% (minx, maxx, miny, maxy)
318	raise RuntimeError, msg
319
320	self.AtA = AtA
321
322
323	def fit(self, point_coordinates_or_filename=None, z=None,
324	verbose=False,
325	point_origin=None,
326	attribute_name=None,
327	max_read_lines=500):
328	"""Fit a smooth surface to given 1d array of data points z.
329
330	The smooth surface is computed at each vertex in the underlying
331	mesh using the formula given in the module doc string.
332
333	Inputs:
334	point_coordinates: The co-ordinates of the data points.
335	List of coordinate pairs [x, y] of
336	data points or an nx2 numeric array or a Geospatial_data object
337	or points file filename
338	z: Single 1d vector or array of data at the point_coordinates.
339
340	"""
341
342	# Use blocking to load in the point info
343	if type(point_coordinates_or_filename) == types.StringType:
344	msg = "Don't set a point origin when reading from a file"
345	assert point_origin is None, msg
346	filename = point_coordinates_or_filename
347
348	G_data = Geospatial_data(filename,
349	max_read_lines=max_read_lines,
350	load_file_now=False,
351	verbose=verbose)
352
353	for i, geo_block in enumerate(G_data):
354	if verbose is True and 0 == i%200:
355	# The time this will take
356	# is dependant on the # of Triangles
357
358	print 'Processing Block %d' %i
359	# FIXME (Ole): It would be good to say how many blocks
360	# there are here. But this is no longer necessary
361	# for pts files as they are reported in geospatial_data
362	# I suggest deleting this verbose output and make
363	# Geospatial_data more informative for txt files.
364	#
365	# I still think so (12/12/7, Ole).
366
367
368
369	# Build the array
370
371	points = geo_block.get_data_points(absolute=True)
372	z = geo_block.get_attributes(attribute_name=attribute_name)
373	self.build_fit_subset(points, z, verbose=verbose)
374
375	# FIXME(Ole): I thought this test would make sense here
376	# See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py
377	# Committed 11 March 2009
378	msg = 'Matrix AtA was not built'
379	assert self.AtA is not None, msg
380
381	#print 'Matrix was built OK'
382
383
384	point_coordinates = None
385	else:
386	point_coordinates = point_coordinates_or_filename
387
388	if point_coordinates is None:
389	if verbose: print 'Warning: no data points in fit'
390	msg = 'No interpolation matrix.'
391	assert self.AtA is not None, msg
392	assert self.Atz is not None
393
394	# FIXME (DSG) - do a message
395	else:
396	point_coordinates = ensure_absolute(point_coordinates,
397	geo_reference=point_origin)
398	# if isinstance(point_coordinates,Geospatial_data) and z is None:
399	# z will come from the geo-ref
400	self.build_fit_subset(point_coordinates, z, verbose)
401
402	# Check sanity
403	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
404	n = self.point_count
405	if n<m and self.alpha == 0.0:
406	msg = 'ERROR (least_squares): Too few data points\n'
407	msg += 'There are only %d data points and alpha == 0. ' %n
408	msg += 'Need at least %d\n' %m
409	msg += 'Alternatively, set smoothing parameter alpha to a small '
410	msg += 'positive value,\ne.g. 1.0e-3.'
411	raise TooFewPointsError(msg)
412
413	self._build_coefficient_matrix_B(verbose)
414	loners = self.mesh.get_lone_vertices()
415	# FIXME - make this as error message.
416	# test with
417	# Not_yet_test_smooth_att_to_mesh_with_excess_verts.
418	if len(loners)>0:
419	msg = 'WARNING: (least_squares): \nVertices with no triangles\n'
420	msg += 'All vertices should be part of a triangle.\n'
421	msg += 'In the future this will be inforced.\n'
422	msg += 'The following vertices are not part of a triangle;\n'
423	msg += str(loners)
424	print msg
425	#raise VertsWithNoTrianglesError(msg)
426
427
428	return conjugate_gradient(self.B, self.Atz, self.Atz,
429	imax=2*len(self.Atz) )
430
431
432	def build_fit_subset(self, point_coordinates, z=None, attribute_name=None,
433	verbose=False):
434	"""Fit a smooth surface to given 1d array of data points z.
435
436	The smooth surface is computed at each vertex in the underlying
437	mesh using the formula given in the module doc string.
438
439	Inputs:
440	point_coordinates: The co-ordinates of the data points.
441	List of coordinate pairs [x, y] of
442	data points or an nx2 numeric array or a Geospatial_data object
443	z: Single 1d vector or array of data at the point_coordinates.
444	attribute_name: Used to get the z values from the
445	geospatial object if no attribute_name is specified,
446	it's a bit of a lucky dip as to what attributes you get.
447	If there is only one attribute it will be that one.
448
449	"""
450
451	# FIXME(DSG-DSG): Check that the vert and point coords
452	# have the same zone.
453	if isinstance(point_coordinates,Geospatial_data):
454	point_coordinates = point_coordinates.get_data_points( \
455	absolute = True)
456
457	# Convert input to numeric arrays
458	if z is not None:
459	z = ensure_numeric(z, num.float)
460	else:
461	msg = 'z not specified'
462	assert isinstance(point_coordinates,Geospatial_data), msg
463	z = point_coordinates.get_attributes(attribute_name)
464
465	point_coordinates = ensure_numeric(point_coordinates, num.float)
466	self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
467
468
469	############################################################################
470
471	def fit_to_mesh(point_coordinates, # this can also be a points file name
472	vertex_coordinates=None,
473	triangles=None,
474	mesh=None,
475	point_attributes=None,
476	alpha=DEFAULT_ALPHA,
477	verbose=False,
478	mesh_origin=None,
479	data_origin=None,
480	max_read_lines=None,
481	attribute_name=None,
482	use_cache=False):
483	"""Wrapper around internal function _fit_to_mesh for use with caching.
484
485	"""
486
487	args = (point_coordinates, )
488	kwargs = {'vertex_coordinates': vertex_coordinates,
489	'triangles': triangles,
490	'mesh': mesh,
491	'point_attributes': point_attributes,
492	'alpha': alpha,
493	'verbose': verbose,
494	'mesh_origin': mesh_origin,
495	'data_origin': data_origin,
496	'max_read_lines': max_read_lines,
497	'attribute_name': attribute_name
498	}
499
500	if use_cache is True:
501	if isinstance(point_coordinates, basestring):
502	# We assume that point_coordinates is the name of a .csv/.txt
503	# file which must be passed onto caching as a dependency
504	# (in case it has changed on disk)
505	dep = [point_coordinates]
506	else:
507	dep = None
508
509
510	#from caching import myhash
511	#import copy
512	#print args
513	#print kwargs
514	#print 'hashing:'
515	#print 'args', myhash( (args, kwargs) )
516	#print 'again', myhash( copy.deepcopy( (args, kwargs)) )
517
518	#print 'mesh hash', myhash( kwargs['mesh'] )
519
520	#print '-------------------------'
521	#print 'vertices hash', myhash( kwargs['mesh'].nodes )
522	#print 'triangles hash', myhash( kwargs['mesh'].triangles )
523	#print '-------------------------'
524
525	#for key in mesh.__dict__:
526	# print key, myhash(mesh.__dict__[key])
527
528	#for key in mesh.quantities.keys():
529	# print key, myhash(mesh.quantities[key])
530
531	#import sys; sys.exit()
532
533	return cache(_fit_to_mesh,
534	args, kwargs,
535	verbose=verbose,
536	compression=False,
537	dependencies=dep)
538	else:
539	return apply(_fit_to_mesh,
540	args, kwargs)
541
542	def _fit_to_mesh(point_coordinates, # this can also be a points file name
543	vertex_coordinates=None,
544	triangles=None,
545	mesh=None,
546	point_attributes=None,
547	alpha=DEFAULT_ALPHA,
548	verbose=False,
549	mesh_origin=None,
550	data_origin=None,
551	max_read_lines=None,
552	attribute_name=None):
553	"""
554	Fit a smooth surface to a triangulation,
555	given data points with attributes.
556
557
558	Inputs:
559	vertex_coordinates: List of coordinate pairs [xi, eta] of
560	points constituting a mesh (or an m x 2 numeric array or
561	a geospatial object)
562	Points may appear multiple times
563	(e.g. if vertices have discontinuities)
564
565	triangles: List of 3-tuples (or a numeric array) of
566	integers representing indices of all vertices in the mesh.
567
568	point_coordinates: List of coordinate pairs [x, y] of data points
569	(or an nx2 numeric array). This can also be a .csv/.txt/.pts
570	file name.
571
572	alpha: Smoothing parameter.
573
574	mesh_origin: A geo_reference object or 3-tuples consisting of
575	UTM zone, easting and northing.
576	If specified vertex coordinates are assumed to be
577	relative to their respective origins.
578
579	point_attributes: Vector or array of data at the
580	point_coordinates.
581
582	"""
583
584	if mesh is None:
585	# FIXME(DSG): Throw errors if triangles or vertex_coordinates
586	# are None
587
588	#Convert input to numeric arrays
589	triangles = ensure_numeric(triangles, num.int)
590	vertex_coordinates = ensure_absolute(vertex_coordinates,
591	geo_reference = mesh_origin)
592
593	if verbose: print 'FitInterpolate: Building mesh'
594	mesh = Mesh(vertex_coordinates, triangles)
595	mesh.check_integrity()
596
597
598	interp = Fit(mesh=mesh,
599	verbose=verbose,
600	alpha=alpha)
601
602	vertex_attributes = interp.fit(point_coordinates,
603	point_attributes,
604	point_origin=data_origin,
605	max_read_lines=max_read_lines,
606	attribute_name=attribute_name,
607	verbose=verbose)
608
609
610	# Add the value checking stuff that's in least squares.
611	# Maybe this stuff should get pushed down into Fit.
612	# at least be a method of Fit.
613	# Or intigrate it into the fit method, saving teh max and min's
614	# as att's.
615
616	return vertex_attributes
617
618
619	#def _fit(args, *kwargs):
620	# """Private function for use with caching. Reason is that classes
621	# may change their byte code between runs which is annoying.
622	# """
623	#
624	# return Fit(args, *kwargs)
625
626
627	def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
628	alpha=DEFAULT_ALPHA, verbose= False,
629	expand_search = False,
630	precrop = False,
631	display_errors = True):
632	"""
633	Given a mesh file (tsh) and a point attribute file, fit
634	point attributes to the mesh and write a mesh file with the
635	results.
636
637	Note: the points file needs titles. If you want anuga to use the tsh file,
638	make sure the title is elevation.
639
640	NOTE: Throws IOErrors, for a variety of file problems.
641
642	"""
643
644	from load_mesh.loadASCII import import_mesh_file, \
645	export_mesh_file, concatinate_attributelist
646
647
648	try:
649	mesh_dict = import_mesh_file(mesh_file)
650	except IOError,e:
651	if display_errors:
652	print "Could not load bad file. ", e
653	raise IOError #Could not load bad mesh file.
654
655	vertex_coordinates = mesh_dict['vertices']
656	triangles = mesh_dict['triangles']
657	if isinstance(mesh_dict['vertex_attributes'], num.ndarray):
658	old_point_attributes = mesh_dict['vertex_attributes'].tolist()
659	else:
660	old_point_attributes = mesh_dict['vertex_attributes']
661
662	if isinstance(mesh_dict['vertex_attribute_titles'], num.ndarray):
663	old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
664	else:
665	old_title_list = mesh_dict['vertex_attribute_titles']
666
667	if verbose: print 'tsh file %s loaded' %mesh_file
668
669	# load in the points file
670	try:
671	geo = Geospatial_data(point_file, verbose=verbose)
672	except IOError,e:
673	if display_errors:
674	print "Could not load bad file. ", e
675	raise IOError #Re-raise exception
676
677	point_coordinates = geo.get_data_points(absolute=True)
678	title_list,point_attributes = concatinate_attributelist( \
679	geo.get_all_attributes())
680
681	if mesh_dict.has_key('geo_reference') and \
682	not mesh_dict['geo_reference'] is None:
683	mesh_origin = mesh_dict['geo_reference'].get_origin()
684	else:
685	mesh_origin = None
686
687	if verbose: print "points file loaded"
688	if verbose: print "fitting to mesh"
689	f = fit_to_mesh(point_coordinates,
690	vertex_coordinates,
691	triangles,
692	None,
693	point_attributes,
694	alpha = alpha,
695	verbose = verbose,
696	data_origin = None,
697	mesh_origin = mesh_origin)
698	if verbose: print "finished fitting to mesh"
699
700	# convert array to list of lists
701	new_point_attributes = f.tolist()
702	#FIXME have this overwrite attributes with the same title - DSG
703	#Put the newer attributes last
704	if old_title_list <> []:
705	old_title_list.extend(title_list)
706	#FIXME can this be done a faster way? - DSG
707	for i in range(len(old_point_attributes)):
708	old_point_attributes[i].extend(new_point_attributes[i])
709	mesh_dict['vertex_attributes'] = old_point_attributes
710	mesh_dict['vertex_attribute_titles'] = old_title_list
711	else:
712	mesh_dict['vertex_attributes'] = new_point_attributes
713	mesh_dict['vertex_attribute_titles'] = title_list
714
715	if verbose: print "exporting to file ", mesh_output_file
716
717	try:
718	export_mesh_file(mesh_output_file, mesh_dict)
719	except IOError,e:
720	if display_errors:
721	print "Could not write file. ", e
722	raise IOError

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

Download in other formats: