Context Navigation

source: inundation/fit_interpolate/interpolate.py @ 2852

Last change on this file since 2852 was 2852, checked in by ole, 19 years ago
Nailed the problem with tests failing when fixing duplicate timestep issue.
File size: 25.1 KB

Line
1	"""Least squares interpolation.
2
3	Implements a least-squares interpolation.
4
5	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
6	Geoscience Australia, 2004.
7
8	DESIGN ISSUES
9	* what variables should be global?
10	- if there are no global vars functions can be moved around alot easier
11
12	* The public interface
13	__init__
14	interpolate
15	interpolate_block
16
17	"""
18
19	import time
20	import os
21	from warnings import warn
22
23	from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
24	ArrayType, allclose, take, NewAxis, arange
25
26	from caching.caching import cache
27	from pyvolution.mesh import Mesh
28	from utilities.sparse import Sparse, Sparse_CSR
29	from utilities.cg_solve import conjugate_gradient, VectorShapeError
30	from coordinate_transforms.geo_reference import Geo_reference
31	from pyvolution.quad import build_quadtree
32	from utilities.numerical_tools import ensure_numeric, mean, INF
33	from utilities.polygon import in_and_outside_polygon
34	from geospatial_data.geospatial_data import Geospatial_data
35	from fit_interpolate.search_functions import search_tree_of_vertices
36	from fit_interpolate.general_fit_interpolate import FitInterpolate
37
38
39
40	class Interpolate (FitInterpolate):
41
42	def __init__(self,
43	vertex_coordinates,
44	triangles,
45	mesh_origin=None,
46	verbose=False,
47	max_vertices_per_cell=30):
48
49
50	""" Build interpolation matrix mapping from
51	function values at vertices to function values at data points
52
53	Inputs:
54
55	vertex_coordinates: List of coordinate pairs [xi, eta] of
56	points constituting a mesh (or an m x 2 Numeric array or
57	a geospatial object)
58	Points may appear multiple times
59	(e.g. if vertices have discontinuities)
60
61	triangles: List of 3-tuples (or a Numeric array) of
62	integers representing indices of all vertices in the mesh.
63
64	mesh_origin: A geo_reference object or 3-tuples consisting of
65	UTM zone, easting and northing.
66	If specified vertex coordinates are assumed to be
67	relative to their respective origins.
68
69	max_vertices_per_cell: Number of vertices in a quad tree cell
70	at which the cell is split into 4.
71
72	Note: Don't supply a vertex coords as a geospatial object and
73	a mesh origin, since geospatial has its own mesh origin.
74	"""
75
76	# Initialise variabels
77	self._A_can_be_reused = False
78	self._point_coordinates = None
79
80	FitInterpolate.__init__(self,
81	vertex_coordinates,
82	triangles,
83	mesh_origin,
84	verbose,
85	max_vertices_per_cell)
86
87	# FIXME: What is a good start_blocking_len value?
88	def interpolate(self, f, point_coordinates = None,
89	start_blocking_len = 500000, verbose=False):
90	"""Interpolate mesh data f to determine values, z, at points.
91
92	f is the data on the mesh vertices.
93
94	The mesh values representing a smooth surface are
95	assumed to be specified in f.
96
97	Inputs:
98	f: Vector or array of data at the mesh vertices.
99	If f is an array, interpolation will be done for each column as
100	per underlying matrix-matrix multiplication
101	point_coordinates: Interpolate mesh data to these positions.
102	List of coordinate pairs [x, y] of
103	data points or an nx2 Numeric array or a Geospatial_data object
104
105	If point_coordinates is absent, the points inputted last time
106	this method was called are used, if possible.
107	start_blocking_len: If the # of points is more or greater than this,
108	start blocking
109
110	Output:
111	Interpolated values at inputted points (z).
112	"""
113	#print "point_coordinates interpolate.interpolate",point_coordinates
114	if isinstance(point_coordinates,Geospatial_data):
115	point_coordinates = point_coordinates.get_data_points( \
116	absolute = True)
117
118	# Can I interpolate, based on previous point_coordinates?
119	if point_coordinates is None:
120	if self._A_can_be_reused is True and \
121	len(self._point_coordinates) < start_blocking_len:
122	z = self._get_point_data_z(f,
123	verbose=verbose)
124	elif self._point_coordinates is not None:
125	# if verbose, give warning
126	if verbose:
127	print 'WARNING: Recalculating A matrix, due to blocking.'
128	point_coordinates = self._point_coordinates
129	else:
130	#There are no good point_coordinates. import sys; sys.exit()
131	msg = 'ERROR (interpolate.py): No point_coordinates inputted'
132	raise Exception(msg)
133
134
135	if point_coordinates is not None:
136	self._point_coordinates = point_coordinates
137	if len(point_coordinates) < start_blocking_len or \
138	start_blocking_len == 0:
139	self._A_can_be_reused = True
140	z = self.interpolate_block(f, point_coordinates,
141	verbose=verbose)
142	else:
143	#Handle blocking
144	self._A_can_be_reused = False
145	start=0
146	# creating a dummy array to concatenate to.
147
148	f = ensure_numeric(f, Float)
149	z = zeros((0,f.shape[1]))
150	for end in range(start_blocking_len
151	,len(point_coordinates)
152	,start_blocking_len):
153	t = self.interpolate_block(f, point_coordinates[start:end],
154	verbose=verbose)
155	#print "t", t
156	#print "z", z
157	z = concatenate((z,t))
158	start = end
159	end = len(point_coordinates)
160	t = self.interpolate_block(f, point_coordinates[start:end],
161	verbose=verbose)
162	z = concatenate((z,t))
163	return z
164
165	def interpolate_block(self, f, point_coordinates = None, verbose=False):
166	"""
167	Call this if you want to control the blocking or make sure blocking
168	doesn't occur.
169
170	Return the point data, z.
171
172	See interpolate for doc info.
173	"""
174	if isinstance(point_coordinates,Geospatial_data):
175	point_coordinates = point_coordinates.get_data_points( \
176	absolute = True)
177	if point_coordinates is not None:
178	self._A =self._build_interpolation_matrix_A(point_coordinates,
179	verbose=verbose)
180	return self._get_point_data_z(f)
181
182	def _get_point_data_z(self, f, verbose=False):
183	"""
184	Return the point data, z.
185
186	Precondition,
187	The _A matrix has been created
188	"""
189	z = self._A * f
190	# Taking into account points outside the mesh.
191	#print "self.outside_poly_indices", self.outside_poly_indices
192	#print "self.inside_poly_indices", self.inside_poly_indices
193	#print "z", z
194	for i in self.outside_poly_indices:
195	z[i] = INF
196	return z
197
198
199	def _build_interpolation_matrix_A(self,
200	point_coordinates,
201	verbose = False):
202	"""Build n x m interpolation matrix, where
203	n is the number of data points and
204	m is the number of basis functions phi_k (one per vertex)
205
206	This algorithm uses a quad tree data structure for fast binning
207	of data points
208	origin is a 3-tuple consisting of UTM zone, easting and northing.
209	If specified coordinates are assumed to be relative to this origin.
210
211	This one will override any data_origin that may be specified in
212	instance interpolation
213
214	Preconditions
215	Point_coordindates and mesh vertices have the same origin.
216	"""
217
218
219
220	#Convert point_coordinates to Numeric arrays, in case it was a list.
221	point_coordinates = ensure_numeric(point_coordinates, Float)
222
223	if verbose: print 'Getting indices inside mesh boundary'
224	#print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon()
225	self.inside_poly_indices, self.outside_poly_indices = \
226	in_and_outside_polygon(point_coordinates,
227	self.mesh.get_boundary_polygon(),
228	closed = True, verbose = verbose)
229	#print "self.inside_poly_indices",self.inside_poly_indices
230	#print "self.outside_poly_indices",self.outside_poly_indices
231	#Build n x m interpolation matrix
232	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
233	n = point_coordinates.shape[0] #Nbr of data points
234
235	if verbose: print 'Number of datapoints: %d' %n
236	if verbose: print 'Number of basis functions: %d' %m
237
238	A = Sparse(n,m)
239
240	#Compute matrix elements for points inside the mesh
241	for i in self.inside_poly_indices:
242	#For each data_coordinate point
243	if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
244	x = point_coordinates[i]
245	element_found, sigma0, sigma1, sigma2, k = \
246	search_tree_of_vertices(self.root, self.mesh, x)
247	#Update interpolation matrix A if necessary
248	if element_found is True:
249	#Assign values to matrix A
250
251	j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
252	j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
253	j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
254
255	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
256	js = [j0,j1,j2]
257
258	for j in js:
259	A[i,j] = sigmas[j]
260	else:
261	msg = 'Could not find triangle for point', x
262	raise Exception(msg)
263	return A
264
265	class Interpolation_function:
266	"""Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
267	which is interpolated from time series defined at vertices of
268	triangular mesh (such as those stored in sww files)
269
270	Let m be the number of vertices, n the number of triangles
271	and p the number of timesteps.
272
273	Mandatory input
274	time: px1 array of monotonously increasing times (Float)
275	quantities: Dictionary of arrays or 1 array (Float)
276	The arrays must either have dimensions pxm or mx1.
277	The resulting function will be time dependent in
278	the former case while it will be constan with
279	respect to time in the latter case.
280
281	Optional input:
282	quantity_names: List of keys into the quantities dictionary
283	vertex_coordinates: mx2 array of coordinates (Float)
284	triangles: nx3 array of indices into vertex_coordinates (Int)
285	interpolation_points: Nx2 array of coordinates to be interpolated to
286	verbose: Level of reporting
287
288
289	The quantities returned by the callable object are specified by
290	the list quantities which must contain the names of the
291	quantities to be returned and also reflect the order, e.g. for
292	the shallow water wave equation, on would have
293	quantities = ['stage', 'xmomentum', 'ymomentum']
294
295	The parameter interpolation_points decides at which points interpolated
296	quantities are to be computed whenever object is called.
297	If None, return average value
298	"""
299
300
301
302	def __init__(self,
303	time,
304	quantities,
305	quantity_names = None,
306	vertex_coordinates = None,
307	triangles = None,
308	interpolation_points = None,
309	verbose = False):
310	"""Initialise object and build spatial interpolation if required
311	"""
312
313	from Numeric import array, zeros, Float, alltrue, concatenate,\
314	reshape, ArrayType
315
316
317	#from util import mean, ensure_numeric
318	from config import time_format
319	import types
320
321
322	#Check temporal info
323	time = ensure_numeric(time)
324	msg = 'Time must be a monotonuosly '
325	msg += 'increasing sequence %s' %time
326	assert alltrue(time[1:] - time[:-1] >= 0 ), msg
327
328
329	#Check if quantities is a single array only
330	if type(quantities) != types.DictType:
331	quantities = ensure_numeric(quantities)
332	quantity_names = ['Attribute']
333
334	#Make it a dictionary
335	quantities = {quantity_names[0]: quantities}
336
337
338	#Use keys if no names are specified
339	if quantity_names is None:
340	quantity_names = quantities.keys()
341
342
343	#Check spatial info
344	if vertex_coordinates is None:
345	self.spatial = False
346	else:
347	vertex_coordinates = ensure_numeric(vertex_coordinates)
348
349	assert triangles is not None, 'Triangles array must be specified'
350	triangles = ensure_numeric(triangles)
351	self.spatial = True
352
353
354	#Save for use with statistics
355	self.quantity_names = quantity_names
356	self.quantities = quantities
357	self.vertex_coordinates = vertex_coordinates
358	self.interpolation_points = interpolation_points
359	self.T = time[:] # Time assumed to be relative to starttime
360	self.index = 0 # Initial time index
361	self.precomputed_values = {}
362
363	#Precomputed spatial interpolation if requested
364	if interpolation_points is not None:
365	if self.spatial is False:
366	raise 'Triangles and vertex_coordinates must be specified'
367
368	try:
369	self.interpolation_points = ensure_numeric(interpolation_points)
370	except:
371	msg = 'Interpolation points must be an N x 2 Numeric array '+\
372	'or a list of points\n'
373	msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
374	'...')
375	raise msg
376
377
378	m = len(self.interpolation_points)
379	p = len(self.T)
380
381	for name in quantity_names:
382	self.precomputed_values[name] = zeros((p, m), Float)
383
384	#Build interpolator
385	interpol = Interpolate(vertex_coordinates,
386	triangles,
387	#point_coordinates = \
388	#self.interpolation_points,
389	#alpha = 0,
390	verbose = verbose)
391
392	if verbose: print 'Interpolate'
393	for i, t in enumerate(self.T):
394	#Interpolate quantities at this timestep
395	if verbose and i%((p+10)/10)==0:
396	print ' time step %d of %d' %(i, p)
397
398	for name in quantity_names:
399	if len(quantities[name].shape) == 2:
400	result = interpol.interpolate(quantities[name][i,:],
401	point_coordinates = \
402	self.interpolation_points)
403	else:
404	#Assume no time dependency
405	result = interpol.interpolate(quantities[name][:],
406	point_coordinates = \
407	self.interpolation_points)
408
409	self.precomputed_values[name][i, :] = result
410
411	#Report
412	if verbose:
413	print self.statistics()
414	#self.print_statistics()
415
416	else:
417	#Store quantitites as is
418	for name in quantity_names:
419	self.precomputed_values[name] = quantities[name]
420
421
422	def __repr__(self):
423	#return 'Interpolation function (spatio-temporal)'
424	return self.statistics()
425
426
427	def __call__(self, t, point_id = None, x = None, y = None):
428	"""Evaluate f(t), f(t, point_id) or f(t, x, y)
429
430	Inputs:
431	t: time - Model time. Must lie within existing timesteps
432	point_id: index of one of the preprocessed points.
433	x, y: Overrides location, point_id ignored
434
435	If spatial info is present and all of x,y,point_id
436	are None an exception is raised
437
438	If no spatial info is present, point_id and x,y arguments are ignored
439	making f a function of time only.
440
441
442	FIXME: point_id could also be a slice
443	FIXME: What if x and y are vectors?
444	FIXME: What about f(x,y) without t?
445	"""
446
447	from math import pi, cos, sin, sqrt
448	from Numeric import zeros, Float
449	from util import mean
450
451	if self.spatial is True:
452	if point_id is None:
453	if x is None or y is None:
454	msg = 'Either point_id or x and y must be specified'
455	raise Exception(msg)
456	else:
457	if self.interpolation_points is None:
458	msg = 'Interpolation_function must be instantiated ' +\
459	'with a list of interpolation points before parameter ' +\
460	'point_id can be used'
461	raise Exception(msg)
462
463	msg = 'Time interval [%.16f:%.16f]' %(self.T[0], self.T[-1])
464	msg += ' does not match model time: %.16f\n' %t
465	if t < self.T[0]: raise Exception(msg)
466	if t > self.T[-1]: raise Exception(msg)
467
468	oldindex = self.index #Time index
469
470	#Find current time slot
471	while t > self.T[self.index]: self.index += 1
472	while t < self.T[self.index]: self.index -= 1
473
474	if t == self.T[self.index]:
475	#Protect against case where t == T[-1] (last time)
476	# - also works in general when t == T[i]
477	ratio = 0
478	else:
479	#t is now between index and index+1
480	ratio = (t - self.T[self.index])/\
481	(self.T[self.index+1] - self.T[self.index])
482
483	#Compute interpolated values
484	q = zeros(len(self.quantity_names), Float)
485	#print "self.precomputed_values", self.precomputed_values
486	for i, name in enumerate(self.quantity_names):
487	Q = self.precomputed_values[name]
488
489	if self.spatial is False:
490	#If there is no spatial info
491	assert len(Q.shape) == 1
492
493	Q0 = Q[self.index]
494	if ratio > 0: Q1 = Q[self.index+1]
495
496	else:
497	if x is not None and y is not None:
498	#Interpolate to x, y
499
500	raise 'x,y interpolation not yet implemented'
501	else:
502	#Use precomputed point
503	Q0 = Q[self.index, point_id]
504	if ratio > 0:
505	Q1 = Q[self.index+1, point_id]
506
507	#Linear temporal interpolation
508	if ratio > 0:
509	if Q0 == INF and Q1 == INF:
510	q[i] = Q0
511	else:
512	q[i] = Q0 + ratio*(Q1 - Q0)
513	else:
514	q[i] = Q0
515
516
517	#Return vector of interpolated values
518	#if len(q) == 1:
519	# return q[0]
520	#else:
521	# return q
522
523
524	#Return vector of interpolated values
525	#FIXME:
526	if self.spatial is True:
527	return q
528	else:
529	#Replicate q according to x and y
530	#This is e.g used for Wind_stress
531	if x is None or y is None:
532	return q
533	else:
534	try:
535	N = len(x)
536	except:
537	return q
538	else:
539	from Numeric import ones, Float
540	#x is a vector - Create one constant column for each value
541	N = len(x)
542	assert len(y) == N, 'x and y must have same length'
543	res = []
544	for col in q:
545	res.append(col*ones(N, Float))
546
547	return res
548
549
550	def statistics(self):
551	"""Output statistics about interpolation_function
552	"""
553
554	vertex_coordinates = self.vertex_coordinates
555	interpolation_points = self.interpolation_points
556	quantity_names = self.quantity_names
557	quantities = self.quantities
558	precomputed_values = self.precomputed_values
559
560	x = vertex_coordinates[:,0]
561	y = vertex_coordinates[:,1]
562
563	str = '------------------------------------------------\n'
564	str += 'Interpolation_function (spatio-temporal) statistics:\n'
565	str += ' Extent:\n'
566	str += ' x in [%f, %f], len(x) == %d\n'\
567	%(min(x), max(x), len(x))
568	str += ' y in [%f, %f], len(y) == %d\n'\
569	%(min(y), max(y), len(y))
570	str += ' t in [%f, %f], len(t) == %d\n'\
571	%(min(self.T), max(self.T), len(self.T))
572	str += ' Quantities:\n'
573	for name in quantity_names:
574	q = quantities[name][:].flat
575	str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
576
577	if interpolation_points is not None:
578	str += ' Interpolation points (xi, eta):'\
579	' number of points == %d\n' %interpolation_points.shape[0]
580	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
581	max(interpolation_points[:,0]))
582	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
583	max(interpolation_points[:,1]))
584	str += ' Interpolated quantities (over all timesteps):\n'
585
586	for name in quantity_names:
587	q = precomputed_values[name][:].flat
588	str += ' %s at interpolation points in [%f, %f]\n'\
589	%(name, min(q), max(q))
590	str += '------------------------------------------------\n'
591
592	return str
593
594	def interpolate_sww(sww_file, time, interpolation_points,
595	quantity_names = None, verbose = False):
596	"""
597	obsolete.
598	use file_function in utils
599	"""
600	#open sww file
601	x, y, volumes, time, quantities = read_sww(sww_file)
602	print "x",x
603	print "y",y
604
605	print "time", time
606	print "quantities", quantities
607
608	#Add the x and y together
609	vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1)
610
611	#Will return the quantity values at the specified times and locations
612	interp = Interpolation_interface(
613	time,
614	quantities,
615	quantity_names = quantity_names,
616	vertex_coordinates = vertex_coordinates,
617	triangles = volumes,
618	interpolation_points = interpolation_points,
619	verbose = verbose)
620
621
622	def read_sww(file_name):
623	"""
624	obsolete - Nothing should be calling this
625
626	Read in an sww file.
627
628	Input;
629	file_name - the sww file
630
631	Output;
632	x - Vector of x values
633	y - Vector of y values
634	z - Vector of bed elevation
635	volumes - Array. Each row has 3 values, representing
636	the vertices that define the volume
637	time - Vector of the times where there is stage information
638	stage - array with respect to time and vertices (x,y)
639	"""
640
641	#FIXME Have this reader as part of data_manager?
642
643	from Scientific.IO.NetCDF import NetCDFFile
644	import tempfile
645	import sys
646	import os
647
648	#Check contents
649	#Get NetCDF
650
651	# see if the file is there. Throw a QUIET IO error if it isn't
652	# I don't think I could implement the above
653
654	#throws prints to screen if file not present
655	junk = tempfile.mktemp(".txt")
656	fd = open(junk,'w')
657	stdout = sys.stdout
658	sys.stdout = fd
659	fid = NetCDFFile(file_name, 'r')
660	sys.stdout = stdout
661	fd.close()
662	#clean up
663	os.remove(junk)
664
665	# Get the variables
666	x = fid.variables['x'][:]
667	y = fid.variables['y'][:]
668	volumes = fid.variables['volumes'][:]
669	time = fid.variables['time'][:]
670
671	keys = fid.variables.keys()
672	keys.remove("x")
673	keys.remove("y")
674	keys.remove("volumes")
675	keys.remove("time")
676	#Turn NetCDF objects into Numeric arrays
677	quantities = {}
678	for name in keys:
679	quantities[name] = fid.variables[name][:]
680
681	fid.close()
682	return x, y, volumes, time, quantities
683
684	#-------------------------------------------------------------
685	if __name__ == "__main__":
686	a = [0.0, 0.0]
687	b = [0.0, 2.0]
688	c = [2.0,0.0]
689	points = [a, b, c]
690	vertices = [ [1,0,2] ] #bac
691
692	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
693
694	interp = Interpolate(points, vertices) #, data)
695	A = interp._build_interpolation_matrix_A(data, verbose=True)
696	A = A.todense()
697	print "A",A
698	assert allclose(A, [[1./3, 1./3, 1./3]])
699	print "finished"

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

Download in other formats: