Context Navigation

source: inundation/pyvolution/least_squares.py @ 3072

Last change on this file since 3072 was 2884, checked in by ole, 19 years ago
Replaced reference to f.T with f.get_time() and retested
File size: 52.2 KB

Rev	Line
[1911]	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
[2663]	20	import exceptions
	21	class ShapeError(exceptions.Exception): pass
	22	class FittingError(exceptions.Exception): pass
[1911]	23
[2663]	24
[1911]	25	#from general_mesh import General_mesh
[2778]	26	from Numeric import zeros, array, Float, Int, transpose, concatenate, ArrayType, NewAxis
[2502]	27	from pyvolution.mesh import Mesh
[1911]	28
[2778]	29	from Numeric import dot, zeros, take, compress, array, Float, Int, transpose, concatenate, ArrayType
[2503]	30	from utilities.sparse import Sparse, Sparse_CSR
	31	from utilities.cg_solve import conjugate_gradient, VectorShapeError
[2516]	32	from utilities.numerical_tools import ensure_numeric, mean, gradient
[1911]	33
[2516]	34
[1911]	35	from coordinate_transforms.geo_reference import Geo_reference
	36
	37	import time
	38
	39	DEFAULT_ALPHA = 0.001
	40
	41	def 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"
[2447]	113	if verbose: print "fitting to mesh"
[1911]	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
[2447]	142	if verbose: print "exporting to file ", mesh_output_file
[1911]	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
	151	def fit_to_mesh(vertex_coordinates,
	152	triangles,
	153	point_coordinates,
	154	point_attributes,
	155	alpha = DEFAULT_ALPHA,
	156	verbose = False,
[2668]	157	acceptable_overshoot = 1.01,
[1911]	158	expand_search = False,
	159	data_origin = None,
	160	mesh_origin = None,
[2447]	161	precrop = False,
	162	use_cache = False):
[1911]	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
[2666]	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
[1911]	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	"""
[2447]	195	if use_cache is True:
	196	from 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)
[2347]	219
[1911]	220	vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
[2663]	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]
[2665]	230	Z = ensure_numeric(point_attributes)
	231	if len(Z.shape) == 1:
	232	Z = Z[:, NewAxis]
	233
[2663]	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
[2583]	245
[2663]	246	#Vertex coordinates
	247	Xi = vertex_coordinates[:,0]
	248	Eta = vertex_coordinates[:,1]
[2665]	249	Zeta = ensure_numeric(vertex_attributes)
	250	if len(Zeta.shape) == 1:
	251	Zeta = Zeta[:, NewAxis]
[2663]	252
[2665]	253	for i in range(Zeta.shape[1]): #For each attribute
	254	zeta = Zeta[:,i]
	255	z = Z[:,i]
[2666]	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
[2665]	263
[2666]	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,
[2665]	269	min(zeta), max(zeta))
[2666]	270	msg += 'If greater range is needed, increase the value of '
[2668]	271	msg += 'acceptable_fit_overshoot (currently %.2f).\n' %(acceptable_overshoot)
[2663]	272
	273
[2666]	274	offending_vertices = (zeta > upper_limit or zeta < lower_limit)
[2665]	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)
[2663]	280
[2665]	281	msg += 'Offending locations:\n %s' %(offending_coordinates)
	282
	283	raise FittingError, msg
[2663]	284
	285
	286
[2665]	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))
[2585]	295	print
[2663]	296
[2665]	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
[2583]	304
[2665]	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
[1911]	311	return vertex_attributes
	312
	313
	314
	315	def 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
[2447]	370	def _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)
[1911]	376
[2447]	377
[1911]	378	class 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
[1979]	433	if mesh_origin is None:
[1911]	434	geo = None
	435	else:
	436	geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
[2778]	437
	438
[2502]	439	self.mesh = Mesh(vertex_coordinates, triangles,
	440	geo_reference = geo)
[2683]	441
	442	if verbose: print 'Checking mesh integrity'
[1911]	443	self.mesh.check_integrity()
[2683]	444
	445	if verbose: print 'Mesh integrity checked'
	446
[1911]	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)
[2683]	477	if verbose: print 'Finished interpolation'
[1911]	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
[2347]	535
[1911]	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
[1941]	543	from pyvolution.quad import build_quadtree
[2583]	544	from utilities.polygon import inside_polygon
[1911]	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
[2347]	598
	599
[1911]	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)
[1941]	642	#root.show()
	643	self.expanded_quad_searches = []
[1911]	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
[2778]	655
[1911]	656	element_found, sigma0, sigma1, sigma2, k = \
[2778]	657	self.search_triangles_of_vertices(candidate_vertices, x)
[1941]	658	first_expansion = True
[1911]	659	while not element_found and is_more_elements and expand_search:
[1975]	660	#if verbose: print 'Expanding search'
[1942]	661	if first_expansion == True:
	662	self.expanded_quad_searches.append(1)
[1944]	663	first_expansion = False
[1941]	664	else:
[1944]	665	end = len(self.expanded_quad_searches) - 1
[1941]	666	assert end >= 0
[1942]	667	self.expanded_quad_searches[end] += 1
[1911]	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
[1941]	679
[1911]	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):
[2778]	705
	706
[1911]	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
[1941]	715	#print "$ candidate_vertices", candidate_vertices
[1911]	716	#For all vertices in same cell as point x
	717	for v in candidate_vertices:
[1975]	718	#FIXME (DSG-DSG): this catches verts with no triangle.
	719	#Currently pmesh is producing these.
	720	#this should be stopped,
[1979]	721	if self.mesh.vertexlist[v] is None:
[1975]	722	continue
[1911]	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] += (a0a0 + b0b0)*area
	922	self.D[v1,v1] += (a1a1 + b1b1)*area
	923	self.D[v2,v2] += (a2a2 + b2b2)*area
	924
	925	#Compute contributions for basis functions sharing edges
	926	e01 = (a0a1 + b0b1)*area
	927	self.D[v0,v1] += e01
	928	self.D[v1,v0] += e01
	929
	930	e12 = (a1a2 + b1b2)*area
	931	self.D[v1,v2] += e12
	932	self.D[v2,v1] += e12
	933
	934	e20 = (a2a0 + b2b0)*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	"""
[2750]	951
[1911]	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.'
[2689]	976	raise Exception(msg)
[1911]	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	"""
[2764]	1030	print "obsolete in least_squares, use fit_interpolate.interpolate"
[1911]	1031	return self.A * f
	1032
	1033	def cull_outsiders(self, f):
	1034	pass
	1035
	1036
	1037
	1038
	1039	class 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
[2633]	1052	the former case while it will be constant with
[1911]	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
[1941]	1075
[1911]	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 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
[2884]	1134	self.time = time[:] # Time assumed to be relative to starttime
[1911]	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)
[2884]	1156	p = len(self.time)
[1911]	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'
[2884]	1171	for i, t in enumerate(self.time):
[1911]	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:
[2884]	1201	# self.values[name] = zeros(len(self.time), Float)
[1911]	1202	#
	1203	# if verbose: print 'Compute mean values'
[2884]	1204	# for i, t in enumerate(self.time):
	1205	# if verbose: print ' time step %d of %d' %(i, len(self.time))
[1911]	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
[2526]	1239	from utilities.numerical_tools import mean
[1911]	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'
[2689]	1245	raise Exception(msg)
[1911]	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'
[2689]	1251	raise Exception(msg)
[1911]	1252
	1253
[2884]	1254	msg = 'Time interval [%s:%s]' %(self.time[0], self.time[-1])
[1911]	1255	msg += ' does not match model time: %s\n' %t
[2884]	1256	if t < self.time[0]: raise Exception(msg)
	1257	if t > self.time[-1]: raise Exception(msg)
[1911]	1258
	1259	oldindex = self.index #Time index
	1260
	1261	#Find current time slot
[2884]	1262	while t > self.time[self.index]: self.index += 1
	1263	while t < self.time[self.index]: self.index -= 1
[1911]	1264
[2884]	1265	if t == self.time[self.index]:
[1911]	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
[2884]	1271	ratio = (t - self.time[self.index])/\
	1272	(self.time[self.index+1] - self.time[self.index])
[1911]	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
[1979]	1318	if x is None or y is None:
[1911]	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
[2884]	1335
	1336	def get_time(self):
	1337	"""Return model time as a vector of timesteps
	1338	"""
	1339	return self.time
[1911]	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'\
[2884]	1363	%(min(self.time), max(self.time), len(self.time))
[1911]	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'\
[2884]	1395	# %(min(self.time), max(self.time), len(self.time))
[1911]	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	#-------------------------------------------------------------
	1416	if __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.

Download in other formats: