Context Navigation

source: inundation/ga/storm_surge/pyvolution/test_least_squares.py @ 999

Last change on this file since 999 was 997, checked in by ole, 20 years ago
Implemented precropping facility in least squares Also optimised inside_polygon
File size: 36.0 KB

Rev	Line
[301]	1	#!/usr/bin/env python
	2
[494]	3	#TEST
	4
[301]	5	import unittest
	6	from math import sqrt
	7
[466]	8
[301]	9	from least_squares import *
[353]	10	from Numeric import allclose, array, transpose
[301]	11
	12	def distance(x, y):
	13	return sqrt( sum( (array(x)-array(y))**2 ))
[331]	14
	15	def linear_function(point):
	16	point = array(point)
	17	return point[:,0]+point[:,1]
	18
[301]	19
	20	class TestCase(unittest.TestCase):
[331]	21
[454]	22	def setUp(self):
	23	pass
[301]	24
[454]	25	def tearDown(self):
	26	pass
[301]	27
	28	def test_datapoint_at_centroid(self):
	29	a = [0.0, 0.0]
	30	b = [0.0, 2.0]
	31	c = [2.0,0.0]
	32	points = [a, b, c]
[309]	33	vertices = [ [1,0,2] ] #bac
[301]	34
	35	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	36
[321]	37	interp = Interpolation(points, vertices, data)
[424]	38	assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
[301]	39
	40
[608]	41	def test_quad_tree(self):
	42	p0 = [-10.0, -10.0]
	43	p1 = [20.0, -10.0]
	44	p2 = [-10.0, 20.0]
	45	p3 = [10.0, 50.0]
	46	p4 = [30.0, 30.0]
	47	p5 = [50.0, 10.0]
	48	p6 = [40.0, 60.0]
	49	p7 = [60.0, 40.0]
	50	p8 = [-66.0, 20.0]
	51	p9 = [10.0, -66.0]
	52
	53	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	54	triangles = [ [0, 1, 2],
	55	[3, 2, 4],
	56	[4, 2, 1],
	57	[4, 1, 5],
	58	[3, 4, 6],
	59	[6, 4, 7],
	60	[7, 4, 5],
	61	[8, 0, 2],
	62	[0, 9, 1]]
[301]	63
[608]	64	data = [ [4,4] ]
[883]	65	interp = Interpolation(points, triangles, data, alpha = 0.0,
	66	max_points_per_cell = 4)
[608]	67	#print "PDSG - interp.get_A()", interp.get_A()
	68	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	69	0., 0. , 0., 0., 0., 0.]]
[611]	70	assert allclose(interp.get_A(), answer)
	71	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	72	#print "PDSG - interp.get_A()", interp.get_A()
	73	answer = [ [ 0.0, 0.0, 0.0, 0.,
	74	0., 0. , 0., 0., 0., 0.]]
	75	assert allclose(interp.get_A(), answer)
	76
[608]	77
[611]	78	#point outside of quad tree root cell
	79	interp.set_point_coordinates([[-70, -70]])
	80	#print "PDSG - interp.get_A()", interp.get_A()
	81	answer = [ [ 0.0, 0.0, 0.0, 0.,
	82	0., 0. , 0., 0., 0., 0.]]
	83	assert allclose(interp.get_A(), answer)
[705]	84
[883]	85	def test_expand_search(self):
	86	p0 = [-10.0, -10.0]
	87	p1 = [20.0, -10.0]
	88	p2 = [-10.0, 20.0]
	89	p3 = [10.0, 50.0]
	90	p4 = [30.0, 30.0]
	91	p5 = [50.0, 10.0]
	92	p6 = [40.0, 60.0]
	93	p7 = [60.0, 40.0]
	94	p8 = [-66.0, 20.0]
	95	p9 = [10.0, -66.0]
	96
	97	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	98	triangles = [ [0, 1, 2],
	99	[3, 2, 4],
	100	[4, 2, 1],
	101	[4, 1, 5],
	102	[3, 4, 6],
	103	[6, 4, 7],
	104	[7, 4, 5],
	105	[8, 0, 2],
	106	[0, 9, 1]]
	107
	108	data = [ [4,4],
	109	[-30,10],
	110	[-20,0],
	111	[-20,10],
	112	[0,30],
	113	[10,-40],
	114	[10,-30],
	115	[10,-20],
	116	[10,10],
	117	[10,20],
	118	[10,30],
	119	[10,40],
	120	[20,10],
	121	[25,45],
	122	[30,0],
	123	[30,10],
	124	[30,30],
	125	[30,40],
	126	[30,50],
	127	[40,10],
	128	[40,30],
	129	[40,40],
	130	[40,50],
	131	[50,20],
	132	[50,30],
	133	[50,40],
	134	[50,50],
	135	[30,0],
	136	[-20,-20]]
	137	point_attributes = [ -400000,
	138	10,
	139	10,
	140	10,
	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	99]
	166
	167	interp = Interpolation(points, triangles, data,
	168	alpha=0.0, expand_search=False, #verbose = True, #False,
	169	max_points_per_cell = 4)
	170	calc = interp.fit_points(point_attributes, )
	171	#print "calc",calc
	172
	173	# the point at 4,4 is ignored. An expanded search has to be done
	174	# to fine which triangel it's in.
	175	# An expanded search isn't done to find that the last point
	176	# isn't in the mesh. But this isn't tested.
	177	answer= [ 10,
	178	10,
	179	10,
	180	10,
	181	10,
	182	10,
	183	10,
	184	10,
	185	10,
	186	10]
	187	assert allclose(calc, answer)
	188
[705]	189	def test_quad_treeII(self):
	190	p0 = [-66.0, 14.0]
	191	p1 = [14.0, -66.0]
	192	p2 = [14.0, 14.0]
	193	p3 = [60.0, 20.0]
	194	p4 = [10.0, 60.0]
	195	p5 = [60.0, 60.0]
[611]	196
[705]	197	points = [p0, p1, p2, p3, p4, p5]
	198	triangles = [ [0, 1, 2],
	199	[3, 2, 1],
	200	[0, 2, 4],
	201	[0, 2, 4],
	202	[4, 2, 5],
	203	[5, 2, 3]]
[608]	204
[883]	205	data = [ [-26.0,-26.0] ]
	206	interp = Interpolation(points, triangles, data, alpha = 0.0,
	207	max_points_per_cell = 4)
[705]	208	#print "PDSG - interp.get_A()", interp.get_A()
	209	answer = [ [ 0.5, 0.5, 0.0, 0.,
	210	0., 0.]]
	211	assert allclose(interp.get_A(), answer)
	212	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	213	#print "PDSG -30,-30 - interp.get_A()", interp.get_A()
	214	answer = [ [ 0.0, 0.0, 0.0, 0.,
	215	0., 0.]]
	216	assert allclose(interp.get_A(), answer)
	217
	218
	219	#point outside of quad tree root cell
	220	interp.set_point_coordinates([[-70, -70]])
	221	#print "PDSG -70,-70 interp.get_A()", interp.get_A()
	222	answer = [ [ 0.0, 0.0, 0.0, 0.,
	223	0., 0. ]]
	224	assert allclose(interp.get_A(), answer)
	225
	226
[301]	227	def test_datapoints_at_vertices(self):
[454]	228	"""Test that data points coinciding with vertices yield a diagonal matrix
	229	"""
	230
[301]	231	a = [0.0, 0.0]
	232	b = [0.0, 2.0]
	233	c = [2.0,0.0]
	234	points = [a, b, c]
[309]	235	vertices = [ [1,0,2] ] #bac
[301]	236
	237	data = points #Use data at vertices
	238
[321]	239	interp = Interpolation(points, vertices, data)
[424]	240	assert allclose(interp.get_A(), [[1., 0., 0.],
[331]	241	[0., 1., 0.],
	242	[0., 0., 1.]])
[301]	243
	244
	245
	246	def test_datapoints_on_edge_midpoints(self):
[454]	247	"""Try datapoints midway on edges -
	248	each point should affect two matrix entries equally
	249	"""
	250
[301]	251	a = [0.0, 0.0]
	252	b = [0.0, 2.0]
	253	c = [2.0,0.0]
	254	points = [a, b, c]
[309]	255	vertices = [ [1,0,2] ] #bac
[301]	256
	257	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	258
[321]	259	interp = Interpolation(points, vertices, data)
[301]	260
[424]	261	assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
[331]	262	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	263	[0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
[301]	264
	265
	266	def test_datapoints_on_edges(self):
[454]	267	"""Try datapoints on edges -
	268	each point should affect two matrix entries in proportion
	269	"""
[301]	270
	271	a = [0.0, 0.0]
	272	b = [0.0, 2.0]
	273	c = [2.0,0.0]
	274	points = [a, b, c]
[309]	275	vertices = [ [1,0,2] ] #bac
[301]	276
	277	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	278
[321]	279	interp = Interpolation(points, vertices, data)
[301]	280
[424]	281	assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
[331]	282	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	283	[0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
[301]	284
	285	def test_arbitrary_datapoints(self):
[454]	286	"""Try arbitrary datapoints
	287	"""
[301]	288
	289	from Numeric import sum
	290
	291	a = [0.0, 0.0]
	292	b = [0.0, 2.0]
	293	c = [2.0,0.0]
	294	points = [a, b, c]
[309]	295	vertices = [ [1,0,2] ] #bac
[301]	296
	297	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	298
[321]	299	interp = Interpolation(points, vertices, data)
[705]	300	#print "interp.get_A()", interp.get_A()
[424]	301	assert allclose(sum(interp.get_A(), axis=1), 1.0)
[301]	302
[997]	303	def test_arbitrary_datapoints_some_outside(self):
	304	"""Try arbitrary datapoints one outside the triangle.
	305	That one should be ignored
	306	"""
[301]	307
[997]	308	from Numeric import sum
	309
	310	a = [0.0, 0.0]
	311	b = [0.0, 2.0]
	312	c = [2.0,0.0]
	313	points = [a, b, c]
	314	vertices = [ [1,0,2] ] #bac
	315
	316	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
	317
	318
	319	interp = Interpolation(points, vertices, data, precrop = True)
	320	assert allclose(sum(interp.get_A(), axis=1), 1.0)
	321
	322	interp = Interpolation(points, vertices, data, precrop = False)
	323	assert allclose(sum(interp.get_A(), axis=1), [1,1,1,0])
	324
	325
	326
	327	# this causes a memory error in scipy.sparse
[441]	328	def test_more_triangles(self):
[454]	329
[309]	330	a = [-1.0, 0.0]
	331	b = [3.0, 4.0]
	332	c = [4.0,1.0]
	333	d = [-3.0, 2.0] #3
	334	e = [-1.0,-2.0]
	335	f = [1.0, -2.0] #5
[301]	336
[309]	337	points = [a, b, c, d,e,f]
	338	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
[301]	339
	340	#Data points
[454]	341	data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
[424]	342	interp = Interpolation(points, triangles, data_points)
[452]	343
[568]	344	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
[331]	345	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
[309]	346	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
[331]	347	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	348	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
[568]	349	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
[301]	350
[567]	351
	352	A = interp.get_A()
	353	for i in range(A.shape[0]):
	354	for j in range(A.shape[1]):
	355	if not allclose(A[i,j], answer[i][j]):
	356	print i,j,':',A[i,j], answer[i][j]
	357
	358
[454]	359	assert allclose(interp.get_A(), answer)
[331]	360
	361
	362
[424]	363
[430]	364	def test_smooth_attributes_to_mesh(self):
[309]	365	a = [0.0, 0.0]
	366	b = [0.0, 5.0]
	367	c = [5.0, 0.0]
	368	points = [a, b, c]
	369	triangles = [ [1,0,2] ] #bac
	370
	371	d1 = [1.0, 1.0]
	372	d2 = [1.0, 3.0]
	373	d3 = [3.0,1.0]
	374	z1 = 2
	375	z2 = 4
	376	z3 = 4
[331]	377	data_coords = [d1, d2, d3]
[301]	378
[436]	379	interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
[309]	380	z = [z1, z2, z3]
[331]	381	f = interp.fit(z)
[309]	382	answer = [0, 5., 5.]
[454]	383
	384	#print "f\n",f
	385	#print "answer\n",answer
[301]	386
[485]	387	assert allclose(f, answer, atol=1e-7)
[331]	388
[454]	389
[331]	390	def test_smooth_att_to_meshII(self):
[441]	391
[316]	392	a = [0.0, 0.0]
	393	b = [0.0, 5.0]
	394	c = [5.0, 0.0]
	395	points = [a, b, c]
	396	triangles = [ [1,0,2] ] #bac
	397
	398	d1 = [1.0, 1.0]
	399	d2 = [1.0, 2.0]
	400	d3 = [3.0,1.0]
[331]	401	data_coords = [d1, d2, d3]
	402	z = linear_function(data_coords)
	403	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
	404	f = interp.fit(z)
	405	answer = linear_function(points)
	406
[316]	407	assert allclose(f, answer)
	408
[328]	409	def test_smooth_attributes_to_meshIII(self):
[454]	410
[316]	411	a = [-1.0, 0.0]
	412	b = [3.0, 4.0]
	413	c = [4.0,1.0]
	414	d = [-3.0, 2.0] #3
	415	e = [-1.0,-2.0]
	416	f = [1.0, -2.0] #5
	417
	418	vertices = [a, b, c, d,e,f]
[485]	419	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
[316]	420
	421	point_coords = [[-2.0, 2.0],
	422	[-1.0, 1.0],
	423	[0.0,2.0],
	424	[1.0, 1.0],
	425	[2.0, 1.0],
	426	[0.0,0.0],
	427	[1.0, 0.0],
	428	[0.0, -1.0],
	429	[-0.2,-0.5],
	430	[-0.9, -1.5],
	431	[0.5, -1.9],
	432	[3.0,1.0]]
[452]	433
[331]	434	z = linear_function(point_coords)
	435	interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
[452]	436
	437	#print 'z',z
[331]	438	f = interp.fit(z)
	439	answer = linear_function(vertices)
[454]	440	#print "f\n",f
	441	#print "answer\n",answer
[316]	442	assert allclose(f, answer)
	443
[328]	444
	445	def test_smooth_attributes_to_meshIV(self):
[454]	446	""" Testing 2 attributes smoothed to the mesh
	447	"""
	448
[328]	449	a = [0.0, 0.0]
	450	b = [0.0, 5.0]
	451	c = [5.0, 0.0]
	452	points = [a, b, c]
	453	triangles = [ [1,0,2] ] #bac
	454
	455	d1 = [1.0, 1.0]
	456	d2 = [1.0, 3.0]
[485]	457	d3 = [3.0, 1.0]
	458	z1 = [2, 4]
	459	z2 = [4, 8]
	460	z3 = [4, 8]
	461	data_coords = [d1, d2, d3]
[328]	462
[331]	463	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
[328]	464	z = [z1, z2, z3]
[436]	465	f = interp.fit_points(z)
[328]	466	answer = [[0,0], [5., 10.], [5., 10.]]
	467	assert allclose(f, answer)
	468
	469	def test_interpolate_attributes_to_points(self):
[316]	470	v0 = [0.0, 0.0]
	471	v1 = [0.0, 5.0]
	472	v2 = [5.0, 0.0]
	473
	474	vertices = [v0, v1, v2]
	475	triangles = [ [1,0,2] ] #bac
	476
	477	d0 = [1.0, 1.0]
	478	d1 = [1.0, 2.0]
[485]	479	d2 = [3.0, 1.0]
[316]	480	point_coords = [ d0, d1, d2]
	481
[321]	482	interp = Interpolation(vertices, triangles, point_coords)
[331]	483	f = linear_function(vertices)
	484	z = interp.interpolate(f)
	485	answer = linear_function(point_coords)
	486
	487
[316]	488	assert allclose(z, answer)
[331]	489
[316]	490
[328]	491	def test_interpolate_attributes_to_pointsII(self):
[316]	492	a = [-1.0, 0.0]
	493	b = [3.0, 4.0]
[485]	494	c = [4.0, 1.0]
[316]	495	d = [-3.0, 2.0] #3
[485]	496	e = [-1.0, -2.0]
[316]	497	f = [1.0, -2.0] #5
	498
	499	vertices = [a, b, c, d,e,f]
[485]	500	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
[316]	501
	502
	503	point_coords = [[-2.0, 2.0],
	504	[-1.0, 1.0],
[485]	505	[0.0, 2.0],
[316]	506	[1.0, 1.0],
	507	[2.0, 1.0],
[485]	508	[0.0, 0.0],
[316]	509	[1.0, 0.0],
	510	[0.0, -1.0],
[485]	511	[-0.2, -0.5],
[316]	512	[-0.9, -1.5],
	513	[0.5, -1.9],
[485]	514	[3.0, 1.0]]
[316]	515
[321]	516	interp = Interpolation(vertices, triangles, point_coords)
[331]	517	f = linear_function(vertices)
	518	z = interp.interpolate(f)
	519	answer = linear_function(point_coords)
[316]	520	#print "z",z
	521	#print "answer",answer
[321]	522	assert allclose(z, answer)
[353]	523
	524	def test_interpolate_attributes_to_pointsIII(self):
[854]	525	"""Test linear interpolation of known values at vertices to
	526	new points inside a triangle
	527	"""
[858]	528	a = [0.0, 0.0]
	529	b = [0.0, 5.0]
	530	c = [5.0, 0.0]
	531	d = [5.0, 5.0]
[353]	532
[858]	533	vertices = [a, b, c, d]
	534	triangles = [ [1,0,2], [2,3,0] ] #bac, cdb
	535
	536	#Points within triangle 1
[353]	537	d0 = [1.0, 1.0]
	538	d1 = [1.0, 2.0]
[854]	539	d2 = [3.0, 1.0]
[858]	540
	541	#Point within triangle 2
	542	d3 = [4.0, 3.0]
	543
	544	#Points on common edge
	545	d4 = [2.5, 2.5]
	546	d5 = [4.0, 1.0]
	547
	548	#Point on common vertex
	549	d6 = [0., 5.]
	550
[353]	551
[858]	552	point_coords = [d0, d1, d2, d3, d4, d5, d6]
	553
[353]	554	interp = Interpolation(vertices, triangles, point_coords)
[854]	555
[858]	556	#Known values at vertices
[854]	557	#Functions are x+y, x+2y, 2x+y, x-y-5
[858]	558	f = [ [0., 0., 0., -5.], # (0,0)
	559	[5., 10., 5., -10.], # (0,5)
	560	[5., 5., 10.0, 0.], # (5,0)
	561	[10., 15., 15., -5.]] # (5,5)
	562
[353]	563	z = interp.interpolate(f)
[858]	564	answer = [ [2., 3., 3., -5.], # (1,1)
	565	[3., 5., 4., -6.], # (1,2)
	566	[4., 5., 7., -3.], # (3,1)
	567	[7., 10., 11., -4.], # (4,3)
	568	[5., 7.5, 7.5, -5.], # (2.5, 2.5)
	569	[5., 6., 9., -2.], # (4,1)
	570	[5., 10., 5., -10.]] # (0,5)
[321]	571
[353]	572	#print "***********"
	573	#print "z",z
	574	#print "answer",answer
	575	#print "***********"
	576
	577	assert allclose(z, answer)
	578
	579	def test_interpolate_attributes_to_pointsIV(self):
	580	a = [-1.0, 0.0]
	581	b = [3.0, 4.0]
[485]	582	c = [4.0, 1.0]
[353]	583	d = [-3.0, 2.0] #3
[485]	584	e = [-1.0, -2.0]
[353]	585	f = [1.0, -2.0] #5
	586
	587	vertices = [a, b, c, d,e,f]
[485]	588	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
[353]	589
	590
	591	point_coords = [[-2.0, 2.0],
	592	[-1.0, 1.0],
[485]	593	[0.0, 2.0],
[353]	594	[1.0, 1.0],
	595	[2.0, 1.0],
[485]	596	[0.0, 0.0],
[353]	597	[1.0, 0.0],
	598	[0.0, -1.0],
[485]	599	[-0.2, -0.5],
[353]	600	[-0.9, -1.5],
	601	[0.5, -1.9],
[485]	602	[3.0, 1.0]]
[321]	603
[353]	604	interp = Interpolation(vertices, triangles, point_coords)
	605	f = array([linear_function(vertices),2*linear_function(vertices) ])
	606	f = transpose(f)
	607	#print "f",f
	608	z = interp.interpolate(f)
[854]	609	answer = [linear_function(point_coords),
	610	2*linear_function(point_coords) ]
[353]	611	answer = transpose(answer)
	612	#print "z",z
	613	#print "answer",answer
	614	assert allclose(z, answer)
[328]	615
	616	def test_smooth_attributes_to_mesh_function(self):
[454]	617	""" Testing 2 attributes smoothed to the mesh
	618	"""
[331]	619
[321]	620	a = [0.0, 0.0]
	621	b = [0.0, 5.0]
	622	c = [5.0, 0.0]
	623	points = [a, b, c]
	624	triangles = [ [1,0,2] ] #bac
	625
	626	d1 = [1.0, 1.0]
	627	d2 = [1.0, 3.0]
[485]	628	d3 = [3.0, 1.0]
	629	z1 = [2, 4]
	630	z2 = [4, 8]
	631	z3 = [4, 8]
	632	data_coords = [d1, d2, d3]
[333]	633	z = [z1, z2, z3]
[321]	634
[333]	635	f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
[331]	636	answer = [[0, 0], [5., 10.], [5., 10.]]
	637
[321]	638	assert allclose(f, answer)
[331]	639
	640
[820]	641
[836]	642	def test_pts2rectangular(self):
[331]	643
[820]	644	import time, os
	645	FN = 'xyatest' + str(time.time()) + '.xya'
	646	fid = open(FN, 'w')
[841]	647	fid.write(' %s \n' %('elevation'))
	648	fid.write('%f %f %f\n' %(1,1,2) )
	649	fid.write('%f %f %f\n' %(1,3,4) )
	650	fid.write('%f %f %f\n' %(3,1,4) )
[820]	651	fid.close()
	652
[834]	653	points, triangles, boundary, attributes =\
[836]	654	pts2rectangular(FN, 4, 4, format = 'asc')
[820]	655
[841]	656
[820]	657	data_coords = [ [1,1], [1,3], [3,1] ]
[841]	658	z = [2, 4, 4]
[820]	659
	660	ref = fit_to_mesh(points, triangles, data_coords, z)
[841]	661
	662	#print attributes
	663	#print ref
[820]	664	assert allclose(attributes, ref)
	665
[835]	666	os.remove(FN)
[820]	667
	668
[331]	669	#Tests of smoothing matrix
	670	def test_smoothing_matrix_one_triangle(self):
	671	from Numeric import dot
	672	a = [0.0, 0.0]
	673	b = [0.0, 2.0]
	674	c = [2.0,0.0]
	675	points = [a, b, c]
[321]	676
[331]	677	vertices = [ [1,0,2] ] #bac
	678
	679	interp = Interpolation(points, vertices)
	680
[424]	681	assert allclose(interp.get_D(), [[1, -0.5, -0.5],
[331]	682	[-0.5, 0.5, 0],
	683	[-0.5, 0, 0.5]])
	684
	685	#Define f(x,y) = x
	686	f = array([0,0,2]) #Value at global vertex 2
	687
	688	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	689	# int 1 dx dy = area = 2
[424]	690	assert dot(dot(f, interp.get_D()), f) == 2
[331]	691
	692	#Define f(x,y) = y
	693	f = array([0,2,0]) #Value at global vertex 1
	694
	695	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	696	# int 1 dx dy = area = 2
[424]	697	assert dot(dot(f, interp.get_D()), f) == 2
[331]	698
	699	#Define f(x,y) = x+y
	700	f = array([0,2,2]) #Values at global vertex 1 and 2
	701
	702	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	703	# int 2 dx dy = 2*area = 4
[424]	704	assert dot(dot(f, interp.get_D()), f) == 4
[331]	705
	706
	707
	708	def test_smoothing_matrix_more_triangles(self):
	709	from Numeric import dot
	710
	711	a = [0.0, 0.0]
	712	b = [0.0, 2.0]
	713	c = [2.0,0.0]
	714	d = [0.0, 4.0]
	715	e = [2.0, 2.0]
	716	f = [4.0,0.0]
	717
	718	points = [a, b, c, d, e, f]
	719	#bac, bce, ecf, dbe, daf, dae
	720	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	721
	722	interp = Interpolation(points, vertices)
	723
	724
[424]	725	#assert allclose(interp.get_D(), [[1, -0.5, -0.5],
[331]	726	# [-0.5, 0.5, 0],
	727	# [-0.5, 0, 0.5]])
	728
	729	#Define f(x,y) = x
	730	f = array([0,0,2,0,2,4]) #f evaluated at points a-f
	731
	732	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	733	# int 1 dx dy = total area = 8
[424]	734	assert dot(dot(f, interp.get_D()), f) == 8
[331]	735
	736	#Define f(x,y) = y
	737	f = array([0,2,0,4,2,0]) #f evaluated at points a-f
	738
	739	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	740	# int 1 dx dy = area = 8
[424]	741	assert dot(dot(f, interp.get_D()), f) == 8
[331]	742
	743	#Define f(x,y) = x+y
	744	f = array([0,2,2,4,4,4]) #f evaluated at points a-f
	745
	746	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	747	# int 2 dx dy = 2*area = 16
[424]	748	assert dot(dot(f, interp.get_D()), f) == 16
[331]	749
	750
	751	def test_fit_and_interpolation(self):
	752	from mesh import Mesh
	753
	754	a = [0.0, 0.0]
	755	b = [0.0, 2.0]
	756	c = [2.0, 0.0]
	757	d = [0.0, 4.0]
	758	e = [2.0, 2.0]
	759	f = [4.0, 0.0]
	760
	761	points = [a, b, c, d, e, f]
	762	#bac, bce, ecf, dbe, daf, dae
	763	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	764
	765	#Get (enough) datapoints
	766	data_points = [[ 0.66666667, 0.66666667],
	767	[ 1.33333333, 1.33333333],
	768	[ 2.66666667, 0.66666667],
	769	[ 0.66666667, 2.66666667],
	770	[ 0.0, 1.0],
	771	[ 0.0, 3.0],
	772	[ 1.0, 0.0],
	773	[ 1.0, 1.0],
	774	[ 1.0, 2.0],
	775	[ 1.0, 3.0],
	776	[ 2.0, 1.0],
	777	[ 3.0, 0.0],
	778	[ 3.0, 1.0]]
	779
	780	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	781
	782	z = linear_function(data_points)
	783	answer = linear_function(points)
	784
	785	f = interp.fit(z)
	786
[454]	787	#print "f",f
	788	#print "answer",answer
[331]	789	assert allclose(f, answer)
	790
	791	#Map back
	792	z1 = interp.interpolate(f)
[454]	793	#print "z1\n", z1
	794	#print "z\n",z
[331]	795	assert allclose(z, z1)
	796
	797
	798	def test_smoothing_and_interpolation(self):
	799
	800	a = [0.0, 0.0]
	801	b = [0.0, 2.0]
	802	c = [2.0, 0.0]
	803	d = [0.0, 4.0]
	804	e = [2.0, 2.0]
	805	f = [4.0, 0.0]
	806
	807	points = [a, b, c, d, e, f]
	808	#bac, bce, ecf, dbe, daf, dae
	809	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	810
	811	#Get (too few!) datapoints
	812	data_points = [[ 0.66666667, 0.66666667],
	813	[ 1.33333333, 1.33333333],
	814	[ 2.66666667, 0.66666667],
	815	[ 0.66666667, 2.66666667]]
	816
	817	z = linear_function(data_points)
	818	answer = linear_function(points)
	819
	820	#Make interpolator with too few data points and no smoothing
	821	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	822	#Must raise an exception
	823	try:
	824	f = interp.fit(z)
	825	except:
	826	pass
	827
	828	#Now try with smoothing parameter
	829	interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
	830
	831	f = interp.fit(z)
	832	#f will be different from answerr due to smoothing
[454]	833	assert allclose(f, answer,atol=5)
[331]	834
	835	#Map back
	836	z1 = interp.interpolate(f)
	837	assert allclose(z, z1)
[346]	838
[642]	839
	840
	841	def test_fit_and_interpolation_with_new_points(self):
	842	"""Fit a surface to one set of points. Then interpolate that surface
	843	using another set of points.
	844	"""
	845	from mesh import Mesh
	846
	847
	848	#Setup mesh used to represent fitted function
	849	a = [0.0, 0.0]
	850	b = [0.0, 2.0]
	851	c = [2.0, 0.0]
	852	d = [0.0, 4.0]
	853	e = [2.0, 2.0]
	854	f = [4.0, 0.0]
	855
	856	points = [a, b, c, d, e, f]
	857	#bac, bce, ecf, dbe, daf, dae
	858	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	859
	860	#Datapoints to fit from
	861	data_points1 = [[ 0.66666667, 0.66666667],
	862	[ 1.33333333, 1.33333333],
	863	[ 2.66666667, 0.66666667],
	864	[ 0.66666667, 2.66666667],
	865	[ 0.0, 1.0],
	866	[ 0.0, 3.0],
	867	[ 1.0, 0.0],
	868	[ 1.0, 1.0],
[997]	869	[ 15, -17], #Outside mesh
[642]	870	[ 1.0, 2.0],
	871	[ 1.0, 3.0],
	872	[ 2.0, 1.0],
	873	[ 3.0, 0.0],
	874	[ 3.0, 1.0]]
	875
	876	#Fit surface to mesh
[997]	877	interp = Interpolation(points, triangles, data_points1, alpha=0.0,
	878	precrop = True)
[642]	879	z = linear_function(data_points1) #Example z-values
	880	f = interp.fit(z) #Fitted values at vertices
	881
	882
	883
	884	#New datapoints where interpolated values are sought
	885	data_points2 = [[ 0.0, 0.0],
	886	[ 0.5, 0.5],
	887	[ 0.7, 0.7],
[997]	888	[-13, 65], #Outside
[642]	889	[ 1.0, 0.5],
	890	[ 2.0, 0.4],
	891	[ 2.8, 1.2]]
	892
	893
	894
[997]	895	#Build new A matrix based on new points (without precrop)
	896	interp.build_interpolation_matrix_A(data_points2, precrop = False)
	897
[642]	898	#Interpolate using fitted surface
	899	z1 = interp.interpolate(f)
	900
[997]	901	#import Numeric
	902	#data_points2 = Numeric.take(data_points2, interp.point_indices)
	903
	904	#Desired result (OK for points inside)
	905
	906	answer = linear_function(data_points2)
	907	import Numeric
	908	z1 = Numeric.take(z1, [0,1,2,4,5,6])
	909	answer = Numeric.take(answer, [0,1,2,4,5,6])
	910	assert allclose(z1, answer)
	911
	912	#Build new A matrix based on new points (with precrop)
	913	interp.build_interpolation_matrix_A(data_points2, precrop = True)
	914
	915	#Interpolate using fitted surface
	916	z1 = interp.interpolate(f)
	917
	918	import Numeric
	919	data_points2 = Numeric.take(data_points2, interp.point_indices)
	920
[642]	921	#Desired result
	922	answer = linear_function(data_points2)
	923	assert allclose(z1, answer)
	924
	925
	926
[976]	927	def test_fit_and_interpolation_with_different_origins(self):
	928	"""Fit a surface to one set of points. Then interpolate that surface
	929	using another set of points.
	930	This test tests situtaion where points and mesh belong to a different
	931	coordinate system as defined by origin.
	932	"""
	933	from mesh import Mesh
[642]	934
[976]	935	#Setup mesh used to represent fitted function
	936	a = [0.0, 0.0]
	937	b = [0.0, 2.0]
	938	c = [2.0, 0.0]
	939	d = [0.0, 4.0]
	940	e = [2.0, 2.0]
	941	f = [4.0, 0.0]
	942
	943	points = [a, b, c, d, e, f]
	944	#bac, bce, ecf, dbe, daf, dae
	945	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	946
	947	#Datapoints to fit from
	948	data_points1 = [[ 0.66666667, 0.66666667],
	949	[ 1.33333333, 1.33333333],
	950	[ 2.66666667, 0.66666667],
	951	[ 0.66666667, 2.66666667],
	952	[ 0.0, 1.0],
	953	[ 0.0, 3.0],
	954	[ 1.0, 0.0],
	955	[ 1.0, 1.0],
	956	[ 1.0, 2.0],
	957	[ 1.0, 3.0],
	958	[ 2.0, 1.0],
	959	[ 3.0, 0.0],
	960	[ 3.0, 1.0]]
	961
	962
	963	#First check that things are OK when using same origin
	964	mesh_origin = (56, 290000, 618000) #zone, easting, northing
	965	data_origin = (56, 290000, 618000) #zone, easting, northing
	966
	967
	968	#Fit surface to mesh
	969	interp = Interpolation(points, triangles, data_points1,
	970	alpha=0.0,
	971	data_origin = data_origin,
	972	mesh_origin = mesh_origin)
	973
	974	z = linear_function(data_points1) #Example z-values
	975	f = interp.fit(z) #Fitted values at vertices
	976
	977
	978	#New datapoints where interpolated values are sought
	979	data_points2 = [[ 0.0, 0.0],
	980	[ 0.5, 0.5],
	981	[ 0.7, 0.7],
	982	[ 1.0, 0.5],
	983	[ 2.0, 0.4],
	984	[ 2.8, 1.2]]
	985
	986
	987	#Build new A matrix based on new points
	988	interp.build_interpolation_matrix_A(data_points2)
	989
	990	#Interpolate using fitted surface
	991	z1 = interp.interpolate(f)
	992
	993	#Desired result
	994	answer = linear_function(data_points2)
	995	assert allclose(z1, answer)
	996
	997
	998	##############################################
	999
	1000	#Then check situtaion where points are relative to a different
	1001	#origin (same zone, though, until we figure that out (FIXME))
	1002
	1003	mesh_origin = (56, 290000, 618000) #zone, easting, northing
	1004	data_origin = (56, 10000, 10000) #zone, easting, northing
	1005
	1006	#Shift datapoints according to new origin
	1007
	1008	for k in range(len(data_points1)):
	1009	data_points1[k][0] += mesh_origin[1] - data_origin[1]
	1010	data_points1[k][1] += mesh_origin[2] - data_origin[2]
	1011
	1012	for k in range(len(data_points2)):
	1013	data_points2[k][0] += mesh_origin[1] - data_origin[1]
	1014	data_points2[k][1] += mesh_origin[2] - data_origin[2]
	1015
	1016
	1017
	1018	#Fit surface to mesh
	1019	interp = Interpolation(points, triangles, data_points1,
	1020	alpha=0.0,
	1021	data_origin = data_origin,
	1022	mesh_origin = mesh_origin)
	1023
	1024	f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
	1025
	1026	assert allclose(f,f1), 'Fit should have been unaltered'
	1027
	1028
	1029	#Build new A matrix based on new points
	1030	interp.build_interpolation_matrix_A(data_points2)
	1031
	1032	#Interpolate using fitted surface
	1033	z1 = interp.interpolate(f)
	1034	assert allclose(z1, answer)
	1035
	1036
[979]	1037	#########################################################
	1038	#Finally try to relate data_points2 to new origin without
	1039	#rebuilding matrix
[976]	1040
[979]	1041	data_origin = (56, 2000, 2000) #zone, easting, northing
	1042	for k in range(len(data_points2)):
	1043	data_points2[k][0] += 8000
	1044	data_points2[k][1] += 8000
	1045
	1046	#Build new A matrix based on new points
	1047	interp.build_interpolation_matrix_A(data_points2,
	1048	data_origin = data_origin)
	1049
	1050	#Interpolate using fitted surface
	1051	z1 = interp.interpolate(f)
	1052	assert allclose(z1, answer)
	1053
	1054
	1055
	1056
[346]	1057	def test_fit_to_mesh_file(self):
	1058	from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
[889]	1059	export_triangulation_file
[346]	1060	import tempfile
	1061	import os
[331]	1062
[346]	1063	# create a .tsh file, no user outline
	1064	mesh_dic = {}
[968]	1065	mesh_dic['vertices'] = [[0.0, 0.0],
[346]	1066	[0.0, 5.0],
	1067	[5.0, 0.0]]
[968]	1068	mesh_dic['triangles'] = [[0, 2, 1]]
	1069	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
[969]	1070	mesh_dic['triangle_tags'] = [['']]
[968]	1071	mesh_dic['vertex_attributes'] = [[], [], []]
	1072	mesh_dic['vertiex_attribute_titles'] = []
	1073	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
[969]	1074	mesh_dic['segment_tags'] = ['external',
[429]	1075	'external',
	1076	'external']
	1077	mesh_file = tempfile.mktemp(".tsh")
[889]	1078	export_triangulation_file(mesh_file,mesh_dic)
[429]	1079
	1080	# create an .xya file
	1081	point_file = tempfile.mktemp(".xya")
	1082	fd = open(point_file,'w')
	1083	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")
	1084	fd.close()
	1085
	1086	mesh_output_file = "new_trianlge.tsh"
	1087	fit_to_mesh_file(mesh_file,
	1088	point_file,
	1089	mesh_output_file,
	1090	alpha = 0.0)
	1091	# load in the .tsh file we just wrote
	1092	mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
[968]	1093	#print "mesh_dic",mesh_dic
	1094	ans =[[0.0, 0.0],
	1095	[5.0, 10.0],
	1096	[5.0,10.0]]
	1097	assert allclose(mesh_dic['vertex_attributes'],ans)
[429]	1098
[968]	1099	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
[429]	1100	['elevation','stage'],
	1101	'test_fit_to_mesh_file failed')
	1102
	1103	#clean up
	1104	os.remove(mesh_file)
	1105	os.remove(point_file)
[994]	1106	os.remove(mesh_output_file)
[429]	1107
	1108	def test_fit_to_mesh_fileII(self):
	1109	from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
[889]	1110	export_triangulation_file
[429]	1111	import tempfile
	1112	import os
	1113
	1114	# create a .tsh file, no user outline
	1115	mesh_dic = {}
[968]	1116	mesh_dic['vertices'] = [[0.0, 0.0],
[429]	1117	[0.0, 5.0],
	1118	[5.0, 0.0]]
[968]	1119	mesh_dic['triangles'] = [[0, 2, 1]]
	1120	mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
[969]	1121	mesh_dic['triangle_tags'] = [['']]
[968]	1122	mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
	1123	mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
	1124	mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
[969]	1125	mesh_dic['segment_tags'] = ['external',
[346]	1126	'external',
	1127	'external']
	1128	mesh_file = tempfile.mktemp(".tsh")
[889]	1129	export_triangulation_file(mesh_file,mesh_dic)
[346]	1130
	1131	# create an .xya file
	1132	point_file = tempfile.mktemp(".xya")
	1133	fd = open(point_file,'w')
[393]	1134	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")
[346]	1135	fd.close()
	1136
[484]	1137	mesh_output_file = "new_triangle.tsh"
[346]	1138	fit_to_mesh_file(mesh_file,
	1139	point_file,
	1140	mesh_output_file,
	1141	alpha = 0.0)
[352]	1142	# load in the .tsh file we just wrote
[346]	1143	mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
	1144
[968]	1145	assert allclose(mesh_dic['vertex_attributes'],
[393]	1146	[[1.0, 2.0,0.0, 0.0],
	1147	[1.0, 2.0,5.0, 10.0],
	1148	[1.0, 2.0,5.0,10.0]])
[346]	1149
[968]	1150	self.failUnless(mesh_dic['vertex_attribute_titles'] ==
[393]	1151	['density', 'temp','elevation','stage'],
[374]	1152	'test_fit_to_mesh_file failed')
	1153
[346]	1154	#clean up
	1155	os.remove(mesh_file)
[484]	1156	os.remove(mesh_output_file)
[346]	1157	os.remove(point_file)
	1158
[301]	1159	#-------------------------------------------------------------
	1160	if __name__ == "__main__":
[997]	1161	suite = unittest.makeSuite(TestCase,'test')
[642]	1162
[705]	1163	#suite = unittest.makeSuite(TestCase,'test_arbitrary_datapoints')
[475]	1164	runner = unittest.TextTestRunner(verbosity=1)
[301]	1165	runner.run(suite)
	1166
	1167
	1168
	1169
	1170
	1171

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

Download in other formats: