Context Navigation

source: inundation/fit_interpolate/spike_least_squares.py @ 3085

Last change on this file since 3085 was 2897, checked in by duncan, 19 years ago
close to finishing fit, which replaces least squares
File size: 50.3 KB

Rev	Line
[2695]	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
	20	import exceptions
	21	class ShapeError(exceptions.Exception): pass
	22	class FittingError(exceptions.Exception): pass
	23
	24
	25	#from general_mesh import General_mesh
	26	from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, ArrayType, NewAxis
	27	from pyvolution.mesh import Mesh
	28
	29	from Numeric import zeros, take, compress, array, Float, Int, dot, transpose, concatenate, ArrayType
	30	from utilities.sparse import Sparse, Sparse_CSR
	31	from utilities.cg_solve import conjugate_gradient, VectorShapeError
	32	from utilities.numerical_tools import ensure_numeric, mean, gradient
	33
	34
	35	from coordinate_transforms.geo_reference import Geo_reference
	36
	37	import time
	38
	39
	40
	41	DEFAULT_ALPHA = 0.001
	42
	43	def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
	44	alpha=DEFAULT_ALPHA, verbose= False,
	45	expand_search = False,
	46	data_origin = None,
	47	mesh_origin = None,
	48	precrop = False,
	49	display_errors = True):
	50	"""
	51	Given a mesh file (tsh) and a point attribute file (xya), fit
	52	point attributes to the mesh and write a mesh file with the
	53	results.
	54
	55
	56	If data_origin is not None it is assumed to be
	57	a 3-tuple with geo referenced
	58	UTM coordinates (zone, easting, northing)
	59
	60	NOTE: Throws IOErrors, for a variety of file problems.
	61
	62	mesh_origin is the same but refers to the input tsh file.
	63	FIXME: When the tsh format contains it own origin, these parameters can go.
	64	FIXME: And both origins should be obtained from the specified files.
	65	"""
	66
	67	from load_mesh.loadASCII import import_mesh_file, \
	68	import_points_file, export_mesh_file, \
	69	concatinate_attributelist
	70
	71
	72	try:
	73	mesh_dict = import_mesh_file(mesh_file)
	74	except IOError,e:
	75	if display_errors:
	76	print "Could not load bad file. ", e
	77	raise IOError #Re-raise exception
	78
	79	vertex_coordinates = mesh_dict['vertices']
	80	triangles = mesh_dict['triangles']
	81	if type(mesh_dict['vertex_attributes']) == ArrayType:
	82	old_point_attributes = mesh_dict['vertex_attributes'].tolist()
	83	else:
	84	old_point_attributes = mesh_dict['vertex_attributes']
	85
	86	if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
	87	old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
	88	else:
	89	old_title_list = mesh_dict['vertex_attribute_titles']
	90
	91	if verbose: print 'tsh file %s loaded' %mesh_file
	92
	93	# load in the .pts file
	94	try:
	95	point_dict = import_points_file(point_file, verbose=verbose)
	96	except IOError,e:
	97	if display_errors:
	98	print "Could not load bad file. ", e
	99	raise IOError #Re-raise exception
	100
	101	point_coordinates = point_dict['pointlist']
	102	title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
	103
	104	if point_dict.has_key('geo_reference') and not point_dict['geo_reference'] is None:
	105	data_origin = point_dict['geo_reference'].get_origin()
	106	else:
	107	data_origin = (56, 0, 0) #FIXME(DSG-DSG)
	108
	109	if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
	110	mesh_origin = mesh_dict['geo_reference'].get_origin()
	111	else:
	112	mesh_origin = (56, 0, 0) #FIXME(DSG-DSG)
	113
	114	if verbose: print "points file loaded"
	115	if verbose: print "fitting to mesh"
	116	f = fit_to_mesh(vertex_coordinates,
	117	triangles,
	118	point_coordinates,
	119	point_attributes,
	120	alpha = alpha,
	121	verbose = verbose,
	122	expand_search = expand_search,
	123	data_origin = data_origin,
	124	mesh_origin = mesh_origin,
	125	precrop = precrop)
	126	if verbose: print "finished fitting to mesh"
	127
	128	# convert array to list of lists
	129	new_point_attributes = f.tolist()
	130	#FIXME have this overwrite attributes with the same title - DSG
	131	#Put the newer attributes last
	132	if old_title_list <> []:
	133	old_title_list.extend(title_list)
	134	#FIXME can this be done a faster way? - DSG
	135	for i in range(len(old_point_attributes)):
	136	old_point_attributes[i].extend(new_point_attributes[i])
	137	mesh_dict['vertex_attributes'] = old_point_attributes
	138	mesh_dict['vertex_attribute_titles'] = old_title_list
	139	else:
	140	mesh_dict['vertex_attributes'] = new_point_attributes
	141	mesh_dict['vertex_attribute_titles'] = title_list
	142
	143	#FIXME (Ole): Remember to output mesh_origin as well
	144	if verbose: print "exporting to file ", mesh_output_file
	145
	146	try:
	147	export_mesh_file(mesh_output_file, mesh_dict)
	148	except IOError,e:
	149	if display_errors:
	150	print "Could not write file. ", e
	151	raise IOError
	152
	153	def fit_to_mesh(vertex_coordinates,
	154	triangles,
	155	point_coordinates,
	156	point_attributes,
	157	alpha = DEFAULT_ALPHA,
	158	verbose = False,
	159	acceptable_overshoot = 1.01,
	160	expand_search = False,
	161	data_origin = None,
	162	mesh_origin = None,
	163	precrop = False,
	164	use_cache = False):
	165	"""
	166	Fit a smooth surface to a triangulation,
	167	given data points with attributes.
	168
	169
	170	Inputs:
	171
	172	vertex_coordinates: List of coordinate pairs [xi, eta] of points
	173	constituting mesh (or a an m x 2 Numeric array)
	174
	175	triangles: List of 3-tuples (or a Numeric array) of
	176	integers representing indices of all vertices in the mesh.
	177
	178	point_coordinates: List of coordinate pairs [x, y] of data points
	179	(or an nx2 Numeric array)
	180
	181	alpha: Smoothing parameter.
	182
	183	acceptable overshoot: controls the allowed factor by which fitted values
	184	may exceed the value of input data. The lower limit is defined
	185	as min(z) - acceptable_overshoot*delta z and upper limit
	186	as max(z) + acceptable_overshoot*delta z
	187
	188
	189	point_attributes: Vector or array of data at the point_coordinates.
	190
	191	data_origin and mesh_origin are 3-tuples consisting of
	192	UTM zone, easting and northing. If specified
	193	point coordinates and vertex coordinates are assumed to be
	194	relative to their respective origins.
	195
	196	"""
	197
	198	if use_cache is True:
	199	from caching.caching import cache
	200	interp = cache(_interpolation,
	201	(vertex_coordinates,
	202	triangles,
	203	point_coordinates),
	204	{'alpha': alpha,
	205	'verbose': verbose,
	206	'expand_search': expand_search,
	207	'data_origin': data_origin,
	208	'mesh_origin': mesh_origin,
	209	'precrop': precrop},
	210	verbose = verbose)
	211
	212	else:
	213	interp = Interpolation(vertex_coordinates,
	214	triangles,
	215	point_coordinates,
	216	alpha = alpha,
	217	verbose = verbose,
	218	expand_search = expand_search,
	219	data_origin = data_origin,
	220	mesh_origin = mesh_origin,
	221	precrop = precrop)
	222
	223	vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
	224
	225
	226	#Sanity check
	227	point_coordinates = ensure_numeric(point_coordinates)
	228	vertex_coordinates = ensure_numeric(vertex_coordinates)
	229
	230	#Data points
	231	X = point_coordinates[:,0]
	232	Y = point_coordinates[:,1]
	233	Z = ensure_numeric(point_attributes)
	234	if len(Z.shape) == 1:
	235	Z = Z[:, NewAxis]
	236
	237
	238	#Data points inside mesh boundary
	239	indices = interp.point_indices
	240	if indices is not None:
	241	Xc = take(X, indices)
	242	Yc = take(Y, indices)
	243	Zc = take(Z, indices)
	244	else:
	245	Xc = X
	246	Yc = Y
	247	Zc = Z
	248
	249	#Vertex coordinates
	250	Xi = vertex_coordinates[:,0]
	251	Eta = vertex_coordinates[:,1]
	252	Zeta = ensure_numeric(vertex_attributes)
	253	if len(Zeta.shape) == 1:
	254	Zeta = Zeta[:, NewAxis]
	255
	256	for i in range(Zeta.shape[1]): #For each attribute
	257	zeta = Zeta[:,i]
	258	z = Z[:,i]
	259	zc = Zc[:,i]
	260
	261	max_zc = max(zc)
	262	min_zc = min(zc)
	263	delta_zc = max_zc-min_zc
	264	upper_limit = max_zc + delta_zc*acceptable_overshoot
	265	lower_limit = min_zc - delta_zc*acceptable_overshoot
	266
	267
	268	if max(zeta) > upper_limit or min(zeta) < lower_limit:
	269	msg = 'Least sqares produced values outside the allowed '
	270	msg += 'range [%f, %f].\n' %(lower_limit, upper_limit)
	271	msg += 'z in [%f, %f], zeta in [%f, %f].\n' %(min_zc, max_zc,
	272	min(zeta), max(zeta))
	273	msg += 'If greater range is needed, increase the value of '
	274	msg += 'acceptable_fit_overshoot (currently %.2f).\n' %(acceptable_overshoot)
	275
	276
	277	offending_vertices = (zeta > upper_limit or zeta < lower_limit)
	278	Xi_c = compress(offending_vertices, Xi)
	279	Eta_c = compress(offending_vertices, Eta)
	280	offending_coordinates = concatenate((Xi_c[:, NewAxis],
	281	Eta_c[:, NewAxis]),
	282	axis=1)
	283
	284	msg += 'Offending locations:\n %s' %(offending_coordinates)
	285
	286	raise FittingError, msg
	287
	288
	289
	290	if verbose:
	291	print '+------------------------------------------------'
	292	print 'Least squares statistics'
	293	print '+------------------------------------------------'
	294	print 'points: %d points' %(len(z))
	295	print ' x in [%f, %f]'%(min(X), max(X))
	296	print ' y in [%f, %f]'%(min(Y), max(Y))
	297	print ' z in [%f, %f]'%(min(z), max(z))
	298	print
	299
	300	if indices is not None:
	301	print 'Cropped points: %d points' %(len(zc))
	302	print ' x in [%f, %f]'%(min(Xc), max(Xc))
	303	print ' y in [%f, %f]'%(min(Yc), max(Yc))
	304	print ' z in [%f, %f]'%(min(zc), max(zc))
	305	print
	306
	307
	308	print 'Mesh: %d vertices' %(len(zeta))
	309	print ' xi in [%f, %f]'%(min(Xi), max(Xi))
	310	print ' eta in [%f, %f]'%(min(Eta), max(Eta))
	311	print ' zeta in [%f, %f]'%(min(zeta), max(zeta))
	312	print '+------------------------------------------------'
	313
	314	return vertex_attributes
	315
	316
	317
	318	def pts2rectangular(pts_name, M, N, alpha = DEFAULT_ALPHA,
	319	verbose = False, reduction = 1):
	320	"""Fits attributes from pts file to MxN rectangular mesh
	321
	322	Read pts file and create rectangular mesh of resolution MxN such that
	323	it covers all points specified in pts file.
	324
	325	FIXME: This may be a temporary function until we decide on
	326	netcdf formats etc
	327
	328	FIXME: Uses elevation hardwired
	329	"""
	330
	331	import mesh_factory
	332	from load_mesh.loadASCII import import_points_file
	333
	334	if verbose: print 'Read pts'
	335	points_dict = import_points_file(pts_name)
	336	#points, attributes = util.read_xya(pts_name)
	337
	338	#Reduce number of points a bit
	339	points = points_dict['pointlist'][::reduction]
	340	elevation = points_dict['attributelist']['elevation'] #Must be elevation
	341	elevation = elevation[::reduction]
	342
	343	if verbose: print 'Got %d data points' %len(points)
	344
	345	if verbose: print 'Create mesh'
	346	#Find extent
	347	max_x = min_x = points[0][0]
	348	max_y = min_y = points[0][1]
	349	for point in points[1:]:
	350	x = point[0]
	351	if x > max_x: max_x = x
	352	if x < min_x: min_x = x
	353	y = point[1]
	354	if y > max_y: max_y = y
	355	if y < min_y: min_y = y
	356
	357	#Create appropriate mesh
	358	vertex_coordinates, triangles, boundary =\
	359	mesh_factory.rectangular(M, N, max_x-min_x, max_y-min_y,
	360	(min_x, min_y))
	361
	362	#Fit attributes to mesh
	363	vertex_attributes = fit_to_mesh(vertex_coordinates,
	364	triangles,
	365	points,
	366	elevation, alpha=alpha, verbose=verbose)
	367
	368
	369
	370	return vertex_coordinates, triangles, boundary, vertex_attributes
	371
	372
	373	def _interpolation(args, *kwargs):
	374	"""Private function for use with caching. Reason is that classes
	375	may change their byte code between runs which is annoying.
	376	"""
	377
	378	return Interpolation(args, *kwargs)
	379
	380
	381	class Interpolation:
	382
	383	def __init__(self,
	384	vertex_coordinates,
	385	triangles,
[2781]	386	z,
[2695]	387	point_coordinates = None,
	388	alpha = None,
	389	verbose = False,
	390	expand_search = True,
	391	interp_only = False,
	392	max_points_per_cell = 30,
	393	mesh_origin = None,
	394	data_origin = None,
	395	precrop = False):
	396
	397
	398	""" Build interpolation matrix mapping from
	399	function values at vertices to function values at data points
	400
	401	Inputs:
	402
	403	vertex_coordinates: List of coordinate pairs [xi, eta] of
	404	points constituting mesh (or a an m x 2 Numeric array)
	405	Points may appear multiple times
	406	(e.g. if vertices have discontinuities)
	407
	408	triangles: List of 3-tuples (or a Numeric array) of
	409	integers representing indices of all vertices in the mesh.
	410
	411	point_coordinates: List of coordinate pairs [x, y] of
	412	data points (or an nx2 Numeric array)
	413	If point_coordinates is absent, only smoothing matrix will
	414	be built
	415
	416	alpha: Smoothing parameter
	417
	418	data_origin and mesh_origin are 3-tuples consisting of
	419	UTM zone, easting and northing. If specified
	420	point coordinates and vertex coordinates are assumed to be
	421	relative to their respective origins.
[2781]	422
	423	z: Single 1d vector or array of data at the point_coordinates.
[2695]	424
	425	"""
	426
	427	#Convert input to Numeric arrays
	428	triangles = ensure_numeric(triangles, Int)
	429	vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
	430
	431	#Build underlying mesh
	432	if verbose: print 'Building mesh'
	433	#self.mesh = General_mesh(vertex_coordinates, triangles,
	434	#FIXME: Trying the normal mesh while testing precrop,
	435	# The functionality of boundary_polygon is needed for that
	436
	437	#FIXME - geo ref does not have to go into mesh.
	438	# Change the point co-ords to conform to the
	439	# mesh co-ords early in the code
	440	if mesh_origin is None:
	441	geo = None
	442	else:
	443	geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
	444	self.mesh = Mesh(vertex_coordinates, triangles,
	445	geo_reference = geo)
	446
	447	if verbose: print 'Checking mesh integrity'
	448	self.mesh.check_integrity()
	449
	450	if verbose: print 'Mesh integrity checked'
	451
	452	self.data_origin = data_origin
	453
	454	self.point_indices = None
	455
	456	#Smoothing parameter
	457	if alpha is None:
	458	self.alpha = DEFAULT_ALPHA
	459	else:
	460	self.alpha = alpha
	461
	462
	463	if point_coordinates is not None:
	464	if verbose: print 'Building interpolation matrix'
	465	self.build_interpolation_matrix_A(point_coordinates,
[2781]	466	z,
[2695]	467	verbose = verbose,
	468	expand_search = expand_search,
	469	interp_only = interp_only,
	470	max_points_per_cell =\
	471	max_points_per_cell,
	472	data_origin = data_origin,
	473	precrop = precrop)
	474	#Build coefficient matrices
	475	if interp_only == False:
	476	self.build_coefficient_matrix_B(point_coordinates,
	477	verbose = verbose,
	478	expand_search = expand_search,
	479	max_points_per_cell =\
	480	max_points_per_cell,
	481	data_origin = data_origin,
	482	precrop = precrop)
	483	if verbose: print 'Finished interpolation'
	484
	485	def set_point_coordinates(self, point_coordinates,
	486	data_origin = None,
	487	verbose = False,
	488	precrop = True):
	489	"""
	490	A public interface to setting the point co-ordinates.
	491	"""
	492	if point_coordinates is not None:
	493	if verbose: print 'Building interpolation matrix'
	494	self.build_interpolation_matrix_A(point_coordinates,
[2781]	495	z,
[2695]	496	verbose = verbose,
	497	data_origin = data_origin,
	498	precrop = precrop)
	499	self.build_coefficient_matrix_B(point_coordinates, data_origin)
	500
	501	def build_coefficient_matrix_B(self, point_coordinates=None,
	502	verbose = False, expand_search = True,
	503	max_points_per_cell=30,
	504	data_origin = None,
	505	precrop = False):
	506	"""Build final coefficient matrix"""
	507
	508
	509	if self.alpha <> 0:
	510	if verbose: print 'Building smoothing matrix'
	511	self.build_smoothing_matrix_D()
	512
	513	if point_coordinates is not None:
	514	if self.alpha <> 0:
	515	self.B = self.AtA + self.alpha*self.D
	516	else:
	517	self.B = self.AtA
	518
	519	#Convert self.B matrix to CSR format for faster matrix vector
	520	self.B = Sparse_CSR(self.B)
	521
	522	def build_interpolation_matrix_A(self, point_coordinates,
[2781]	523	z,
[2695]	524	verbose = False, expand_search = True,
	525	max_points_per_cell=30,
	526	data_origin = None,
	527	precrop = False,
	528	interp_only = False):
	529	"""Build n x m interpolation matrix, where
	530	n is the number of data points and
	531	m is the number of basis functions phi_k (one per vertex)
	532
	533	This algorithm uses a quad tree data structure for fast binning of data points
	534	origin is a 3-tuple consisting of UTM zone, easting and northing.
	535	If specified coordinates are assumed to be relative to this origin.
	536
	537	This one will override any data_origin that may be specified in
	538	interpolation instance
	539
	540	"""
	541
	542	from pyvolution.quad import build_quadtree
	543	from utilities.polygon import inside_polygon
	544
[2781]	545	#Convert input to Numeric arrays
	546	z = ensure_numeric(z, Float)
	547
[2695]	548	#Convert input to Numeric arrays just in case.
	549	point_coordinates = ensure_numeric(point_coordinates, Float)
	550
	551	#Build n x m interpolation matrix
	552	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
	553	n = point_coordinates.shape[0] #Nbr of data points
[2781]	554	if len(z.shape) > 1:
	555	att_num = z.shape[1]
	556	else:
	557	att_num = 0
	558
	559	assert z.shape[0] == point_coordinates.shape[0]
	560	if verbose: print 'Number of attributes: %d' %att_num
[2695]	561	if verbose: print 'Number of datapoints: %d' %n
	562	if verbose: print 'Number of basis functions: %d' %m
	563
	564	self.A = Sparse(n,m)
	565	self.AtA = Sparse(m,m)
[2781]	566	self.Atz = zeros((m,att_num), Float)
[2897]	567	print "self.Atz",self.Atz
	568	self.Atz[0] = 100
	569	print "self.Atz",self.Atz
[2695]	570
[2750]	571	#Build quad tree of vertices
[2695]	572	root = build_quadtree(self.mesh,
	573	max_points_per_cell = max_points_per_cell)
	574	#root.show()
	575	self.expanded_quad_searches = []
	576	#Compute matrix elements
	577	for i in range(n):
	578	#For each data_coordinate point
	579
	580	if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
	581	x = point_coordinates[i]
	582
	583	#Find vertices near x
	584	candidate_vertices = root.search(x[0], x[1])
	585	is_more_elements = True
	586
	587	element_found, sigma0, sigma1, sigma2, k = \
	588	self.search_triangles_of_vertices(candidate_vertices, x)
	589	first_expansion = True
	590	while not element_found and is_more_elements and expand_search:
	591	#if verbose: print 'Expanding search'
	592	if first_expansion == True:
	593	self.expanded_quad_searches.append(1)
	594	first_expansion = False
	595	else:
	596	end = len(self.expanded_quad_searches) - 1
	597	assert end >= 0
	598	self.expanded_quad_searches[end] += 1
	599	candidate_vertices, branch = root.expand_search()
	600	if branch == []:
	601	# Searching all the verts from the root cell that haven't
	602	# been searched. This is the last try
	603	element_found, sigma0, sigma1, sigma2, k = \
	604	self.search_triangles_of_vertices(candidate_vertices, x)
	605	is_more_elements = False
	606	else:
	607	element_found, sigma0, sigma1, sigma2, k = \
	608	self.search_triangles_of_vertices(candidate_vertices, x)
	609
	610
	611	#Update interpolation matrix A if necessary
	612	if element_found is True:
	613	#Assign values to matrix A
	614
	615	j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
	616	j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
	617	j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2
	618
	619	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
	620	js = [j0,j1,j2]
	621
	622	for j in js:
	623	self.A[i,j] = sigmas[j]
[2781]	624	#self.Atz[n-i-1] += sigmas[j]*z[n-j-1]
	625	#print "i",i
	626	#print "j",j
	627	#print "self.Atz",self.Atz
	628	#print "z",z
	629	#a=self.Atz[j]
	630	#b=z[i]
	631	self.Atz[j] += sigmas[j]*z[i]
	632	#print "m-i-1",m-i-1
	633	#print "n-j-1",n-j-1
	634	#print "self.A[i,j]",self.A[i,j]
	635	#print "z[n-j-1]",z[n-j-1]
	636	#print "z[i]",z[i]
	637	#print "",
	638	#print "self.Atz",self.Atz
[2695]	639	for k in js:
	640	if interp_only == False:
	641	self.AtA[j,k] += sigmas[j]*sigmas[k]
	642	else:
	643	pass
	644	#Ok if there is no triangle for datapoint
	645	#(as in brute force version)
	646	#raise 'Could not find triangle for point', x
	647
	648
	649
	650	def search_triangles_of_vertices(self, candidate_vertices, x):
	651	#Find triangle containing x:
	652	element_found = False
	653
	654	# This will be returned if element_found = False
	655	sigma2 = -10.0
	656	sigma0 = -10.0
	657	sigma1 = -10.0
	658	k = -10.0
	659	#print "$ candidate_vertices", candidate_vertices
	660	#For all vertices in same cell as point x
	661	for v in candidate_vertices:
	662	#FIXME (DSG-DSG): this catches verts with no triangle.
	663	#Currently pmesh is producing these.
	664	#this should be stopped,
	665	if self.mesh.vertexlist[v] is None:
	666	continue
	667	#for each triangle id (k) which has v as a vertex
	668	for k, _ in self.mesh.vertexlist[v]:
	669
	670	#Get the three vertex_points of candidate triangle
	671	xi0 = self.mesh.get_vertex_coordinate(k, 0)
	672	xi1 = self.mesh.get_vertex_coordinate(k, 1)
	673	xi2 = self.mesh.get_vertex_coordinate(k, 2)
	674
	675	#print "PDSG - k", k
	676	#print "PDSG - xi0", xi0
	677	#print "PDSG - xi1", xi1
	678	#print "PDSG - xi2", xi2
	679	#print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
	680	# % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
	681
	682	#Get the three normals
	683	n0 = self.mesh.get_normal(k, 0)
	684	n1 = self.mesh.get_normal(k, 1)
	685	n2 = self.mesh.get_normal(k, 2)
	686
	687
	688	#Compute interpolation
	689	sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
	690	sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
	691	sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
	692
	693	#print "PDSG - sigma0", sigma0
	694	#print "PDSG - sigma1", sigma1
	695	#print "PDSG - sigma2", sigma2
	696
	697	#FIXME: Maybe move out to test or something
	698	epsilon = 1.0e-6
	699	assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
	700
	701	#Check that this triangle contains the data point
	702
	703	#Sigmas can get negative within
	704	#machine precision on some machines (e.g nautilus)
	705	#Hence the small eps
	706	eps = 1.0e-15
	707	if sigma0 >= -eps and sigma1 >= -eps and sigma2 >= -eps:
	708	element_found = True
	709	break
	710
	711	if element_found is True:
	712	#Don't look for any other triangle
	713	break
	714	return element_found, sigma0, sigma1, sigma2, k
	715
	716
	717
	718	def build_interpolation_matrix_A_brute(self, point_coordinates):
	719	"""Build n x m interpolation matrix, where
	720	n is the number of data points and
	721	m is the number of basis functions phi_k (one per vertex)
	722
	723	This is the brute force which is too slow for large problems,
[2781]	724	but could be used for testing
[2695]	725	"""
	726
	727
	728	#Convert input to Numeric arrays
	729	point_coordinates = ensure_numeric(point_coordinates, Float)
	730
	731	#Build n x m interpolation matrix
	732	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
	733	n = point_coordinates.shape[0] #Nbr of data points
	734
	735	self.A = Sparse(n,m)
	736	self.AtA = Sparse(m,m)
	737
	738	#Compute matrix elements
	739	for i in range(n):
	740	#For each data_coordinate point
	741
	742	x = point_coordinates[i]
	743	element_found = False
	744	k = 0
	745	while not element_found and k < len(self.mesh):
	746	#For each triangle (brute force)
	747	#FIXME: Real algorithm should only visit relevant triangles
	748
	749	#Get the three vertex_points
	750	xi0 = self.mesh.get_vertex_coordinate(k, 0)
	751	xi1 = self.mesh.get_vertex_coordinate(k, 1)
	752	xi2 = self.mesh.get_vertex_coordinate(k, 2)
	753
	754	#Get the three normals
	755	n0 = self.mesh.get_normal(k, 0)
	756	n1 = self.mesh.get_normal(k, 1)
	757	n2 = self.mesh.get_normal(k, 2)
	758
	759	#Compute interpolation
	760	sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
	761	sigma0 = dot((x-xi1), n0)/dot((xi0-xi1), n0)
	762	sigma1 = dot((x-xi2), n1)/dot((xi1-xi2), n1)
	763
	764	#FIXME: Maybe move out to test or something
	765	epsilon = 1.0e-6
	766	assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
	767
	768	#Check that this triangle contains data point
	769	if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
	770	element_found = True
	771	#Assign values to matrix A
	772
	773	j0 = self.mesh.triangles[k,0] #Global vertex id
	774	#self.A[i, j0] = sigma0
	775
	776	j1 = self.mesh.triangles[k,1] #Global vertex id
	777	#self.A[i, j1] = sigma1
	778
	779	j2 = self.mesh.triangles[k,2] #Global vertex id
	780	#self.A[i, j2] = sigma2
	781
	782	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
	783	js = [j0,j1,j2]
	784
	785	for j in js:
	786	self.A[i,j] = sigmas[j]
	787	for k in js:
	788	self.AtA[j,k] += sigmas[j]*sigmas[k]
	789	k = k+1
	790
	791
	792
	793	def get_A(self):
	794	return self.A.todense()
	795
	796	def get_B(self):
	797	return self.B.todense()
	798
	799	def get_D(self):
	800	return self.D.todense()
	801
	802	#FIXME: Remember to re-introduce the 1/n factor in the
	803	#interpolation term
	804
	805	def build_smoothing_matrix_D(self):
	806	"""Build m x m smoothing matrix, where
	807	m is the number of basis functions phi_k (one per vertex)
	808
	809	The smoothing matrix is defined as
	810
	811	D = D1 + D2
	812
	813	where
	814
	815	[D1]_{k,l} = \int_\Omega
	816	\frac{\partial \phi_k}{\partial x}
	817	\frac{\partial \phi_l}{\partial x}\,
	818	dx dy
	819
	820	[D2]_{k,l} = \int_\Omega
	821	\frac{\partial \phi_k}{\partial y}
	822	\frac{\partial \phi_l}{\partial y}\,
	823	dx dy
	824
	825
	826	The derivatives \frac{\partial \phi_k}{\partial x},
	827	\frac{\partial \phi_k}{\partial x} for a particular triangle
	828	are obtained by computing the gradient a_k, b_k for basis function k
	829	"""
	830
	831	#FIXME: algorithm might be optimised by computing local 9x9
	832	#"element stiffness matrices:
	833
	834	m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)
	835
	836	self.D = Sparse(m,m)
	837
	838	#For each triangle compute contributions to D = D1+D2
	839	for i in range(len(self.mesh)):
	840
	841	#Get area
	842	area = self.mesh.areas[i]
	843
	844	#Get global vertex indices
	845	v0 = self.mesh.triangles[i,0]
	846	v1 = self.mesh.triangles[i,1]
	847	v2 = self.mesh.triangles[i,2]
	848
	849	#Get the three vertex_points
	850	xi0 = self.mesh.get_vertex_coordinate(i, 0)
	851	xi1 = self.mesh.get_vertex_coordinate(i, 1)
	852	xi2 = self.mesh.get_vertex_coordinate(i, 2)
	853
	854	#Compute gradients for each vertex
	855	a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
	856	1, 0, 0)
	857
	858	a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
	859	0, 1, 0)
	860
	861	a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
	862	0, 0, 1)
	863
	864	#Compute diagonal contributions
	865	self.D[v0,v0] += (a0a0 + b0b0)*area
	866	self.D[v1,v1] += (a1a1 + b1b1)*area
	867	self.D[v2,v2] += (a2a2 + b2b2)*area
	868
	869	#Compute contributions for basis functions sharing edges
	870	e01 = (a0a1 + b0b1)*area
	871	self.D[v0,v1] += e01
	872	self.D[v1,v0] += e01
	873
	874	e12 = (a1a2 + b1b2)*area
	875	self.D[v1,v2] += e12
	876	self.D[v2,v1] += e12
	877
	878	e20 = (a2a0 + b2b0)*area
	879	self.D[v2,v0] += e20
	880	self.D[v0,v2] += e20
	881
	882
	883	def fit(self, z):
	884	"""Fit a smooth surface to given 1d array of data points z.
	885
	886	The smooth surface is computed at each vertex in the underlying
	887	mesh using the formula given in the module doc string.
	888
	889	Pre Condition:
	890	self.A, self.AtA and self.B have been initialised
	891
	892	Inputs:
	893	z: Single 1d vector or array of data at the point_coordinates.
	894	"""
	895
	896	#Convert input to Numeric arrays
	897	z = ensure_numeric(z, Float)
	898
	899	if len(z.shape) > 1 :
	900	raise VectorShapeError, 'Can only deal with 1d data vector'
	901
	902	if self.point_indices is not None:
	903	#Remove values for any points that were outside mesh
	904	z = take(z, self.point_indices)
	905
	906	#Compute right hand side based on data
	907	#FIXME (DSG-DsG): could Sparse_CSR be used here? Use this format
	908	# after a matrix is built, before calcs.
[2781]	909	Atz = self.Atz
	910	#print "self.get_A()", self.get_A()
	911	#print "self.Atz",self.Atz
	912	#print 'z', z
[2695]	913	Atz = self.A.trans_mult(z)
[2781]	914	#print "self.A.trans_mult(z)",self.A.trans_mult(z)
[2695]	915
	916
	917	#Check sanity
	918	n, m = self.A.shape
	919	if n<m and self.alpha == 0.0:
	920	msg = 'ERROR (least_squares): Too few data points\n'
	921	msg += 'There are only %d data points and alpha == 0. ' %n
	922	msg += 'Need at least %d\n' %m
	923	msg += 'Alternatively, set smoothing parameter alpha to a small '
	924	msg += 'positive value,\ne.g. 1.0e-3.'
	925	raise Exception(msg)
	926
	927
	928
	929	return conjugate_gradient(self.B, Atz, Atz, imax=2*len(Atz) )
	930	#FIXME: Should we store the result here for later use? (ON)
	931
	932
	933	def fit_points(self, z, verbose=False):
	934	"""Like fit, but more robust when each point has two or more attributes
	935	FIXME (Ole): The name fit_points doesn't carry any meaning
	936	for me. How about something like fit_multiple or fit_columns?
	937	"""
	938
	939	try:
	940	if verbose: print 'Solving penalised least_squares problem'
	941	return self.fit(z)
	942	except VectorShapeError, e:
	943	# broadcasting is not supported.
	944
	945	#Convert input to Numeric arrays
	946	z = ensure_numeric(z, Float)
	947
	948	#Build n x m interpolation matrix
	949	m = self.mesh.coordinates.shape[0] #Number of vertices
	950	n = z.shape[1] #Number of data points
	951
	952	f = zeros((m,n), Float) #Resulting columns
	953
	954	for i in range(z.shape[1]):
	955	f[:,i] = self.fit(z[:,i])
	956
	957	return f
	958
	959
	960	def interpolate(self, f):
	961	"""Evaluate smooth surface f at data points implied in self.A.
	962
	963	The mesh values representing a smooth surface are
	964	assumed to be specified in f. This argument could,
	965	for example have been obtained from the method self.fit()
	966
	967	Pre Condition:
	968	self.A has been initialised
	969
	970	Inputs:
	971	f: Vector or array of data at the mesh vertices.
	972	If f is an array, interpolation will be done for each column as
	973	per underlying matrix-matrix multiplication
	974
	975	Output:
	976	Interpolated values at data points implied in self.A
	977
	978	"""
	979
	980	return self.A * f
	981
	982	def cull_outsiders(self, f):
	983	pass
	984
	985
	986
	987
	988	class Interpolation_function:
	989	"""Interpolation_function - creates callable object f(t, id) or f(t,x,y)
	990	which is interpolated from time series defined at vertices of
	991	triangular mesh (such as those stored in sww files)
	992
	993	Let m be the number of vertices, n the number of triangles
	994	and p the number of timesteps.
	995
	996	Mandatory input
	997	time: px1 array of monotonously increasing times (Float)
	998	quantities: Dictionary of arrays or 1 array (Float)
	999	The arrays must either have dimensions pxm or mx1.
	1000	The resulting function will be time dependent in
	1001	the former case while it will be constant with
	1002	respect to time in the latter case.
	1003
	1004	Optional input:
	1005	quantity_names: List of keys into the quantities dictionary
	1006	vertex_coordinates: mx2 array of coordinates (Float)
	1007	triangles: nx3 array of indices into vertex_coordinates (Int)
	1008	interpolation_points: Nx2 array of coordinates to be interpolated to
	1009	verbose: Level of reporting
	1010
	1011
	1012	The quantities returned by the callable object are specified by
	1013	the list quantities which must contain the names of the
	1014	quantities to be returned and also reflect the order, e.g. for
	1015	the shallow water wave equation, on would have
	1016	quantities = ['stage', 'xmomentum', 'ymomentum']
	1017
	1018	The parameter interpolation_points decides at which points interpolated
	1019	quantities are to be computed whenever object is called.
	1020	If None, return average value
	1021	"""
	1022
	1023
	1024
	1025	def __init__(self,
	1026	time,
	1027	quantities,
	1028	quantity_names = None,
	1029	vertex_coordinates = None,
	1030	triangles = None,
	1031	interpolation_points = None,
	1032	verbose = False):
	1033	"""Initialise object and build spatial interpolation if required
	1034	"""
	1035
	1036	from Numeric import array, zeros, Float, alltrue, concatenate,\
	1037	reshape, ArrayType
	1038
	1039
	1040	from config import time_format
	1041	import types
	1042
	1043
	1044
	1045	#Check temporal info
	1046	time = ensure_numeric(time)
	1047	msg = 'Time must be a monotonuosly '
	1048	msg += 'increasing sequence %s' %time
	1049	assert alltrue(time[1:] - time[:-1] >= 0 ), msg
	1050
	1051
	1052	#Check if quantities is a single array only
	1053	if type(quantities) != types.DictType:
	1054	quantities = ensure_numeric(quantities)
	1055	quantity_names = ['Attribute']
	1056
	1057	#Make it a dictionary
	1058	quantities = {quantity_names[0]: quantities}
	1059
	1060
	1061	#Use keys if no names are specified
	1062	if quantity_names is None:
	1063	quantity_names = quantities.keys()
	1064
	1065
	1066	#Check spatial info
	1067	if vertex_coordinates is None:
	1068	self.spatial = False
	1069	else:
	1070	vertex_coordinates = ensure_numeric(vertex_coordinates)
	1071
	1072	assert triangles is not None, 'Triangles array must be specified'
	1073	triangles = ensure_numeric(triangles)
	1074	self.spatial = True
	1075
	1076
	1077
	1078	#Save for use with statistics
	1079	self.quantity_names = quantity_names
	1080	self.quantities = quantities
	1081	self.vertex_coordinates = vertex_coordinates
	1082	self.interpolation_points = interpolation_points
[2884]	1083	self.time = time[:] # Time assumed to be relative to starttime
[2695]	1084	self.index = 0 # Initial time index
	1085	self.precomputed_values = {}
	1086
	1087
	1088
	1089	#Precomputed spatial interpolation if requested
	1090	if interpolation_points is not None:
	1091	if self.spatial is False:
	1092	raise 'Triangles and vertex_coordinates must be specified'
	1093
	1094	try:
	1095	self.interpolation_points = ensure_numeric(interpolation_points)
	1096	except:
	1097	msg = 'Interpolation points must be an N x 2 Numeric array '+\
	1098	'or a list of points\n'
	1099	msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\
	1100	'...')
	1101	raise msg
	1102
	1103
	1104	m = len(self.interpolation_points)
[2884]	1105	p = len(self.time)
[2695]	1106
	1107	for name in quantity_names:
	1108	self.precomputed_values[name] = zeros((p, m), Float)
	1109
	1110	#Build interpolator
	1111	interpol = Interpolation(vertex_coordinates,
	1112	triangles,
	1113	point_coordinates = \
	1114	self.interpolation_points,
	1115	alpha = 0,
	1116	precrop = False,
	1117	verbose = verbose)
	1118
	1119	if verbose: print 'Interpolate'
[2884]	1120	for i, t in enumerate(self.time):
[2695]	1121	#Interpolate quantities at this timestep
	1122	if verbose and i%((p+10)/10)==0:
	1123	print ' time step %d of %d' %(i, p)
	1124
	1125	for name in quantity_names:
	1126	if len(quantities[name].shape) == 2:
	1127	result = interpol.interpolate(quantities[name][i,:])
	1128	else:
	1129	#Assume no time dependency
	1130	result = interpol.interpolate(quantities[name][:])
	1131
	1132	self.precomputed_values[name][i, :] = result
	1133
	1134
	1135
	1136	#Report
	1137	if verbose:
	1138	print self.statistics()
	1139	#self.print_statistics()
	1140
	1141	else:
	1142	#Store quantitites as is
	1143	for name in quantity_names:
	1144	self.precomputed_values[name] = quantities[name]
	1145
	1146
	1147	#else:
	1148	# #Return an average, making this a time series
	1149	# for name in quantity_names:
[2884]	1150	# self.values[name] = zeros(len(self.time), Float)
[2695]	1151	#
	1152	# if verbose: print 'Compute mean values'
[2884]	1153	# for i, t in enumerate(self.time):
	1154	# if verbose: print ' time step %d of %d' %(i, len(self.time))
[2695]	1155	# for name in quantity_names:
	1156	# self.values[name][i] = mean(quantities[name][i,:])
	1157
	1158
	1159
	1160
	1161	def __repr__(self):
	1162	#return 'Interpolation function (spatio-temporal)'
	1163	return self.statistics()
	1164
	1165
	1166	def __call__(self, t, point_id = None, x = None, y = None):
	1167	"""Evaluate f(t), f(t, point_id) or f(t, x, y)
	1168
	1169	Inputs:
	1170	t: time - Model time. Must lie within existing timesteps
	1171	point_id: index of one of the preprocessed points.
	1172	x, y: Overrides location, point_id ignored
	1173
	1174	If spatial info is present and all of x,y,point_id
	1175	are None an exception is raised
	1176
	1177	If no spatial info is present, point_id and x,y arguments are ignored
	1178	making f a function of time only.
	1179
	1180
	1181	FIXME: point_id could also be a slice
	1182	FIXME: What if x and y are vectors?
	1183	FIXME: What about f(x,y) without t?
	1184	"""
	1185
	1186	from math import pi, cos, sin, sqrt
	1187	from Numeric import zeros, Float
	1188	from utilities.numerical_tools import mean
	1189
	1190	if self.spatial is True:
	1191	if point_id is None:
	1192	if x is None or y is None:
	1193	msg = 'Either point_id or x and y must be specified'
	1194	raise Exception(msg)
	1195	else:
	1196	if self.interpolation_points is None:
	1197	msg = 'Interpolation_function must be instantiated ' +\
	1198	'with a list of interpolation points before parameter ' +\
	1199	'point_id can be used'
	1200	raise Exception(msg)
	1201
	1202
[2884]	1203	msg = 'Time interval [%s:%s]' %(self.time[0], self.time[-1])
[2695]	1204	msg += ' does not match model time: %s\n' %t
[2884]	1205	if t < self.time[0]: raise Exception(msg)
	1206	if t > self.time[-1]: raise Exception(msg)
[2695]	1207
	1208	oldindex = self.index #Time index
	1209
	1210	#Find current time slot
[2884]	1211	while t > self.time[self.index]: self.index += 1
	1212	while t < self.time[self.index]: self.index -= 1
[2695]	1213
[2884]	1214	if t == self.time[self.index]:
[2695]	1215	#Protect against case where t == T[-1] (last time)
	1216	# - also works in general when t == T[i]
	1217	ratio = 0
	1218	else:
	1219	#t is now between index and index+1
[2884]	1220	ratio = (t - self.time[self.index])/\
	1221	(self.time[self.index+1] - self.time[self.index])
[2695]	1222
	1223	#Compute interpolated values
	1224	q = zeros(len(self.quantity_names), Float)
	1225
	1226	for i, name in enumerate(self.quantity_names):
	1227	Q = self.precomputed_values[name]
	1228
	1229	if self.spatial is False:
	1230	#If there is no spatial info
	1231	assert len(Q.shape) == 1
	1232
	1233	Q0 = Q[self.index]
	1234	if ratio > 0: Q1 = Q[self.index+1]
	1235
	1236	else:
	1237	if x is not None and y is not None:
	1238	#Interpolate to x, y
	1239
	1240	raise 'x,y interpolation not yet implemented'
	1241	else:
	1242	#Use precomputed point
	1243	Q0 = Q[self.index, point_id]
	1244	if ratio > 0: Q1 = Q[self.index+1, point_id]
	1245
	1246	#Linear temporal interpolation
	1247	if ratio > 0:
	1248	q[i] = Q0 + ratio*(Q1 - Q0)
	1249	else:
	1250	q[i] = Q0
	1251
	1252
	1253	#Return vector of interpolated values
	1254	#if len(q) == 1:
	1255	# return q[0]
	1256	#else:
	1257	# return q
	1258
	1259
	1260	#Return vector of interpolated values
	1261	#FIXME:
	1262	if self.spatial is True:
	1263	return q
	1264	else:
	1265	#Replicate q according to x and y
	1266	#This is e.g used for Wind_stress
	1267	if x is None or y is None:
	1268	return q
	1269	else:
	1270	try:
	1271	N = len(x)
	1272	except:
	1273	return q
	1274	else:
	1275	from Numeric import ones, Float
	1276	#x is a vector - Create one constant column for each value
	1277	N = len(x)
	1278	assert len(y) == N, 'x and y must have same length'
	1279	res = []
	1280	for col in q:
	1281	res.append(col*ones(N, Float))
	1282
	1283	return res
	1284
	1285
	1286	def statistics(self):
	1287	"""Output statistics about interpolation_function
	1288	"""
	1289
	1290	vertex_coordinates = self.vertex_coordinates
	1291	interpolation_points = self.interpolation_points
	1292	quantity_names = self.quantity_names
	1293	quantities = self.quantities
	1294	precomputed_values = self.precomputed_values
	1295
	1296	x = vertex_coordinates[:,0]
	1297	y = vertex_coordinates[:,1]
	1298
	1299	str = '------------------------------------------------\n'
	1300	str += 'Interpolation_function (spatio-temporal) statistics:\n'
	1301	str += ' Extent:\n'
	1302	str += ' x in [%f, %f], len(x) == %d\n'\
	1303	%(min(x), max(x), len(x))
	1304	str += ' y in [%f, %f], len(y) == %d\n'\
	1305	%(min(y), max(y), len(y))
	1306	str += ' t in [%f, %f], len(t) == %d\n'\
[2884]	1307	%(min(self.time), max(self.time), len(self.time))
[2695]	1308	str += ' Quantities:\n'
	1309	for name in quantity_names:
	1310	q = quantities[name][:].flat
	1311	str += ' %s in [%f, %f]\n' %(name, min(q), max(q))
	1312
	1313	if interpolation_points is not None:
	1314	str += ' Interpolation points (xi, eta):'\
	1315	' number of points == %d\n' %interpolation_points.shape[0]
	1316	str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]),
	1317	max(interpolation_points[:,0]))
	1318	str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]),
	1319	max(interpolation_points[:,1]))
	1320	str += ' Interpolated quantities (over all timesteps):\n'
	1321
	1322	for name in quantity_names:
	1323	q = precomputed_values[name][:].flat
	1324	str += ' %s at interpolation points in [%f, %f]\n'\
	1325	%(name, min(q), max(q))
	1326	str += '------------------------------------------------\n'
	1327
	1328	return str
	1329
	1330	#FIXME: Delete
	1331	#print '------------------------------------------------'
	1332	#print 'Interpolation_function statistics:'
	1333	#print ' Extent:'
	1334	#print ' x in [%f, %f], len(x) == %d'\
	1335	# %(min(x), max(x), len(x))
	1336	#print ' y in [%f, %f], len(y) == %d'\
	1337	# %(min(y), max(y), len(y))
	1338	#print ' t in [%f, %f], len(t) == %d'\
[2884]	1339	# %(min(self.time), max(self.time), len(self.time))
[2695]	1340	#print ' Quantities:'
	1341	#for name in quantity_names:
	1342	# q = quantities[name][:].flat
	1343	# print ' %s in [%f, %f]' %(name, min(q), max(q))
	1344	#print ' Interpolation points (xi, eta):'\
	1345	# ' number of points == %d ' %interpolation_points.shape[0]
	1346	#print ' xi in [%f, %f]' %(min(interpolation_points[:,0]),
	1347	# max(interpolation_points[:,0]))
	1348	#print ' eta in [%f, %f]' %(min(interpolation_points[:,1]),
	1349	# max(interpolation_points[:,1]))
	1350	#print ' Interpolated quantities (over all timesteps):'
	1351	#
	1352	#for name in quantity_names:
	1353	# q = precomputed_values[name][:].flat
	1354	# print ' %s at interpolation points in [%f, %f]'\
	1355	# %(name, min(q), max(q))
	1356	#print '------------------------------------------------'
	1357
	1358
	1359	#-------------------------------------------------------------
	1360	if __name__ == "__main__":
	1361	"""
	1362	Load in a mesh and data points with attributes.
	1363	Fit the attributes to the mesh.
	1364	Save a new mesh file.
	1365	"""
	1366	import os, sys
	1367	usage = "usage: %s mesh_input.tsh point.xya mesh_output.tsh [expand\|no_expand][vervose\|non_verbose] [alpha] [display_errors\|no_display_errors]"\
	1368	%os.path.basename(sys.argv[0])
	1369
	1370	if len(sys.argv) < 4:
	1371	print usage
	1372	else:
	1373	mesh_file = sys.argv[1]
	1374	point_file = sys.argv[2]
	1375	mesh_output_file = sys.argv[3]
	1376
	1377	expand_search = False
	1378	if len(sys.argv) > 4:
	1379	if sys.argv[4][0] == "e" or sys.argv[4][0] == "E":
	1380	expand_search = True
	1381	else:
	1382	expand_search = False
	1383
	1384	verbose = False
	1385	if len(sys.argv) > 5:
	1386	if sys.argv[5][0] == "n" or sys.argv[5][0] == "N":
	1387	verbose = False
	1388	else:
	1389	verbose = True
	1390
	1391	if len(sys.argv) > 6:
	1392	alpha = sys.argv[6]
	1393	else:
	1394	alpha = DEFAULT_ALPHA
	1395
	1396	# This is used more for testing
	1397	if len(sys.argv) > 7:
	1398	if sys.argv[7][0] == "n" or sys.argv[5][0] == "N":
	1399	display_errors = False
	1400	else:
	1401	display_errors = True
	1402
	1403	t0 = time.time()
	1404	try:
	1405	fit_to_mesh_file(mesh_file,
	1406	point_file,
	1407	mesh_output_file,
	1408	alpha,
	1409	verbose= verbose,
	1410	expand_search = expand_search,
	1411	display_errors = display_errors)
	1412	except IOError,e:
	1413	import sys; sys.exit(1)
	1414
	1415	print 'That took %.2f seconds' %(time.time()-t0)
	1416

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

Download in other formats: