Context Navigation

source: inundation/pyvolution/test_least_squares.py @ 2945

Last change on this file since 2945 was 2778, checked in by steve, 19 years ago
On our 64 bit machine, ran into problem in pmesh/mesh.py where seg[0], seg[1], triangle, and neighbor where not seen as integers. Had to wrap them in int(..)
File size: 58.9 KB

Rev	Line
[2526]	1	#!/usr/bin/env python
	2
	3	#TEST
	4	import sys
	5	import unittest
	6	from math import sqrt
	7
	8
[2753]	9	from pyvolution.least_squares import *
[2526]	10	from Numeric import allclose, array, transpose
	11
	12	from coordinate_transforms.geo_reference import Geo_reference
	13
	14	def distance(x, y):
	15	return sqrt( sum( (array(x)-array(y))**2 ))
	16
	17	def linear_function(point):
	18	point = array(point)
	19	return point[:,0]+point[:,1]
	20
	21
	22	class Test_Least_Squares(unittest.TestCase):
	23
	24	def setUp(self):
	25	pass
	26
	27	def tearDown(self):
	28	pass
	29
	30	def test_datapoint_at_centroid(self):
	31	a = [0.0, 0.0]
	32	b = [0.0, 2.0]
	33	c = [2.0,0.0]
	34	points = [a, b, c]
	35	vertices = [ [1,0,2] ] #bac
	36
	37	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	38
[2778]	39
[2526]	40	interp = Interpolation(points, vertices, data)
	41	assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
	42
	43
	44	def test_quad_tree(self):
	45	p0 = [-10.0, -10.0]
	46	p1 = [20.0, -10.0]
	47	p2 = [-10.0, 20.0]
	48	p3 = [10.0, 50.0]
	49	p4 = [30.0, 30.0]
	50	p5 = [50.0, 10.0]
	51	p6 = [40.0, 60.0]
	52	p7 = [60.0, 40.0]
	53	p8 = [-66.0, 20.0]
	54	p9 = [10.0, -66.0]
	55
	56	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	57	triangles = [ [0, 1, 2],
	58	[3, 2, 4],
	59	[4, 2, 1],
	60	[4, 1, 5],
	61	[3, 4, 6],
	62	[6, 4, 7],
	63	[7, 4, 5],
	64	[8, 0, 2],
	65	[0, 9, 1]]
	66
	67	data = [ [4,4] ]
	68	interp = Interpolation(points, triangles, data, alpha = 0.0,
	69	max_points_per_cell = 4)
	70	#print "PDSG - interp.get_A()", interp.get_A()
	71	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	72	0., 0. , 0., 0., 0., 0.]]
	73	assert allclose(interp.get_A(), answer)
	74	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	75	#print "PDSG - interp.get_A()", interp.get_A()
	76	answer = [ [ 0.0, 0.0, 0.0, 0.,
	77	0., 0. , 0., 0., 0., 0.]]
	78	assert allclose(interp.get_A(), answer)
	79
	80
	81	#point outside of quad tree root cell
	82	interp.set_point_coordinates([[-70, -70]])
	83	#print "PDSG - interp.get_A()", interp.get_A()
	84	answer = [ [ 0.0, 0.0, 0.0, 0.,
	85	0., 0. , 0., 0., 0., 0.]]
	86	assert allclose(interp.get_A(), answer)
	87
	88	def test_expand_search(self):
	89	p0 = [-10.0, -10.0]
	90	p1 = [20.0, -10.0]
	91	p2 = [-10.0, 20.0]
	92	p3 = [10.0, 50.0]
	93	p4 = [30.0, 30.0]
	94	p5 = [50.0, 10.0]
	95	p6 = [40.0, 60.0]
	96	p7 = [60.0, 40.0]
	97	p8 = [-66.0, 20.0]
	98	p9 = [10.0, -66.0]
	99
	100	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	101	triangles = [ [0, 1, 2],
	102	[3, 2, 4],
	103	[4, 2, 1],
	104	[4, 1, 5],
	105	[3, 4, 6],
	106	[6, 4, 7],
	107	[7, 4, 5],
	108	[8, 0, 2],
	109	[0, 9, 1]]
	110
	111	data = [ [4,4],
	112	[-30,10],
	113	[-20,0],
	114	[-20,10],
	115	[0,30],
	116	[10,-40],
	117	[10,-30],
	118	[10,-20],
	119	[10,10],
	120	[10,20],
	121	[10,30],
	122	[10,40],
	123	[20,10],
	124	[25,45],
	125	[30,0],
	126	[30,10],
	127	[30,30],
	128	[30,40],
	129	[30,50],
	130	[40,10],
	131	[40,30],
	132	[40,40],
	133	[40,50],
	134	[50,20],
	135	[50,30],
	136	[50,40],
	137	[50,50],
	138	[30,0],
	139	[-20,-20]]
	140	point_attributes = [ -400000,
	141	10,
	142	10,
	143	10,
	144	10,
	145	10,
	146	10,
	147	10,
	148	10,
	149	10,
	150	10,
	151	10,
	152	10,
	153	10,
	154	10,
	155	10,
	156	10,
	157	10,
	158	10,
	159	10,
	160	10,
	161	10,
	162	10,
	163	10,
	164	10,
	165	10,
	166	10,
	167	10,
	168	99]
	169
	170	interp = Interpolation(points, triangles, data,
	171	alpha=0.0, expand_search=False, #verbose = True, #False,
	172	max_points_per_cell = 4)
	173	calc = interp.fit_points(point_attributes, )
	174	#print "calc",calc
	175
	176	# the point at 4,4 is ignored. An expanded search has to be done
	177	# to fine which triangel it's in.
	178	# An expanded search isn't done to find that the last point
	179	# isn't in the mesh. But this isn't tested.
	180	answer= [ 10,
	181	10,
	182	10,
	183	10,
	184	10,
	185	10,
	186	10,
	187	10,
	188	10,
	189	10]
	190	assert allclose(calc, answer)
	191
	192	def test_quad_treeII(self):
	193	p0 = [-66.0, 14.0]
	194	p1 = [14.0, -66.0]
	195	p2 = [14.0, 14.0]
	196	p3 = [60.0, 20.0]
	197	p4 = [10.0, 60.0]
	198	p5 = [60.0, 60.0]
	199
	200	points = [p0, p1, p2, p3, p4, p5]
	201	triangles = [ [0, 1, 2],
	202	[3, 2, 1],
	203	[0, 2, 4],
	204	[4, 2, 5],
	205	[5, 2, 3]]
	206
	207	data = [ [-26.0,-26.0] ]
	208	interp = Interpolation(points, triangles, data, alpha = 0.0,
	209	max_points_per_cell = 4)
	210	#print "PDSG - interp.get_A()", interp.get_A()
	211	answer = [ [ 0.5, 0.5, 0.0, 0.,
	212	0., 0.]]
	213	assert allclose(interp.get_A(), answer)
	214	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	215	#print "PDSG -30,-30 - interp.get_A()", interp.get_A()
	216	answer = [ [ 0.0, 0.0, 0.0, 0.,
	217	0., 0.]]
	218	assert allclose(interp.get_A(), answer)
	219
	220
	221	#point outside of quad tree root cell
	222	interp.set_point_coordinates([[-70, -70]])
	223	#print "PDSG -70,-70 interp.get_A()", interp.get_A()
	224	answer = [ [ 0.0, 0.0, 0.0, 0.,
	225	0., 0. ]]
	226	assert allclose(interp.get_A(), answer)
	227
	228
	229	def test_datapoints_at_vertices(self):
	230	"""Test that data points coinciding with vertices yield a diagonal matrix
	231	"""
	232
	233	a = [0.0, 0.0]
	234	b = [0.0, 2.0]
	235	c = [2.0,0.0]
	236	points = [a, b, c]
	237	vertices = [ [1,0,2] ] #bac
	238
	239	data = points #Use data at vertices
	240
	241	interp = Interpolation(points, vertices, data)
	242	assert allclose(interp.get_A(), [[1., 0., 0.],
	243	[0., 1., 0.],
	244	[0., 0., 1.]])
	245
	246
	247
	248	def test_datapoints_on_edge_midpoints(self):
	249	"""Try datapoints midway on edges -
	250	each point should affect two matrix entries equally
	251	"""
	252
	253	a = [0.0, 0.0]
	254	b = [0.0, 2.0]
	255	c = [2.0,0.0]
	256	points = [a, b, c]
	257	vertices = [ [1,0,2] ] #bac
	258
	259	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	260
	261	interp = Interpolation(points, vertices, data)
	262
	263	assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
	264	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	265	[0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
	266
	267
	268	def test_datapoints_on_edges(self):
	269	"""Try datapoints on edges -
	270	each point should affect two matrix entries in proportion
	271	"""
	272
	273	a = [0.0, 0.0]
	274	b = [0.0, 2.0]
	275	c = [2.0,0.0]
	276	points = [a, b, c]
	277	vertices = [ [1,0,2] ] #bac
	278
	279	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	280
	281	interp = Interpolation(points, vertices, data)
	282
	283	assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
	284	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	285	[0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
	286
	287	def test_arbitrary_datapoints(self):
	288	"""Try arbitrary datapoints
	289	"""
	290
	291	from Numeric import sum
	292
	293	a = [0.0, 0.0]
	294	b = [0.0, 2.0]
	295	c = [2.0,0.0]
	296	points = [a, b, c]
	297	vertices = [ [1,0,2] ] #bac
	298
	299	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	300
	301	interp = Interpolation(points, vertices, data)
	302	#print "interp.get_A()", interp.get_A()
	303	assert allclose(sum(interp.get_A(), axis=1), 1.0)
	304
	305	def test_arbitrary_datapoints_some_outside(self):
	306	"""Try arbitrary datapoints one outside the triangle.
	307	That one should be ignored
	308	"""
	309
	310	from Numeric import sum
	311
	312	a = [0.0, 0.0]
	313	b = [0.0, 2.0]
	314	c = [2.0,0.0]
	315	points = [a, b, c]
	316	vertices = [ [1,0,2] ] #bac
	317
	318	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
	319
	320
	321	interp = Interpolation(points, vertices, data, precrop = True)
	322	assert allclose(sum(interp.get_A(), axis=1), 1.0)
	323
	324	interp = Interpolation(points, vertices, data, precrop = False)
	325	assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
	326
	327
	328
	329	# this causes a memory error in scipy.sparse
	330	def test_more_triangles(self):
	331
	332	a = [-1.0, 0.0]
	333	b = [3.0, 4.0]
	334	c = [4.0,1.0]
	335	d = [-3.0, 2.0] #3
	336	e = [-1.0,-2.0]
	337	f = [1.0, -2.0] #5
	338
	339	points = [a, b, c, d,e,f]
	340	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
	341
	342	#Data points
	343	data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
	344	interp = Interpolation(points, triangles, data_points)
	345
	346	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
	347	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	348	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	349	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	350	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	351	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
	352
	353
	354	A = interp.get_A()
	355	for i in range(A.shape[0]):
	356	for j in range(A.shape[1]):
	357	if not allclose(A[i,j], answer[i][j]):
	358	print i,j,':',A[i,j], answer[i][j]
	359
	360
	361	assert allclose(interp.get_A(), answer)
	362
	363
	364
	365
	366	def test_smooth_attributes_to_mesh(self):
	367	a = [0.0, 0.0]
	368	b = [0.0, 5.0]
	369	c = [5.0, 0.0]
	370	points = [a, b, c]
	371	triangles = [ [1,0,2] ] #bac
	372
	373	d1 = [1.0, 1.0]
	374	d2 = [1.0, 3.0]
	375	d3 = [3.0,1.0]
	376	z1 = 2
	377	z2 = 4
	378	z3 = 4
	379	data_coords = [d1, d2, d3]
	380
	381	interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
	382	z = [z1, z2, z3]
	383	f = interp.fit(z)
	384	answer = [0, 5., 5.]
	385
	386	#print "f\n",f
	387	#print "answer\n",answer
	388
	389	assert allclose(f, answer, atol=1e-7)
	390
	391
	392	def test_smooth_att_to_meshII(self):
	393
	394	a = [0.0, 0.0]
	395	b = [0.0, 5.0]
	396	c = [5.0, 0.0]
	397	points = [a, b, c]
	398	triangles = [ [1,0,2] ] #bac
	399
	400	d1 = [1.0, 1.0]
	401	d2 = [1.0, 2.0]
	402	d3 = [3.0,1.0]
	403	data_coords = [d1, d2, d3]
	404	z = linear_function(data_coords)
	405	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
	406	f = interp.fit(z)
	407	answer = linear_function(points)
	408
	409	assert allclose(f, answer)
	410
	411	def test_smooth_attributes_to_meshIII(self):
	412
	413	a = [-1.0, 0.0]
	414	b = [3.0, 4.0]
	415	c = [4.0,1.0]
	416	d = [-3.0, 2.0] #3
	417	e = [-1.0,-2.0]
	418	f = [1.0, -2.0] #5
	419
	420	vertices = [a, b, c, d,e,f]
	421	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	422
	423	point_coords = [[-2.0, 2.0],
	424	[-1.0, 1.0],
	425	[0.0,2.0],
	426	[1.0, 1.0],
	427	[2.0, 1.0],
	428	[0.0,0.0],
	429	[1.0, 0.0],
	430	[0.0, -1.0],
	431	[-0.2,-0.5],
	432	[-0.9, -1.5],
	433	[0.5, -1.9],
	434	[3.0,1.0]]
	435
	436	z = linear_function(point_coords)
	437	interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
	438
	439	#print 'z',z
	440	f = interp.fit(z)
	441	answer = linear_function(vertices)
	442	#print "f\n",f
	443	#print "answer\n",answer
	444	assert allclose(f, answer)
	445
	446
	447	def test_smooth_attributes_to_meshIV(self):
	448	""" Testing 2 attributes smoothed to the mesh
	449	"""
	450
	451	a = [0.0, 0.0]
	452	b = [0.0, 5.0]
	453	c = [5.0, 0.0]
	454	points = [a, b, c]
	455	triangles = [ [1,0,2] ] #bac
	456
	457	d1 = [1.0, 1.0]
	458	d2 = [1.0, 3.0]
	459	d3 = [3.0, 1.0]
	460	z1 = [2, 4]
	461	z2 = [4, 8]
	462	z3 = [4, 8]
	463	data_coords = [d1, d2, d3]
	464
	465	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
	466	z = [z1, z2, z3]
	467	f = interp.fit_points(z)
	468	answer = [[0,0], [5., 10.], [5., 10.]]
	469	assert allclose(f, answer)
	470
	471	def test_interpolate_attributes_to_points(self):
	472	v0 = [0.0, 0.0]
	473	v1 = [0.0, 5.0]
	474	v2 = [5.0, 0.0]
	475
	476	vertices = [v0, v1, v2]
	477	triangles = [ [1,0,2] ] #bac
	478
	479	d0 = [1.0, 1.0]
	480	d1 = [1.0, 2.0]
	481	d2 = [3.0, 1.0]
	482	point_coords = [ d0, d1, d2]
	483
	484	interp = Interpolation(vertices, triangles, point_coords)
	485	f = linear_function(vertices)
[2750]	486	#z = interp.interpolate(f)
[2526]	487	answer = linear_function(point_coords)
	488
	489
[2750]	490	#assert allclose(z, answer)
[2526]	491
	492
	493	def test_interpolate_attributes_to_points_interp_only(self):
	494	v0 = [0.0, 0.0]
	495	v1 = [0.0, 5.0]
	496	v2 = [5.0, 0.0]
	497
	498	vertices = [v0, v1, v2]
	499	triangles = [ [1,0,2] ] #bac
	500
	501	d0 = [1.0, 1.0]
	502	d1 = [1.0, 2.0]
	503	d2 = [3.0, 1.0]
	504	point_coords = [ d0, d1, d2]
	505
	506	interp = Interpolation(vertices, triangles, point_coords,
	507	interp_only = True)
	508
	509	f = linear_function(vertices)
[2750]	510	#z = interp.interpolate(f)
[2526]	511	answer = linear_function(point_coords)
	512	#print "answer", answer
	513	#print "z", z
	514
[2750]	515	#assert allclose(z, answer)
[2526]	516
	517	def test_interpolate_attributes_to_pointsII(self):
	518	a = [-1.0, 0.0]
	519	b = [3.0, 4.0]
	520	c = [4.0, 1.0]
	521	d = [-3.0, 2.0] #3
	522	e = [-1.0, -2.0]
	523	f = [1.0, -2.0] #5
	524
	525	vertices = [a, b, c, d,e,f]
	526	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	527
	528
	529	point_coords = [[-2.0, 2.0],
	530	[-1.0, 1.0],
	531	[0.0, 2.0],
	532	[1.0, 1.0],
	533	[2.0, 1.0],
	534	[0.0, 0.0],
	535	[1.0, 0.0],
	536	[0.0, -1.0],
	537	[-0.2, -0.5],
	538	[-0.9, -1.5],
	539	[0.5, -1.9],
	540	[3.0, 1.0]]
	541
	542	interp = Interpolation(vertices, triangles, point_coords)
	543	f = linear_function(vertices)
[2750]	544	#z = interp.interpolate(f)
[2526]	545	answer = linear_function(point_coords)
	546	#print "z",z
	547	#print "answer",answer
[2750]	548	#assert allclose(z, answer)
[2526]	549
	550	def test_interpolate_attributes_to_pointsIII(self):
	551	"""Test linear interpolation of known values at vertices to
	552	new points inside a triangle
	553	"""
	554	a = [0.0, 0.0]
	555	b = [0.0, 5.0]
	556	c = [5.0, 0.0]
	557	d = [5.0, 5.0]
	558
	559	vertices = [a, b, c, d]
	560	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	561
	562	#Points within triangle 1
	563	d0 = [1.0, 1.0]
	564	d1 = [1.0, 2.0]
	565	d2 = [3.0, 1.0]
	566
	567	#Point within triangle 2
	568	d3 = [4.0, 3.0]
	569
	570	#Points on common edge
	571	d4 = [2.5, 2.5]
	572	d5 = [4.0, 1.0]
	573
	574	#Point on common vertex
	575	d6 = [0., 5.]
	576
	577	point_coords = [d0, d1, d2, d3, d4, d5, d6]
	578
	579	interp = Interpolation(vertices, triangles, point_coords)
	580
	581	#Known values at vertices
	582	#Functions are x+y, x+2y, 2x+y, x-y-5
	583	f = [ [0., 0., 0., -5.], # (0,0)
	584	[5., 10., 5., -10.], # (0,5)
	585	[5., 5., 10.0, 0.], # (5,0)
	586	[10., 15., 15., -5.]] # (5,5)
	587
[2750]	588	#z = interp.interpolate(f)
[2526]	589	answer = [ [2., 3., 3., -5.], # (1,1)
	590	[3., 5., 4., -6.], # (1,2)
	591	[4., 5., 7., -3.], # (3,1)
	592	[7., 10., 11., -4.], # (4,3)
	593	[5., 7.5, 7.5, -5.], # (2.5, 2.5)
	594	[5., 6., 9., -2.], # (4,1)
	595	[5., 10., 5., -10.]] # (0,5)
	596
	597	#print "***********"
	598	#print "z",z
	599	#print "answer",answer
	600	#print "***********"
	601
	602	#Should an error message be returned if points are outside
	603	# of the mesh? Not currently.
	604
[2750]	605	#assert allclose(z, answer)
[2526]	606
	607
	608	def test_interpolate_point_outside_of_mesh(self):
	609	"""Test linear interpolation of known values at vertices to
	610	new points inside a triangle
	611	"""
	612	a = [0.0, 0.0]
	613	b = [0.0, 5.0]
	614	c = [5.0, 0.0]
	615	d = [5.0, 5.0]
	616
	617	vertices = [a, b, c, d]
	618	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	619
	620	#Far away point
	621	d7 = [-1., -1.]
	622
	623	point_coords = [ d7]
	624
	625	interp = Interpolation(vertices, triangles, point_coords)
	626
	627	#Known values at vertices
	628	#Functions are x+y, x+2y, 2x+y, x-y-5
	629	f = [ [0., 0., 0., -5.], # (0,0)
	630	[5., 10., 5., -10.], # (0,5)
	631	[5., 5., 10.0, 0.], # (5,0)
	632	[10., 15., 15., -5.]] # (5,5)
	633
	634
[2750]	635	#z = interp.interpolate(f)
	636	#answer = [ [0., 0., 0., 0.]] # (-1,-1)
	637
[2526]	638	#print "***********"
	639	#print "z",z
	640	#print "answer",answer
	641	#print "***********"
	642
	643	#Should an error message be returned if points are outside
	644	# of the mesh? Not currently.
	645
[2750]	646	#assert allclose(z, answer)
[2526]	647
	648	def test_interpolate_attributes_to_pointsIV(self):
	649	a = [-1.0, 0.0]
	650	b = [3.0, 4.0]
	651	c = [4.0, 1.0]
	652	d = [-3.0, 2.0] #3
	653	e = [-1.0, -2.0]
	654	f = [1.0, -2.0] #5
	655
	656	vertices = [a, b, c, d,e,f]
	657	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	658
	659
	660	point_coords = [[-2.0, 2.0],
	661	[-1.0, 1.0],
	662	[0.0, 2.0],
	663	[1.0, 1.0],
	664	[2.0, 1.0],
	665	[0.0, 0.0],
	666	[1.0, 0.0],
	667	[0.0, -1.0],
	668	[-0.2, -0.5],
	669	[-0.9, -1.5],
	670	[0.5, -1.9],
	671	[3.0, 1.0]]
	672
	673	interp = Interpolation(vertices, triangles, point_coords)
	674	f = array([linear_function(vertices),2*linear_function(vertices) ])
	675	f = transpose(f)
	676	#print "f",f
[2750]	677	#z = interp.interpolate(f)
[2526]	678	answer = [linear_function(point_coords),
	679	2*linear_function(point_coords) ]
	680	answer = transpose(answer)
	681	#print "z",z
	682	#print "answer",answer
[2750]	683	#assert allclose(z, answer)
[2526]	684
	685	def test_smooth_attributes_to_mesh_function(self):
	686	""" Testing 2 attributes smoothed to the mesh
	687	"""
	688
	689	a = [0.0, 0.0]
	690	b = [0.0, 5.0]
	691	c = [5.0, 0.0]
	692	points = [a, b, c]
	693	triangles = [ [1,0,2] ] #bac
	694
	695	d1 = [1.0, 1.0]
	696	d2 = [1.0, 3.0]
	697	d3 = [3.0, 1.0]
	698	z1 = [2, 4]
	699	z2 = [4, 8]
	700	z3 = [4, 8]
	701	data_coords = [d1, d2, d3]
	702	z = [z1, z2, z3]
	703
	704	f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
	705	answer = [[0, 0], [5., 10.], [5., 10.]]
	706
	707	assert allclose(f, answer)
	708
	709
	710
	711	def test_pts2rectangular(self):
	712
	713	import time, os
	714	FN = 'xyatest' + str(time.time()) + '.xya'
	715	fid = open(FN, 'w')
	716	fid.write(' %s \n' %('elevation'))
	717	fid.write('%f %f %f\n' %(1,1,2) )
	718	fid.write('%f %f %f\n' %(1,3,4) )
	719	fid.write('%f %f %f\n' %(3,1,4) )
	720	fid.close()
	721
	722	points, triangles, boundary, attributes =\
	723	pts2rectangular(FN, 4, 4)
	724
	725
	726	data_coords = [ [1,1], [1,3], [3,1] ]
	727	z = [2, 4, 4]
	728
[2585]	729	ref = fit_to_mesh(points, triangles, data_coords, z, verbose=False)
[2526]	730
	731	#print attributes
	732	#print ref
	733	assert allclose(attributes, ref)
	734
	735	os.remove(FN)
	736
	737
	738	#Tests of smoothing matrix
	739	def test_smoothing_matrix_one_triangle(self):
	740	from Numeric import dot
	741	a = [0.0, 0.0]
	742	b = [0.0, 2.0]
	743	c = [2.0,0.0]
	744	points = [a, b, c]
	745
	746	vertices = [ [1,0,2] ] #bac
	747
	748	interp = Interpolation(points, vertices)
	749
	750	assert allclose(interp.get_D(), [[1, -0.5, -0.5],
	751	[-0.5, 0.5, 0],
	752	[-0.5, 0, 0.5]])
	753
	754	#Define f(x,y) = x
	755	f = array([0,0,2]) #Value at global vertex 2
	756
	757	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	758	# int 1 dx dy = area = 2
	759	assert dot(dot(f, interp.get_D()), f) == 2
	760
	761	#Define f(x,y) = y
	762	f = array([0,2,0]) #Value at global vertex 1
	763
	764	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	765	# int 1 dx dy = area = 2
	766	assert dot(dot(f, interp.get_D()), f) == 2
	767
	768	#Define f(x,y) = x+y
	769	f = array([0,2,2]) #Values at global vertex 1 and 2
	770
	771	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	772	# int 2 dx dy = 2*area = 4
	773	assert dot(dot(f, interp.get_D()), f) == 4
	774
	775
	776
	777	def test_smoothing_matrix_more_triangles(self):
	778	from Numeric import dot
	779
	780	a = [0.0, 0.0]
	781	b = [0.0, 2.0]
	782	c = [2.0,0.0]
	783	d = [0.0, 4.0]
	784	e = [2.0, 2.0]
	785	f = [4.0,0.0]
	786
	787	points = [a, b, c, d, e, f]
	788	#bac, bce, ecf, dbe, daf, dae
	789	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	790
	791	interp = Interpolation(points, vertices)
	792
	793
	794	#assert allclose(interp.get_D(), [[1, -0.5, -0.5],
	795	# [-0.5, 0.5, 0],
	796	# [-0.5, 0, 0.5]])
	797
	798	#Define f(x,y) = x
	799	f = array([0,0,2,0,2,4]) #f evaluated at points a-f
	800
	801	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	802	# int 1 dx dy = total area = 8
	803	assert dot(dot(f, interp.get_D()), f) == 8
	804
	805	#Define f(x,y) = y
	806	f = array([0,2,0,4,2,0]) #f evaluated at points a-f
	807
	808	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	809	# int 1 dx dy = area = 8
	810	assert dot(dot(f, interp.get_D()), f) == 8
	811
	812	#Define f(x,y) = x+y
	813	f = array([0,2,2,4,4,4]) #f evaluated at points a-f
	814
	815	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	816	# int 2 dx dy = 2*area = 16
	817	assert dot(dot(f, interp.get_D()), f) == 16
	818
	819
	820	def test_fit_and_interpolation(self):
	821	from mesh import Mesh
	822
	823	a = [0.0, 0.0]
	824	b = [0.0, 2.0]
	825	c = [2.0, 0.0]
	826	d = [0.0, 4.0]
	827	e = [2.0, 2.0]
	828	f = [4.0, 0.0]
	829
	830	points = [a, b, c, d, e, f]
	831	#bac, bce, ecf, dbe, daf, dae
	832	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	833
	834	#Get (enough) datapoints
	835	data_points = [[ 0.66666667, 0.66666667],
	836	[ 1.33333333, 1.33333333],
	837	[ 2.66666667, 0.66666667],
	838	[ 0.66666667, 2.66666667],
	839	[ 0.0, 1.0],
	840	[ 0.0, 3.0],
	841	[ 1.0, 0.0],
	842	[ 1.0, 1.0],
	843	[ 1.0, 2.0],
	844	[ 1.0, 3.0],
	845	[ 2.0, 1.0],
	846	[ 3.0, 0.0],
	847	[ 3.0, 1.0]]
	848
	849	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	850
	851	z = linear_function(data_points)
	852	answer = linear_function(points)
	853
	854	f = interp.fit(z)
	855
	856	#print "f",f
	857	#print "answer",answer
	858	assert allclose(f, answer)
	859
	860	#Map back
[2750]	861	#z1 = interp.interpolate(f)
[2526]	862	#print "z1\n", z1
	863	#print "z\n",z
[2750]	864	#assert allclose(z, z1)
[2526]	865
	866
	867	def test_smoothing_and_interpolation(self):
	868
	869	a = [0.0, 0.0]
	870	b = [0.0, 2.0]
	871	c = [2.0, 0.0]
	872	d = [0.0, 4.0]
	873	e = [2.0, 2.0]
	874	f = [4.0, 0.0]
	875
	876	points = [a, b, c, d, e, f]
	877	#bac, bce, ecf, dbe, daf, dae
	878	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	879
	880	#Get (too few!) datapoints
	881	data_points = [[ 0.66666667, 0.66666667],
	882	[ 1.33333333, 1.33333333],
	883	[ 2.66666667, 0.66666667],
	884	[ 0.66666667, 2.66666667]]
	885
	886	z = linear_function(data_points)
	887	answer = linear_function(points)
	888
	889	#Make interpolator with too few data points and no smoothing
	890	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	891	#Must raise an exception
	892	try:
	893	f = interp.fit(z)
	894	except:
	895	pass
	896
	897	#Now try with smoothing parameter
	898	interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
	899
	900	f = interp.fit(z)
	901	#f will be different from answer due to smoothing
	902	assert allclose(f, answer,atol=5)
	903
	904	#Map back
[2750]	905	#z1 = interp.interpolate(f)
	906	#assert allclose(z, z1)
[2526]	907
	908
	909
	910	def test_fit_and_interpolation_with_new_points(self):
	911	"""Fit a surface to one set of points. Then interpolate that surface
	912	using another set of points.
	913	"""
	914	from mesh import Mesh
	915
	916
	917	#Setup mesh used to represent fitted function
	918	a = [0.0, 0.0]
	919	b = [0.0, 2.0]
	920	c = [2.0, 0.0]
	921	d = [0.0, 4.0]
	922	e = [2.0, 2.0]
	923	f = [4.0, 0.0]
	924
	925	points = [a, b, c, d, e, f]
	926	#bac, bce, ecf, dbe, daf, dae
	927	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	928
	929	#Datapoints to fit from
	930	data_points1 = [[ 0.66666667, 0.66666667],
	931	[ 1.33333333, 1.33333333],
	932	[ 2.66666667, 0.66666667],
	933	[ 0.66666667, 2.66666667],
	934	[ 0.0, 1.0],
	935	[ 0.0, 3.0],
	936	[ 1.0, 0.0],
	937	[ 1.0, 1.0],
	938	[ 15, -17], #Outside mesh
	939	[ 1.0, 2.0],
	940	[ 1.0, 3.0],
	941	[ 2.0, 1.0],
	942	[ 3.0, 0.0],
	943	[ 3.0, 1.0]]
	944
	945	#Fit surface to mesh
	946	interp = Interpolation(points, triangles, data_points1, alpha=0.0,
[2585]	947	precrop = True, verbose=False)
[2526]	948	z = linear_function(data_points1) #Example z-values
	949	f = interp.fit(z) #Fitted values at vertices
	950
	951
	952
	953	#New datapoints where interpolated values are sought
	954	data_points2 = [[ 0.0, 0.0],
	955	[ 0.5, 0.5],
	956	[ 0.7, 0.7],
	957	[-13, 65], #Outside
	958	[ 1.0, 0.5],
	959	[ 2.0, 0.4],
	960	[ 2.8, 1.2]]
	961
	962
	963
	964	#Build new A matrix based on new points (without precrop)
	965	interp.build_interpolation_matrix_A(data_points2, precrop = False)
	966
	967	#Interpolate using fitted surface
[2750]	968	#z1 = interp.interpolate(f)
[2526]	969
	970	#import Numeric
	971	#data_points2 = Numeric.take(data_points2, interp.point_indices)
	972
	973	#Desired result (OK for points inside)
	974
	975	answer = linear_function(data_points2)
	976	import Numeric
[2750]	977	#z1 = Numeric.take(z1, [0,1,2,4,5,6])
[2526]	978	answer = Numeric.take(answer, [0,1,2,4,5,6])
[2750]	979	#assert allclose(z1, answer)
[2526]	980
	981	#Build new A matrix based on new points (with precrop)
	982	interp.build_interpolation_matrix_A(data_points2, precrop = True)
	983
	984	#Interpolate using fitted surface
[2750]	985	#z1 = interp.interpolate(f)
[2526]	986
	987	import Numeric
	988	data_points2 = Numeric.take(data_points2, interp.point_indices)
	989
	990	#Desired result
	991	answer = linear_function(data_points2)
[2750]	992	#assert allclose(z1, answer)
[2526]	993
	994
	995
	996
	997
	998
	999	def test_interpolation_from_discontinuous_vertex_values(self):
	1000	"""test_interpolation_from_discontinuous_vertex_values.
	1001	This will test the format used internally in pyvolution and also
	1002	interpolation from sww files
	1003	"""
	1004
	1005	from mesh import Mesh
	1006
	1007
	1008	#Setup mesh used to represent discontinuous function
	1009	a = [0.0, 0.0]
	1010	b = [0.0, 2.0]
	1011	c = [2.0, 0.0]
	1012	d = [0.0, 4.0]
	1013	e = [2.0, 2.0]
	1014	f = [4.0, 0.0]
	1015
	1016	points = [b, a, c,
	1017	b, c, e,
	1018	e, c, f,
	1019	d, b, e]
	1020
	1021	#bac, bce, ecf, dbe
	1022	triangles = [[0,1,2], [3,4,5], [6,7,8], [9,10,11]]
	1023
	1024
	1025	vertex_values = [0.,0.,0.,1.,1.,1.,2.,2.,2.,7.,3.,3.]
	1026
	1027
	1028
	1029	#New datapoints where interpolated values are sought
	1030	data_points = [[0.0, 0.0], #T0
	1031	[0.5, 0.5], #T0
	1032	[1.5, 1.5], #T1
	1033	[2.5, 0.5], #T2
	1034	[0.0, 3.0], #T3
	1035	[1.0, 2.0], #In between T1 and T3 (T1 is used) FIXME?
	1036	[2.0, 1.0], #In between T1 and T2 (T1 is used) FIXME?
	1037	[1.0, 1.0]] #In between T1 and T0 (T0 is used) FIXME?
	1038
	1039
	1040
	1041
	1042	#Build interpolation matrix
	1043	interp = Interpolation(points, triangles, data_points)
	1044	#, alpha=0.0, precrop = True)
	1045
	1046	#print interp.A.todense()
	1047	#print vertex_values
	1048
	1049	#Interpolate using fitted surface
[2750]	1050	#z = interp.interpolate(vertex_values)
[2526]	1051
	1052	#print z
	1053
[2750]	1054	#assert allclose(z, [0,0,1,2,5,1,1,0])
[2526]	1055
	1056
	1057
	1058
	1059	def test_interpolation_function_time_only(self):
	1060	"""Test spatio-temporal interpolation
	1061	Test that spatio temporal function performs the correct
	1062	interpolations in both time and space
	1063	"""
	1064
	1065
	1066	#Three timesteps
	1067	time = [1.0, 5.0, 6.0]
	1068
	1069
	1070	#One quantity
	1071	Q = zeros( (3,6), Float )
	1072
	1073	#Linear in time and space
	1074	a = [0.0, 0.0]
	1075	b = [0.0, 2.0]
	1076	c = [2.0, 0.0]
	1077	d = [0.0, 4.0]
	1078	e = [2.0, 2.0]
	1079	f = [4.0, 0.0]
	1080
	1081	points = [a, b, c, d, e, f]
	1082
	1083	for i, t in enumerate(time):
	1084	Q[i, :] = t*linear_function(points)
	1085
	1086
	1087	#Check basic interpolation of one quantity using averaging
	1088	#(no interpolation points or spatial info)
	1089	from utilities.numerical_tools import mean
	1090	I = Interpolation_function(time, [mean(Q[0,:]),
	1091	mean(Q[1,:]),
	1092	mean(Q[2,:])])
	1093
	1094
	1095
	1096	#Check temporal interpolation
	1097	for i in [0,1,2]:
	1098	assert allclose(I(time[i]), mean(Q[i,:]))
	1099
	1100	#Midway
	1101	assert allclose(I( (time[0] + time[1])/2 ),
	1102	(I(time[0]) + I(time[1]))/2 )
	1103
	1104	assert allclose(I( (time[1] + time[2])/2 ),
	1105	(I(time[1]) + I(time[2]))/2 )
	1106
	1107	assert allclose(I( (time[0] + time[2])/2 ),
	1108	(I(time[0]) + I(time[2]))/2 )
	1109
	1110	#1/3
	1111	assert allclose(I( (time[0] + time[2])/3 ),
	1112	(I(time[0]) + I(time[2]))/3 )
	1113
	1114
	1115	#Out of bounds checks
	1116	try:
	1117	I(time[0]-1)
	1118	except:
	1119	pass
	1120	else:
	1121	raise 'Should raise exception'
	1122
	1123	try:
	1124	I(time[-1]+1)
	1125	except:
	1126	pass
	1127	else:
	1128	raise 'Should raise exception'
	1129
	1130
	1131
	1132
[2750]	1133	def interpolation_test_interpolation_function_spatial_only(self):
[2526]	1134	"""Test spatio-temporal interpolation with constant time
	1135	"""
	1136
	1137	#Three timesteps
	1138	time = [1.0, 5.0, 6.0]
	1139
	1140
	1141	#Setup mesh used to represent fitted function
	1142	a = [0.0, 0.0]
	1143	b = [0.0, 2.0]
	1144	c = [2.0, 0.0]
	1145	d = [0.0, 4.0]
	1146	e = [2.0, 2.0]
	1147	f = [4.0, 0.0]
	1148
	1149	points = [a, b, c, d, e, f]
	1150	#bac, bce, ecf, dbe
	1151	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1152
	1153
	1154	#New datapoints where interpolated values are sought
	1155	interpolation_points = [[ 0.0, 0.0],
	1156	[ 0.5, 0.5],
	1157	[ 0.7, 0.7],
	1158	[ 1.0, 0.5],
	1159	[ 2.0, 0.4],
	1160	[ 2.8, 1.2]]
	1161
	1162
	1163	#One quantity linear in space
	1164	Q = linear_function(points)
	1165
	1166
	1167	#Check interpolation of one quantity using interpolaton points
	1168	I = Interpolation_function(time, Q,
	1169	vertex_coordinates = points,
	1170	triangles = triangles,
	1171	interpolation_points = interpolation_points,
	1172	verbose = False)
	1173
	1174
	1175	answer = linear_function(interpolation_points)
	1176
	1177	t = time[0]
	1178	for j in range(50): #t in [1, 6]
	1179	for id in range(len(interpolation_points)):
	1180	assert allclose(I(t, id), answer[id])
	1181
	1182	t += 0.1
	1183
	1184
	1185	try:
	1186	I(1)
	1187	except:
	1188	pass
	1189	else:
	1190	raise 'Should raise exception'
	1191
	1192
	1193
[2750]	1194	def interpolation_test_interpolation_function(self):
[2526]	1195	"""Test spatio-temporal interpolation
	1196	Test that spatio temporal function performs the correct
	1197	interpolations in both time and space
	1198	"""
	1199
	1200
	1201	#Three timesteps
	1202	time = [1.0, 5.0, 6.0]
	1203
	1204
	1205	#Setup mesh used to represent fitted function
	1206	a = [0.0, 0.0]
	1207	b = [0.0, 2.0]
	1208	c = [2.0, 0.0]
	1209	d = [0.0, 4.0]
	1210	e = [2.0, 2.0]
	1211	f = [4.0, 0.0]
	1212
	1213	points = [a, b, c, d, e, f]
	1214	#bac, bce, ecf, dbe
	1215	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1216
	1217
	1218	#New datapoints where interpolated values are sought
	1219	interpolation_points = [[ 0.0, 0.0],
	1220	[ 0.5, 0.5],
	1221	[ 0.7, 0.7],
	1222	[ 1.0, 0.5],
	1223	[ 2.0, 0.4],
	1224	[ 2.8, 1.2]]
	1225
	1226
	1227	#One quantity
	1228	Q = zeros( (3,6), Float )
	1229
	1230	#Linear in time and space
	1231	for i, t in enumerate(time):
	1232	Q[i, :] = t*linear_function(points)
	1233
	1234
	1235	#Check interpolation of one quantity using interpolaton points)
	1236	I = Interpolation_function(time, Q,
	1237	vertex_coordinates = points,
	1238	triangles = triangles,
	1239	interpolation_points = interpolation_points,
	1240	verbose = False)
	1241
	1242
	1243	answer = linear_function(interpolation_points)
	1244
	1245	t = time[0]
	1246	for j in range(50): #t in [1, 6]
	1247	for id in range(len(interpolation_points)):
	1248	assert allclose(I(t, id), t*answer[id])
	1249
	1250	t += 0.1
	1251
	1252
	1253	try:
	1254	I(1)
	1255	except:
	1256	pass
	1257	else:
	1258	raise 'Should raise exception'
	1259
	1260	#
	1261	#interpolation_points = [[ 0.0, 0.0],
	1262	# [ 0.5, 0.5],
	1263	# [ 0.7, 0.7],
	1264	# [-13, 65], #Outside
	1265	# [ 1.0, 0.5],
	1266	# [ 2.0, 0.4],
	1267	# [ 2.8, 1.2]]
	1268	#
	1269	#try:
	1270	# I = Interpolation_function(time, Q,
	1271	# vertex_coordinates = points,
	1272	# triangles = triangles,
	1273	# interpolation_points = interpolation_points,
	1274	# verbose = False)
	1275	#except:
	1276	# pass
	1277	#else:
	1278	# raise 'Should raise exception'
	1279
	1280
	1281
	1282
	1283
	1284	def test_fit_and_interpolation_with_different_origins(self):
	1285	"""Fit a surface to one set of points. Then interpolate that surface
	1286	using another set of points.
	1287	This test tests situtaion where points and mesh belong to a different
	1288	coordinate system as defined by origin.
	1289	"""
	1290	from mesh import Mesh
	1291
	1292	#Setup mesh used to represent fitted function
	1293	a = [0.0, 0.0]
	1294	b = [0.0, 2.0]
	1295	c = [2.0, 0.0]
	1296	d = [0.0, 4.0]
	1297	e = [2.0, 2.0]
	1298	f = [4.0, 0.0]
	1299
	1300	points = [a, b, c, d, e, f]
	1301	#bac, bce, ecf, dbe, daf, dae
	1302	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1303
	1304	#Datapoints to fit from
	1305	data_points1 = [[ 0.66666667, 0.66666667],
	1306	[ 1.33333333, 1.33333333],
	1307	[ 2.66666667, 0.66666667],
	1308	[ 0.66666667, 2.66666667],
	1309	[ 0.0, 1.0],
	1310	[ 0.0, 3.0],
	1311	[ 1.0, 0.0],
	1312	[ 1.0, 1.0],
	1313	[ 1.0, 2.0],
	1314	[ 1.0, 3.0],
	1315	[ 2.0, 1.0],
	1316	[ 3.0, 0.0],
	1317	[ 3.0, 1.0]]
	1318
	1319
	1320	#First check that things are OK when using same origin
	1321	mesh_origin = (56, 290000, 618000) #zone, easting, northing
	1322	data_origin = (56, 290000, 618000) #zone, easting, northing
	1323
	1324
	1325	#Fit surface to mesh
	1326	interp = Interpolation(points, triangles, data_points1,
	1327	alpha=0.0,
	1328	data_origin = data_origin,
	1329	mesh_origin = mesh_origin)
	1330
	1331	z = linear_function(data_points1) #Example z-values
	1332	f = interp.fit(z) #Fitted values at vertices
	1333
	1334
	1335	#New datapoints where interpolated values are sought
	1336	data_points2 = [[ 0.0, 0.0],
	1337	[ 0.5, 0.5],
	1338	[ 0.7, 0.7],
	1339	[ 1.0, 0.5],
	1340	[ 2.0, 0.4],
	1341	[ 2.8, 1.2]]
	1342
	1343
	1344	#Build new A matrix based on new points
	1345	interp.build_interpolation_matrix_A(data_points2)
	1346
	1347	#Interpolate using fitted surface
[2750]	1348	#z1 = interp.interpolate(f)
[2526]	1349
	1350	#Desired result
[2750]	1351	#answer = linear_function(data_points2)
	1352	#assert allclose(z1, answer)
[2526]	1353
	1354
	1355	##############################################
	1356
	1357	#Then check situation where points are relative to a different
	1358	#origin (same zone, though, until we figure that out (FIXME))
	1359
	1360	mesh_origin = (56, 290000, 618000) #zone, easting, northing
	1361	data_origin = (56, 10000, 10000) #zone, easting, northing
	1362
	1363	#Shift datapoints according to new origin
	1364
	1365	for k in range(len(data_points1)):
	1366	data_points1[k][0] += mesh_origin[1] - data_origin[1]
	1367	data_points1[k][1] += mesh_origin[2] - data_origin[2]
	1368
	1369	for k in range(len(data_points2)):
	1370	data_points2[k][0] += mesh_origin[1] - data_origin[1]
	1371	data_points2[k][1] += mesh_origin[2] - data_origin[2]
	1372
	1373
	1374
	1375	#Fit surface to mesh
	1376	interp = Interpolation(points, triangles, data_points1,
	1377	alpha=0.0,
	1378	data_origin = data_origin,
	1379	mesh_origin = mesh_origin)
	1380
	1381	f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
	1382
	1383	assert allclose(f,f1), 'Fit should have been unaltered'
	1384
	1385
	1386	#Build new A matrix based on new points
	1387	interp.build_interpolation_matrix_A(data_points2)
	1388
	1389	#Interpolate using fitted surface
[2750]	1390	#z1 = interp.interpolate(f)
	1391	#assert allclose(z1, answer)
[2526]	1392
	1393
	1394	#########################################################
	1395	#Finally try to relate data_points2 to new origin without
	1396	#rebuilding matrix
	1397
	1398	data_origin = (56, 2000, 2000) #zone, easting, northing
	1399	for k in range(len(data_points2)):
	1400	data_points2[k][0] += 8000
	1401	data_points2[k][1] += 8000
	1402
	1403	#Build new A matrix based on new points
	1404	interp.build_interpolation_matrix_A(data_points2,
	1405	data_origin = data_origin)
	1406
	1407	#Interpolate using fitted surface
[2750]	1408	#z1 = interp.interpolate(f)
	1409	#assert allclose(z1, answer)
[2526]	1410
	1411
	1412
	1413	def test_fit_to_mesh_w_georef(self):
	1414	"""Simple check that georef works at the fit_to_mesh level
	1415	"""
	1416
	1417	from coordinate_transforms.geo_reference import Geo_reference
	1418
	1419	#Mesh
	1420	vertex_coordinates = [[0.76, 0.76],
	1421	[0.76, 5.76],
	1422	[5.76, 0.76]]
	1423	triangles = [[0,2,1]]
	1424
	1425	mesh_geo = Geo_reference(56,-0.76,-0.76)
	1426
	1427
	1428	#Data
	1429	data_points = [[ 201.0, 401.0],
	1430	[ 201.0, 403.0],
	1431	[ 203.0, 401.0]]
	1432
	1433	z = [2, 4, 4]
	1434
	1435	data_geo = Geo_reference(56,-200,-400)
	1436
	1437	#Fit
	1438	zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z,
	1439	data_origin = data_geo.get_origin(),
	1440	mesh_origin = mesh_geo.get_origin(),
	1441	alpha = 0)
	1442	assert allclose( zz, [0,5,5] )
	1443
	1444
	1445	def test_fit_to_mesh_file(self):
	1446	from load_mesh.loadASCII import import_mesh_file, \
	1447	export_mesh_file
	1448	import tempfile
	1449	import os
	1450
	1451	# create a .tsh file, no user outline
	1452	mesh_dic = {}
	1453	mesh_dic['vertices'] = [[0.0, 0.0],
	1454	[0.0, 5.0],
	1455	[5.0, 0.0]]
	1456	mesh_dic['triangles'] = [[0, 2, 1]]
	1457	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1458	mesh_dic['triangle_tags'] = ['']
	1459	mesh_dic['vertex_attributes'] = [[], [], []]
	1460	mesh_dic['vertiex_attribute_titles'] = []
	1461	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1462	mesh_dic['segment_tags'] = ['external',
	1463	'external',
	1464	'external']
	1465	mesh_file = tempfile.mktemp(".tsh")
	1466	export_mesh_file(mesh_file,mesh_dic)
	1467
	1468	# create an .xya file
	1469	point_file = tempfile.mktemp(".xya")
	1470	fd = open(point_file,'w')
	1471	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1472	fd.close()
	1473
	1474	mesh_output_file = tempfile.mktemp(".tsh")
	1475	fit_to_mesh_file(mesh_file,
	1476	point_file,
	1477	mesh_output_file,
	1478	alpha = 0.0)
	1479	# load in the .tsh file we just wrote
	1480	mesh_dic = import_mesh_file(mesh_output_file)
	1481	#print "mesh_dic",mesh_dic
	1482	ans =[[0.0, 0.0],
	1483	[5.0, 10.0],
	1484	[5.0,10.0]]
	1485	assert allclose(mesh_dic['vertex_attributes'],ans)
	1486
	1487	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1488	['elevation','stage'],
	1489	'test_fit_to_mesh_file failed')
	1490
	1491	#clean up
	1492	os.remove(mesh_file)
	1493	os.remove(point_file)
	1494	os.remove(mesh_output_file)
	1495
	1496	def test_fit_to_mesh_file3(self):
	1497	from load_mesh.loadASCII import import_mesh_file, \
	1498	export_mesh_file
	1499	import tempfile
	1500	import os
	1501
	1502	# create a .tsh file, no user outline
	1503	mesh_dic = {}
	1504	mesh_dic['vertices'] = [[0.76, 0.76],
	1505	[0.76, 5.76],
	1506	[5.76, 0.76]]
	1507	mesh_dic['triangles'] = [[0, 2, 1]]
	1508	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1509	mesh_dic['triangle_tags'] = ['']
	1510	mesh_dic['vertex_attributes'] = [[], [], []]
	1511	mesh_dic['vertiex_attribute_titles'] = []
	1512	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1513	mesh_dic['segment_tags'] = ['external',
	1514	'external',
	1515	'external']
	1516	mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
	1517	mesh_file = tempfile.mktemp(".tsh")
	1518	export_mesh_file(mesh_file,mesh_dic)
	1519
	1520	# create an .xya file
	1521	point_file = tempfile.mktemp(".xya")
	1522	fd = open(point_file,'w')
	1523	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1524	fd.close()
	1525
	1526	mesh_output_file = tempfile.mktemp(".tsh")
	1527	fit_to_mesh_file(mesh_file,
	1528	point_file,
	1529	mesh_output_file,
	1530	alpha = 0.0)
	1531	# load in the .tsh file we just wrote
	1532	mesh_dic = import_mesh_file(mesh_output_file)
	1533	#print "mesh_dic",mesh_dic
	1534	ans =[[0.0, 0.0],
	1535	[5.0, 10.0],
	1536	[5.0,10.0]]
	1537	assert allclose(mesh_dic['vertex_attributes'],ans)
	1538
	1539	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1540	['elevation','stage'],
	1541	'test_fit_to_mesh_file failed')
	1542
	1543	#clean up
	1544	os.remove(mesh_file)
	1545	os.remove(point_file)
	1546	os.remove(mesh_output_file)
	1547
	1548	def test_fit_to_mesh_file4(self):
	1549	from load_mesh.loadASCII import import_mesh_file, \
	1550	export_mesh_file
	1551	import tempfile
	1552	import os
	1553
	1554	# create a .tsh file, no user outline
	1555	mesh_dic = {}
	1556	mesh_dic['vertices'] = [[0.76, 0.76],
	1557	[0.76, 5.76],
	1558	[5.76, 0.76]]
	1559	mesh_dic['triangles'] = [[0, 2, 1]]
	1560	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1561	mesh_dic['triangle_tags'] = ['']
	1562	mesh_dic['vertex_attributes'] = [[], [], []]
	1563	mesh_dic['vertiex_attribute_titles'] = []
	1564	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1565	mesh_dic['segment_tags'] = ['external',
	1566	'external',
	1567	'external']
	1568	mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
	1569	mesh_file = tempfile.mktemp(".tsh")
	1570	export_mesh_file(mesh_file,mesh_dic)
	1571
	1572	geo_ref = Geo_reference(56,-200,-400)
	1573	# create an .xya file
	1574	point_file = tempfile.mktemp(".xya")
	1575	fd = open(point_file,'w')
	1576	fd.write("elevation, stage \n 201.0, 401.0,2.,4 \n 201.0, 403.0,4,8 \n 203.0, 401.0,4.,8 \n")
	1577	geo_ref.write_ASCII(fd)
	1578	fd.close()
	1579
	1580	mesh_output_file = tempfile.mktemp(".tsh")
	1581	fit_to_mesh_file(mesh_file,
	1582	point_file,
	1583	mesh_output_file,
	1584	alpha = 0.0)
	1585	# load in the .tsh file we just wrote
	1586	mesh_dic = import_mesh_file(mesh_output_file)
	1587	#print "mesh_dic",mesh_dic
	1588	ans =[[0.0, 0.0],
	1589	[5.0, 10.0],
	1590	[5.0, 10.0]]
	1591	assert allclose(mesh_dic['vertex_attributes'],ans)
	1592
	1593	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1594	['elevation','stage'],
	1595	'test_fit_to_mesh_file failed')
	1596
	1597	#clean up
	1598	os.remove(mesh_file)
	1599	os.remove(point_file)
	1600	os.remove(mesh_output_file)
	1601
	1602	def test_fit_to_mesh_fileII(self):
	1603	from load_mesh.loadASCII import import_mesh_file, \
	1604	export_mesh_file
	1605	import tempfile
	1606	import os
	1607
	1608	# create a .tsh file, no user outline
	1609	mesh_dic = {}
	1610	mesh_dic['vertices'] = [[0.0, 0.0],
	1611	[0.0, 5.0],
	1612	[5.0, 0.0]]
	1613	mesh_dic['triangles'] = [[0, 2, 1]]
	1614	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1615	mesh_dic['triangle_tags'] = ['']
	1616	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1617	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1618	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1619	mesh_dic['segment_tags'] = ['external',
	1620	'external',
	1621	'external']
	1622	mesh_file = tempfile.mktemp(".tsh")
	1623	export_mesh_file(mesh_file,mesh_dic)
	1624
	1625	# create an .xya file
	1626	point_file = tempfile.mktemp(".xya")
	1627	fd = open(point_file,'w')
	1628	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1629	fd.close()
	1630
	1631	mesh_output_file = "new_triangle.tsh"
	1632	fit_to_mesh_file(mesh_file,
	1633	point_file,
	1634	mesh_output_file,
	1635	alpha = 0.0)
	1636	# load in the .tsh file we just wrote
	1637	mesh_dic = import_mesh_file(mesh_output_file)
	1638
	1639	assert allclose(mesh_dic['vertex_attributes'],
	1640	[[1.0, 2.0,0.0, 0.0],
	1641	[1.0, 2.0,5.0, 10.0],
	1642	[1.0, 2.0,5.0,10.0]])
	1643
	1644	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1645	['density', 'temp','elevation','stage'],
	1646	'test_fit_to_mesh_file failed')
	1647
	1648	#clean up
	1649	os.remove(mesh_file)
	1650	os.remove(mesh_output_file)
	1651	os.remove(point_file)
	1652
	1653	def test_fit_to_mesh_file_errors(self):
	1654	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1655	import tempfile
	1656	import os
	1657
	1658	# create a .tsh file, no user outline
	1659	mesh_dic = {}
	1660	mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
	1661	mesh_dic['triangles'] = [[0, 2, 1]]
	1662	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1663	mesh_dic['triangle_tags'] = ['']
	1664	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1665	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1666	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1667	mesh_dic['segment_tags'] = ['external', 'external','external']
	1668	mesh_file = tempfile.mktemp(".tsh")
	1669	export_mesh_file(mesh_file,mesh_dic)
	1670
	1671	# create an .xya file
	1672	point_file = tempfile.mktemp(".xya")
	1673	fd = open(point_file,'w')
	1674	fd.write("elevation stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1675	fd.close()
	1676
	1677	mesh_output_file = "new_triangle.tsh"
	1678	try:
	1679	fit_to_mesh_file(mesh_file, point_file,
	1680	mesh_output_file, display_errors = False)
	1681	except IOError:
	1682	pass
	1683	else:
	1684	#self.failUnless(0 ==1, 'Bad file did not raise error!')
	1685	raise 'Bad file did not raise error!'
	1686
	1687	#clean up
	1688	os.remove(mesh_file)
	1689	os.remove(point_file)
	1690
	1691	def test_fit_to_mesh_file_errorsII(self):
	1692	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1693	import tempfile
	1694	import os
	1695
	1696	# create a .tsh file, no user outline
	1697	mesh_file = tempfile.mktemp(".tsh")
	1698	fd = open(mesh_file,'w')
	1699	fd.write("unit testing a bad .tsh file \n")
	1700	fd.close()
	1701
	1702	# create an .xya file
	1703	point_file = tempfile.mktemp(".xya")
	1704	fd = open(point_file,'w')
	1705	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1706	fd.close()
	1707
	1708	mesh_output_file = "new_triangle.tsh"
	1709	try:
	1710	fit_to_mesh_file(mesh_file, point_file,
	1711	mesh_output_file, display_errors = False)
	1712	except IOError:
	1713	pass
	1714	else:
	1715	raise 'Bad file did not raise error!'
	1716
	1717	#clean up
	1718	os.remove(mesh_file)
	1719	os.remove(point_file)
	1720
	1721	def test_fit_to_mesh_file_errorsIII(self):
	1722	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1723	import tempfile
	1724	import os
	1725
	1726	# create a .tsh file, no user outline
	1727	mesh_dic = {}
	1728	mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
	1729	mesh_dic['triangles'] = [[0, 2, 1]]
	1730	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1731	mesh_dic['triangle_tags'] = ['']
	1732	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1733	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1734	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1735	mesh_dic['segment_tags'] = ['external', 'external','external']
	1736	mesh_file = tempfile.mktemp(".tsh")
	1737	export_mesh_file(mesh_file,mesh_dic)
	1738
	1739	# create an .xya file
	1740	point_file = tempfile.mktemp(".xya")
	1741	fd = open(point_file,'w')
	1742	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1743	fd.close()
	1744
	1745	#This a deliberately illegal filename to invoke the error.
	1746	mesh_output_file = ".../\z\z:ya.tsh"
	1747
	1748	try:
	1749	fit_to_mesh_file(mesh_file, point_file,
	1750	mesh_output_file, display_errors = False)
	1751	except IOError:
	1752	pass
	1753	else:
	1754	raise 'Bad file did not raise error!'
	1755
	1756	#clean up
	1757	os.remove(mesh_file)
	1758	os.remove(point_file)
	1759
	1760	## FIXME? Running from the Comand line isn't in vogue these days
	1761	# The test was breaking when test_all at the inundation level was running
	1762	# was running it.issue - not running the test in this directory
	1763	def Bad_test_fit_to_mesh_file_errorsIV(self):
	1764	import os
	1765	command = '%s least_squares.py q q q e n 0.9 n' %(sys.executable)
	1766	status = os.system(command)
	1767	self.failUnless(status%255 == 1,
	1768	'command prompt least_squares.py failed. Incorect exit status.')
	1769
	1770	def test_fit_to_msh_netcdf_fileII(self):
	1771	from load_mesh.loadASCII import import_mesh_file, export_mesh_file
	1772	import tempfile
	1773	import os
	1774
	1775	# create a .tsh file, no user outline
	1776	mesh_dic = {}
	1777	mesh_dic['vertices'] = [[0.0, 0.0],
	1778	[0.0, 5.0],
	1779	[5.0, 0.0]]
	1780	mesh_dic['triangles'] = [[0, 2, 1]]
	1781	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
	1782	mesh_dic['triangle_tags'] = ['']
	1783	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1784	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1785	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
	1786	mesh_dic['segment_tags'] = ['external',
	1787	'external',
	1788	'external']
	1789	mesh_file = tempfile.mktemp(".msh")
	1790	export_mesh_file(mesh_file,mesh_dic)
	1791
	1792	# create an .xya file
	1793	point_file = tempfile.mktemp(".xya")
	1794	fd = open(point_file,'w')
	1795	fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
	1796	fd.close()
	1797
	1798	mesh_output_file = "new_triangle.msh"
	1799	fit_to_mesh_file(mesh_file,
	1800	point_file,
	1801	mesh_output_file,
	1802	alpha = 0.0)
	1803	# load in the .tsh file we just wrote
	1804	mesh_dic = import_mesh_file(mesh_output_file)
	1805
	1806	assert allclose(mesh_dic['vertex_attributes'],
	1807	[[1.0, 2.0,0.0, 0.0],
	1808	[1.0, 2.0,5.0, 10.0],
	1809	[1.0, 2.0,5.0,10.0]])
	1810
	1811	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
	1812	['density', 'temp','elevation','stage'],
	1813	'test_fit_to_mesh_file failed')
	1814
	1815	#clean up
	1816	os.remove(mesh_file)
	1817	os.remove(mesh_output_file)
	1818	os.remove(point_file)
	1819
[2585]	1820
	1821
	1822	def test_fit_using_fit_to_mesh(self):
	1823	"""Fit a surface to one set of points. Then interpolate that surface
	1824	using another set of points.
	1825	"""
	1826	from mesh import Mesh
	1827
	1828
	1829	#Setup mesh used to represent fitted function
	1830	a = [0.0, 0.0]
	1831	b = [0.0, 2.0]
	1832	c = [2.0, 0.0]
	1833	d = [0.0, 4.0]
	1834	e = [2.0, 2.0]
	1835	f = [4.0, 0.0]
	1836
	1837	points = [a, b, c, d, e, f]
	1838	#bac, bce, ecf, dbe, daf, dae
	1839	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1840
	1841	#Datapoints to fit from
	1842	data_points1 = [[ 0.66666667, 0.66666667],
	1843	[ 1.33333333, 1.33333333],
	1844	[ 2.66666667, 0.66666667],
	1845	[ 0.66666667, 2.66666667],
	1846	[ 0.0, 1.0],
	1847	[ 0.0, 3.0],
	1848	[ 1.0, 0.0],
	1849	[ 1.0, 1.0],
	1850	[ 15, -17], #Outside mesh
	1851	[ 1.0, 2.0],
	1852	[ 1.0, 3.0],
	1853	[ 2.0, 1.0],
	1854	[ 3.0, 0.0],
	1855	[ 3.0, 1.0]]
	1856
	1857	#Fit surface to mesh
	1858	z = linear_function(data_points1) #Example z-values
[2666]	1859	v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0,
	1860	precrop=True, verbose=False)
	1861
[2585]	1862	assert allclose(linear_function(points), v)
	1863
[2666]	1864
[2585]	1865
[2666]	1866	def test_acceptable_overshoot(self):
	1867	"""Fit a surface to one set of points. Then interpolate that surface
	1868	using another set of points.
	1869	Check that exceedance in fitted values are caught.
	1870	"""
	1871	from mesh import Mesh
	1872
	1873
	1874	#Setup mesh used to represent fitted function
	1875	a = [0.0, 0.0]
	1876	b = [0.0, 2.0]
	1877	c = [2.0, 0.0]
	1878	d = [0.0, 4.0]
	1879	e = [2.0, 2.0]
	1880	f = [4.0, 0.0]
	1881
	1882	points = [a, b, c, d, e, f]
	1883	#bac, bce, ecf, dbe, daf, dae
	1884	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1885
	1886	#Datapoints to fit from
	1887	data_points1 = [[ 0.66666667, 0.66666667],
	1888	[ 1.33333333, 1.33333333],
	1889	[ 2.66666667, 0.66666667],
	1890	[ 0.66666667, 2.66666667],
	1891	[ 0.0, 1.0],
	1892	[ 0.0, 3.0],
	1893	[ 1.0, 0.0],
	1894	[ 1.0, 1.0],
	1895	[ 15, -17], #Outside mesh
	1896	[ 1.0, 2.0],
	1897	[ 1.0, 3.0],
	1898	[ 2.0, 1.0],
	1899	[ 3.0, 0.0],
	1900	[ 3.0, 1.0]]
	1901
	1902	#Fit surface to mesh
	1903	z = linear_function(data_points1) #Example z-values
	1904
	1905	try:
	1906	v = fit_to_mesh(points, triangles, data_points1, z, alpha=0.0,
	1907	acceptable_overshoot = 0.2,
	1908	precrop=True, verbose=False)
	1909	except FittingError, e:
	1910	pass
	1911	else:
	1912	raise 'Should have raised exception'
	1913
	1914
	1915	#assert allclose(linear_function(points), v)
	1916
	1917
	1918
[2526]	1919	#-------------------------------------------------------------
	1920	if __name__ == "__main__":
[2665]	1921	#suite = unittest.makeSuite(Test_Least_Squares,'test_smooth_attributes_to_mesh_function')
[2778]	1922	#suite = unittest.makeSuite(Test_Least_Squares,'test_datapoint_at_centroid')
[2666]	1923	suite = unittest.makeSuite(Test_Least_Squares,'test')
[2526]	1924
	1925	#suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
	1926	#suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
	1927	runner = unittest.TextTestRunner(verbosity=1)
	1928	runner.run(suite)
	1929
	1930
	1931
	1932
	1933

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

Download in other formats: