Context Navigation

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

Last change on this file since 761 was 705, checked in by duncan, 20 years ago
fixed bug in least_squares algorithm.
File size: 25.8 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
	29	def test_datapoint_at_centroid(self):
	30	a = [0.0, 0.0]
	31	b = [0.0, 2.0]
	32	c = [2.0,0.0]
	33	points = [a, b, c]
[309]	34	vertices = [ [1,0,2] ] #bac
[301]	35
	36	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	37
[321]	38	interp = Interpolation(points, vertices, data)
[424]	39	assert allclose(interp.get_A(), [[1./3, 1./3, 1./3]])
[301]	40
	41
[608]	42	def test_quad_tree(self):
	43	p0 = [-10.0, -10.0]
	44	p1 = [20.0, -10.0]
	45	p2 = [-10.0, 20.0]
	46	p3 = [10.0, 50.0]
	47	p4 = [30.0, 30.0]
	48	p5 = [50.0, 10.0]
	49	p6 = [40.0, 60.0]
	50	p7 = [60.0, 40.0]
	51	p8 = [-66.0, 20.0]
	52	p9 = [10.0, -66.0]
	53
	54	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	55	triangles = [ [0, 1, 2],
	56	[3, 2, 4],
	57	[4, 2, 1],
	58	[4, 1, 5],
	59	[3, 4, 6],
	60	[6, 4, 7],
	61	[7, 4, 5],
	62	[8, 0, 2],
	63	[0, 9, 1]]
[301]	64
[608]	65	data = [ [4,4] ]
	66
	67	interp = Interpolation(points, triangles, data, alpha = 0.0)
	68	#print "PDSG - interp.get_A()", interp.get_A()
	69	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	70	0., 0. , 0., 0., 0., 0.]]
[611]	71	assert allclose(interp.get_A(), answer)
	72	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	73	#print "PDSG - interp.get_A()", interp.get_A()
	74	answer = [ [ 0.0, 0.0, 0.0, 0.,
	75	0., 0. , 0., 0., 0., 0.]]
	76	assert allclose(interp.get_A(), answer)
	77
[608]	78
[611]	79	#point outside of quad tree root cell
	80	interp.set_point_coordinates([[-70, -70]])
	81	#print "PDSG - interp.get_A()", interp.get_A()
	82	answer = [ [ 0.0, 0.0, 0.0, 0.,
	83	0., 0. , 0., 0., 0., 0.]]
	84	assert allclose(interp.get_A(), answer)
[705]	85
	86	def test_quad_treeII(self):
	87	p0 = [-66.0, 14.0]
	88	p1 = [14.0, -66.0]
	89	p2 = [14.0, 14.0]
	90	p3 = [60.0, 20.0]
	91	p4 = [10.0, 60.0]
	92	p5 = [60.0, 60.0]
[611]	93
[705]	94	points = [p0, p1, p2, p3, p4, p5]
	95	triangles = [ [0, 1, 2],
	96	[3, 2, 1],
	97	[0, 2, 4],
	98	[0, 2, 4],
	99	[4, 2, 5],
	100	[5, 2, 3]]
[608]	101
[705]	102	data = [ [-26.0,-26.0] ]
	103
	104	interp = Interpolation(points, triangles, data, alpha = 0.0)
	105	#print "PDSG - interp.get_A()", interp.get_A()
	106	answer = [ [ 0.5, 0.5, 0.0, 0.,
	107	0., 0.]]
	108	assert allclose(interp.get_A(), answer)
	109	interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	110	#print "PDSG -30,-30 - interp.get_A()", interp.get_A()
	111	answer = [ [ 0.0, 0.0, 0.0, 0.,
	112	0., 0.]]
	113	assert allclose(interp.get_A(), answer)
	114
	115
	116	#point outside of quad tree root cell
	117	interp.set_point_coordinates([[-70, -70]])
	118	#print "PDSG -70,-70 interp.get_A()", interp.get_A()
	119	answer = [ [ 0.0, 0.0, 0.0, 0.,
	120	0., 0. ]]
	121	assert allclose(interp.get_A(), answer)
	122
	123
[301]	124	def test_datapoints_at_vertices(self):
[454]	125	"""Test that data points coinciding with vertices yield a diagonal matrix
	126	"""
	127
[301]	128	a = [0.0, 0.0]
	129	b = [0.0, 2.0]
	130	c = [2.0,0.0]
	131	points = [a, b, c]
[309]	132	vertices = [ [1,0,2] ] #bac
[301]	133
	134	data = points #Use data at vertices
	135
[321]	136	interp = Interpolation(points, vertices, data)
[424]	137	assert allclose(interp.get_A(), [[1., 0., 0.],
[331]	138	[0., 1., 0.],
	139	[0., 0., 1.]])
[301]	140
	141
	142
	143	def test_datapoints_on_edge_midpoints(self):
[454]	144	"""Try datapoints midway on edges -
	145	each point should affect two matrix entries equally
	146	"""
	147
[301]	148	a = [0.0, 0.0]
	149	b = [0.0, 2.0]
	150	c = [2.0,0.0]
	151	points = [a, b, c]
[309]	152	vertices = [ [1,0,2] ] #bac
[301]	153
	154	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	155
[321]	156	interp = Interpolation(points, vertices, data)
[301]	157
[424]	158	assert allclose(interp.get_A(), [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
[331]	159	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	160	[0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
[301]	161
	162
	163	def test_datapoints_on_edges(self):
[454]	164	"""Try datapoints on edges -
	165	each point should affect two matrix entries in proportion
	166	"""
[301]	167
	168	a = [0.0, 0.0]
	169	b = [0.0, 2.0]
	170	c = [2.0,0.0]
	171	points = [a, b, c]
[309]	172	vertices = [ [1,0,2] ] #bac
[301]	173
	174	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	175
[321]	176	interp = Interpolation(points, vertices, data)
[301]	177
[424]	178	assert allclose(interp.get_A(), [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
[331]	179	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	180	[0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
[301]	181
	182	def test_arbitrary_datapoints(self):
[454]	183	"""Try arbitrary datapoints
	184	"""
[301]	185
	186	from Numeric import sum
	187
	188	a = [0.0, 0.0]
	189	b = [0.0, 2.0]
	190	c = [2.0,0.0]
	191	points = [a, b, c]
[309]	192	vertices = [ [1,0,2] ] #bac
[301]	193
	194	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	195
[321]	196	interp = Interpolation(points, vertices, data)
[705]	197	#print "interp.get_A()", interp.get_A()
[424]	198	assert allclose(sum(interp.get_A(), axis=1), 1.0)
[301]	199
	200
[454]	201	# this causes a memory error in scipy.sparce
[441]	202	def test_more_triangles(self):
[454]	203
[309]	204	a = [-1.0, 0.0]
	205	b = [3.0, 4.0]
	206	c = [4.0,1.0]
	207	d = [-3.0, 2.0] #3
	208	e = [-1.0,-2.0]
	209	f = [1.0, -2.0] #5
[301]	210
[309]	211	points = [a, b, c, d,e,f]
	212	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
[301]	213
	214	#Data points
[454]	215	data_points = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
[424]	216	interp = Interpolation(points, triangles, data_points)
[452]	217
[568]	218	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
[331]	219	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
[309]	220	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
[331]	221	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	222	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
[568]	223	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
[301]	224
[567]	225
	226	A = interp.get_A()
	227	for i in range(A.shape[0]):
	228	for j in range(A.shape[1]):
	229	if not allclose(A[i,j], answer[i][j]):
	230	print i,j,':',A[i,j], answer[i][j]
	231
	232
[454]	233	assert allclose(interp.get_A(), answer)
[331]	234
	235
	236
[424]	237
[430]	238	def test_smooth_attributes_to_mesh(self):
[309]	239	a = [0.0, 0.0]
	240	b = [0.0, 5.0]
	241	c = [5.0, 0.0]
	242	points = [a, b, c]
	243	triangles = [ [1,0,2] ] #bac
	244
	245	d1 = [1.0, 1.0]
	246	d2 = [1.0, 3.0]
	247	d3 = [3.0,1.0]
	248	z1 = 2
	249	z2 = 4
	250	z3 = 4
[331]	251	data_coords = [d1, d2, d3]
[301]	252
[436]	253	interp = Interpolation(points, triangles, data_coords, alpha=5.0e-20)
[309]	254	z = [z1, z2, z3]
[331]	255	f = interp.fit(z)
[309]	256	answer = [0, 5., 5.]
[454]	257
	258	#print "f\n",f
	259	#print "answer\n",answer
[301]	260
[485]	261	assert allclose(f, answer, atol=1e-7)
[331]	262
[454]	263
[331]	264	def test_smooth_att_to_meshII(self):
[441]	265
[316]	266	a = [0.0, 0.0]
	267	b = [0.0, 5.0]
	268	c = [5.0, 0.0]
	269	points = [a, b, c]
	270	triangles = [ [1,0,2] ] #bac
	271
	272	d1 = [1.0, 1.0]
	273	d2 = [1.0, 2.0]
	274	d3 = [3.0,1.0]
[331]	275	data_coords = [d1, d2, d3]
	276	z = linear_function(data_coords)
	277	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
	278	f = interp.fit(z)
	279	answer = linear_function(points)
	280
[316]	281	assert allclose(f, answer)
	282
[328]	283	def test_smooth_attributes_to_meshIII(self):
[454]	284
[316]	285	a = [-1.0, 0.0]
	286	b = [3.0, 4.0]
	287	c = [4.0,1.0]
	288	d = [-3.0, 2.0] #3
	289	e = [-1.0,-2.0]
	290	f = [1.0, -2.0] #5
	291
	292	vertices = [a, b, c, d,e,f]
[485]	293	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
[316]	294
	295	point_coords = [[-2.0, 2.0],
	296	[-1.0, 1.0],
	297	[0.0,2.0],
	298	[1.0, 1.0],
	299	[2.0, 1.0],
	300	[0.0,0.0],
	301	[1.0, 0.0],
	302	[0.0, -1.0],
	303	[-0.2,-0.5],
	304	[-0.9, -1.5],
	305	[0.5, -1.9],
	306	[3.0,1.0]]
[452]	307
[331]	308	z = linear_function(point_coords)
	309	interp = Interpolation(vertices, triangles, point_coords, alpha=0.0)
[452]	310
	311	#print 'z',z
[331]	312	f = interp.fit(z)
	313	answer = linear_function(vertices)
[454]	314	#print "f\n",f
	315	#print "answer\n",answer
[316]	316	assert allclose(f, answer)
	317
[328]	318
	319	def test_smooth_attributes_to_meshIV(self):
[454]	320	""" Testing 2 attributes smoothed to the mesh
	321	"""
	322
[328]	323	a = [0.0, 0.0]
	324	b = [0.0, 5.0]
	325	c = [5.0, 0.0]
	326	points = [a, b, c]
	327	triangles = [ [1,0,2] ] #bac
	328
	329	d1 = [1.0, 1.0]
	330	d2 = [1.0, 3.0]
[485]	331	d3 = [3.0, 1.0]
	332	z1 = [2, 4]
	333	z2 = [4, 8]
	334	z3 = [4, 8]
	335	data_coords = [d1, d2, d3]
[328]	336
[331]	337	interp = Interpolation(points, triangles, data_coords, alpha=0.0)
[328]	338	z = [z1, z2, z3]
[436]	339	f = interp.fit_points(z)
[328]	340	answer = [[0,0], [5., 10.], [5., 10.]]
	341	assert allclose(f, answer)
	342
	343	def test_interpolate_attributes_to_points(self):
[316]	344	v0 = [0.0, 0.0]
	345	v1 = [0.0, 5.0]
	346	v2 = [5.0, 0.0]
	347
	348	vertices = [v0, v1, v2]
	349	triangles = [ [1,0,2] ] #bac
	350
	351	d0 = [1.0, 1.0]
	352	d1 = [1.0, 2.0]
[485]	353	d2 = [3.0, 1.0]
[316]	354	point_coords = [ d0, d1, d2]
	355
[321]	356	interp = Interpolation(vertices, triangles, point_coords)
[331]	357	f = linear_function(vertices)
	358	z = interp.interpolate(f)
	359	answer = linear_function(point_coords)
	360
	361
[316]	362	assert allclose(z, answer)
[331]	363
[316]	364
[328]	365	def test_interpolate_attributes_to_pointsII(self):
[316]	366	a = [-1.0, 0.0]
	367	b = [3.0, 4.0]
[485]	368	c = [4.0, 1.0]
[316]	369	d = [-3.0, 2.0] #3
[485]	370	e = [-1.0, -2.0]
[316]	371	f = [1.0, -2.0] #5
	372
	373	vertices = [a, b, c, d,e,f]
[485]	374	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
[316]	375
	376
	377	point_coords = [[-2.0, 2.0],
	378	[-1.0, 1.0],
[485]	379	[0.0, 2.0],
[316]	380	[1.0, 1.0],
	381	[2.0, 1.0],
[485]	382	[0.0, 0.0],
[316]	383	[1.0, 0.0],
	384	[0.0, -1.0],
[485]	385	[-0.2, -0.5],
[316]	386	[-0.9, -1.5],
	387	[0.5, -1.9],
[485]	388	[3.0, 1.0]]
[316]	389
[321]	390	interp = Interpolation(vertices, triangles, point_coords)
[331]	391	f = linear_function(vertices)
	392	z = interp.interpolate(f)
	393	answer = linear_function(point_coords)
[316]	394	#print "z",z
	395	#print "answer",answer
[321]	396	assert allclose(z, answer)
[353]	397
	398	def test_interpolate_attributes_to_pointsIII(self):
	399	v0 = [0.0, 0.0]
	400	v1 = [0.0, 5.0]
	401	v2 = [5.0, 0.0]
	402
	403	vertices = [v0, v1, v2]
	404	triangles = [ [1,0,2] ] #bac
	405
	406	d0 = [1.0, 1.0]
	407	d1 = [1.0, 2.0]
	408	d2 = [3.0,1.0]
	409	point_coords = [ d0, d1, d2]
	410
	411	interp = Interpolation(vertices, triangles, point_coords)
	412	f = [ [0., 0.], [5., 10.], [5., 10.]] #linear_function(vertices)
	413	#print "f",f
	414
	415	z = interp.interpolate(f)
	416	answer = [ [2., 4.], [3., 6.], [4., 8.]]
[321]	417
[353]	418	#print "***********"
	419	#print "z",z
	420	#print "answer",answer
	421	#print "***********"
	422
	423	assert allclose(z, answer)
	424
	425	def test_interpolate_attributes_to_pointsIV(self):
	426	a = [-1.0, 0.0]
	427	b = [3.0, 4.0]
[485]	428	c = [4.0, 1.0]
[353]	429	d = [-3.0, 2.0] #3
[485]	430	e = [-1.0, -2.0]
[353]	431	f = [1.0, -2.0] #5
	432
	433	vertices = [a, b, c, d,e,f]
[485]	434	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
[353]	435
	436
	437	point_coords = [[-2.0, 2.0],
	438	[-1.0, 1.0],
[485]	439	[0.0, 2.0],
[353]	440	[1.0, 1.0],
	441	[2.0, 1.0],
[485]	442	[0.0, 0.0],
[353]	443	[1.0, 0.0],
	444	[0.0, -1.0],
[485]	445	[-0.2, -0.5],
[353]	446	[-0.9, -1.5],
	447	[0.5, -1.9],
[485]	448	[3.0, 1.0]]
[321]	449
[353]	450	interp = Interpolation(vertices, triangles, point_coords)
	451	f = array([linear_function(vertices),2*linear_function(vertices) ])
	452	f = transpose(f)
	453	#print "f",f
	454	z = interp.interpolate(f)
	455	answer = [linear_function(point_coords),
	456
	457	2*linear_function(point_coords) ]
	458	answer = transpose(answer)
	459	#print "z",z
	460	#print "answer",answer
	461	assert allclose(z, answer)
[328]	462
	463	def test_smooth_attributes_to_mesh_function(self):
[454]	464	""" Testing 2 attributes smoothed to the mesh
	465	"""
	466	"""Test multiple attributes
	467	"""
[331]	468
[321]	469	a = [0.0, 0.0]
	470	b = [0.0, 5.0]
	471	c = [5.0, 0.0]
	472	points = [a, b, c]
	473	triangles = [ [1,0,2] ] #bac
	474
	475	d1 = [1.0, 1.0]
	476	d2 = [1.0, 3.0]
[485]	477	d3 = [3.0, 1.0]
	478	z1 = [2, 4]
	479	z2 = [4, 8]
	480	z3 = [4, 8]
	481	data_coords = [d1, d2, d3]
[333]	482	z = [z1, z2, z3]
[321]	483
[333]	484	f = fit_to_mesh(points, triangles, data_coords, z, alpha=0.0)
[331]	485	answer = [[0, 0], [5., 10.], [5., 10.]]
	486
[321]	487	assert allclose(f, answer)
[331]	488
	489
	490
	491	#Tests of smoothing matrix
	492	def test_smoothing_matrix_one_triangle(self):
	493	from Numeric import dot
	494	a = [0.0, 0.0]
	495	b = [0.0, 2.0]
	496	c = [2.0,0.0]
	497	points = [a, b, c]
[321]	498
[331]	499	vertices = [ [1,0,2] ] #bac
	500
	501	interp = Interpolation(points, vertices)
	502
[424]	503	assert allclose(interp.get_D(), [[1, -0.5, -0.5],
[331]	504	[-0.5, 0.5, 0],
	505	[-0.5, 0, 0.5]])
	506
	507	#Define f(x,y) = x
	508	f = array([0,0,2]) #Value at global vertex 2
	509
	510	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	511	# int 1 dx dy = area = 2
[424]	512	assert dot(dot(f, interp.get_D()), f) == 2
[331]	513
	514	#Define f(x,y) = y
	515	f = array([0,2,0]) #Value at global vertex 1
	516
	517	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	518	# int 1 dx dy = area = 2
[424]	519	assert dot(dot(f, interp.get_D()), f) == 2
[331]	520
	521	#Define f(x,y) = x+y
	522	f = array([0,2,2]) #Values at global vertex 1 and 2
	523
	524	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	525	# int 2 dx dy = 2*area = 4
[424]	526	assert dot(dot(f, interp.get_D()), f) == 4
[331]	527
	528
	529
	530	def test_smoothing_matrix_more_triangles(self):
	531	from Numeric import dot
	532
	533	a = [0.0, 0.0]
	534	b = [0.0, 2.0]
	535	c = [2.0,0.0]
	536	d = [0.0, 4.0]
	537	e = [2.0, 2.0]
	538	f = [4.0,0.0]
	539
	540	points = [a, b, c, d, e, f]
	541	#bac, bce, ecf, dbe, daf, dae
	542	vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	543
	544	interp = Interpolation(points, vertices)
	545
	546
[424]	547	#assert allclose(interp.get_D(), [[1, -0.5, -0.5],
[331]	548	# [-0.5, 0.5, 0],
	549	# [-0.5, 0, 0.5]])
	550
	551	#Define f(x,y) = x
	552	f = array([0,0,2,0,2,4]) #f evaluated at points a-f
	553
	554	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	555	# int 1 dx dy = total area = 8
[424]	556	assert dot(dot(f, interp.get_D()), f) == 8
[331]	557
	558	#Define f(x,y) = y
	559	f = array([0,2,0,4,2,0]) #f evaluated at points a-f
	560
	561	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	562	# int 1 dx dy = area = 8
[424]	563	assert dot(dot(f, interp.get_D()), f) == 8
[331]	564
	565	#Define f(x,y) = x+y
	566	f = array([0,2,2,4,4,4]) #f evaluated at points a-f
	567
	568	#Check that int (df/dx)2 + (df/dy)2 dx dy =
	569	# int 2 dx dy = 2*area = 16
[424]	570	assert dot(dot(f, interp.get_D()), f) == 16
[331]	571
	572
	573	def test_fit_and_interpolation(self):
	574	from mesh import Mesh
	575
	576	a = [0.0, 0.0]
	577	b = [0.0, 2.0]
	578	c = [2.0, 0.0]
	579	d = [0.0, 4.0]
	580	e = [2.0, 2.0]
	581	f = [4.0, 0.0]
	582
	583	points = [a, b, c, d, e, f]
	584	#bac, bce, ecf, dbe, daf, dae
	585	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	586
	587	#Get (enough) datapoints
	588	data_points = [[ 0.66666667, 0.66666667],
	589	[ 1.33333333, 1.33333333],
	590	[ 2.66666667, 0.66666667],
	591	[ 0.66666667, 2.66666667],
	592	[ 0.0, 1.0],
	593	[ 0.0, 3.0],
	594	[ 1.0, 0.0],
	595	[ 1.0, 1.0],
	596	[ 1.0, 2.0],
	597	[ 1.0, 3.0],
	598	[ 2.0, 1.0],
	599	[ 3.0, 0.0],
	600	[ 3.0, 1.0]]
	601
	602	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	603
	604	z = linear_function(data_points)
	605	answer = linear_function(points)
	606
	607	f = interp.fit(z)
	608
[454]	609	#print "f",f
	610	#print "answer",answer
[331]	611	assert allclose(f, answer)
	612
	613	#Map back
	614	z1 = interp.interpolate(f)
[454]	615	#print "z1\n", z1
	616	#print "z\n",z
[331]	617	assert allclose(z, z1)
	618
	619
	620	def test_smoothing_and_interpolation(self):
	621
	622	a = [0.0, 0.0]
	623	b = [0.0, 2.0]
	624	c = [2.0, 0.0]
	625	d = [0.0, 4.0]
	626	e = [2.0, 2.0]
	627	f = [4.0, 0.0]
	628
	629	points = [a, b, c, d, e, f]
	630	#bac, bce, ecf, dbe, daf, dae
	631	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	632
	633	#Get (too few!) datapoints
	634	data_points = [[ 0.66666667, 0.66666667],
	635	[ 1.33333333, 1.33333333],
	636	[ 2.66666667, 0.66666667],
	637	[ 0.66666667, 2.66666667]]
	638
	639	z = linear_function(data_points)
	640	answer = linear_function(points)
	641
	642	#Make interpolator with too few data points and no smoothing
	643	interp = Interpolation(points, triangles, data_points, alpha=0.0)
	644	#Must raise an exception
	645	try:
	646	f = interp.fit(z)
	647	except:
	648	pass
	649
	650	#Now try with smoothing parameter
	651	interp = Interpolation(points, triangles, data_points, alpha=1.0e-13)
	652
	653	f = interp.fit(z)
	654	#f will be different from answerr due to smoothing
[454]	655	assert allclose(f, answer,atol=5)
[331]	656
	657	#Map back
	658	z1 = interp.interpolate(f)
	659	assert allclose(z, z1)
[346]	660
[642]	661
	662
	663	def test_fit_and_interpolation_with_new_points(self):
	664	"""Fit a surface to one set of points. Then interpolate that surface
	665	using another set of points.
	666	"""
	667	from mesh import Mesh
	668
	669
	670	#Setup mesh used to represent fitted function
	671	a = [0.0, 0.0]
	672	b = [0.0, 2.0]
	673	c = [2.0, 0.0]
	674	d = [0.0, 4.0]
	675	e = [2.0, 2.0]
	676	f = [4.0, 0.0]
	677
	678	points = [a, b, c, d, e, f]
	679	#bac, bce, ecf, dbe, daf, dae
	680	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	681
	682	#Datapoints to fit from
	683	data_points1 = [[ 0.66666667, 0.66666667],
	684	[ 1.33333333, 1.33333333],
	685	[ 2.66666667, 0.66666667],
	686	[ 0.66666667, 2.66666667],
	687	[ 0.0, 1.0],
	688	[ 0.0, 3.0],
	689	[ 1.0, 0.0],
	690	[ 1.0, 1.0],
	691	[ 1.0, 2.0],
	692	[ 1.0, 3.0],
	693	[ 2.0, 1.0],
	694	[ 3.0, 0.0],
	695	[ 3.0, 1.0]]
	696
	697	#Fit surface to mesh
	698	interp = Interpolation(points, triangles, data_points1, alpha=0.0)
	699	z = linear_function(data_points1) #Example z-values
	700	f = interp.fit(z) #Fitted values at vertices
	701
	702
	703
	704	#New datapoints where interpolated values are sought
	705	data_points2 = [[ 0.0, 0.0],
	706	[ 0.5, 0.5],
	707	[ 0.7, 0.7],
	708	[ 1.0, 0.5],
	709	[ 2.0, 0.4],
	710	[ 2.8, 1.2]]
	711
	712
	713	#Build new A matrix based on new points
	714	interp.build_interpolation_matrix_A(data_points2)
	715
	716	#Interpolate using fitted surface
	717	z1 = interp.interpolate(f)
	718
	719	#Desired result
	720	answer = linear_function(data_points2)
	721	assert allclose(z1, answer)
	722
	723
	724
	725
[346]	726	def test_fit_to_mesh_file(self):
	727	from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
	728	export_trianglulation_file
	729	import tempfile
	730	import os
[331]	731
[346]	732	# create a .tsh file, no user outline
	733	mesh_dic = {}
	734	mesh_dic['generatedpointlist'] = [[0.0, 0.0],
	735	[0.0, 5.0],
	736	[5.0, 0.0]]
	737	mesh_dic['generatedtrianglelist'] = [[0, 2, 1]]
	738	mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
	739	mesh_dic['generatedtriangleattributelist'] = [['']]
[429]	740	mesh_dic['generatedpointattributelist'] = [[], [], []]
	741	mesh_dic['generatedpointattributetitlelist'] = []
	742	mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
	743	mesh_dic['generatedsegmentmarkerlist'] = ['external',
	744	'external',
	745	'external']
	746	mesh_file = tempfile.mktemp(".tsh")
	747	export_trianglulation_file(mesh_file,mesh_dic)
	748
	749	# create an .xya file
	750	point_file = tempfile.mktemp(".xya")
	751	fd = open(point_file,'w')
	752	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")
	753	fd.close()
	754
	755	mesh_output_file = "new_trianlge.tsh"
	756	fit_to_mesh_file(mesh_file,
	757	point_file,
	758	mesh_output_file,
	759	alpha = 0.0)
	760	# load in the .tsh file we just wrote
	761	mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
	762
	763	assert allclose(mesh_dic['generatedpointattributelist'],
	764	[[0.0, 0.0],
	765	[5.0, 10.0],
	766	[5.0,10.0]])
	767
	768	self.failUnless(mesh_dic['generatedpointattributetitlelist'] ==
	769	['elevation','stage'],
	770	'test_fit_to_mesh_file failed')
	771
	772	#clean up
	773	os.remove(mesh_file)
	774	os.remove(point_file)
	775
	776	def test_fit_to_mesh_fileII(self):
	777	from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
	778	export_trianglulation_file
	779	import tempfile
	780	import os
	781
	782	# create a .tsh file, no user outline
	783	mesh_dic = {}
	784	mesh_dic['generatedpointlist'] = [[0.0, 0.0],
	785	[0.0, 5.0],
	786	[5.0, 0.0]]
	787	mesh_dic['generatedtrianglelist'] = [[0, 2, 1]]
	788	mesh_dic['generatedsegmentlist'] = [[0, 1], [2, 0], [1, 2]]
	789	mesh_dic['generatedtriangleattributelist'] = [['']]
[393]	790	mesh_dic['generatedpointattributelist'] = [[1,2], [1,2], [1,2]]
	791	mesh_dic['generatedpointattributetitlelist'] = ['density', 'temp']
[346]	792	mesh_dic['generatedtriangleneighborlist'] = [[-1, -1, -1]]
	793	mesh_dic['generatedsegmentmarkerlist'] = ['external',
	794	'external',
	795	'external']
	796	mesh_file = tempfile.mktemp(".tsh")
	797	export_trianglulation_file(mesh_file,mesh_dic)
	798
	799	# create an .xya file
	800	point_file = tempfile.mktemp(".xya")
	801	fd = open(point_file,'w')
[393]	802	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]	803	fd.close()
	804
[484]	805	mesh_output_file = "new_triangle.tsh"
[346]	806	fit_to_mesh_file(mesh_file,
	807	point_file,
	808	mesh_output_file,
	809	alpha = 0.0)
[352]	810	# load in the .tsh file we just wrote
[346]	811	mesh_dic = mesh_file_to_mesh_dictionary(mesh_output_file)
	812
	813	assert allclose(mesh_dic['generatedpointattributelist'],
[393]	814	[[1.0, 2.0,0.0, 0.0],
	815	[1.0, 2.0,5.0, 10.0],
	816	[1.0, 2.0,5.0,10.0]])
[346]	817
[374]	818	self.failUnless(mesh_dic['generatedpointattributetitlelist'] ==
[393]	819	['density', 'temp','elevation','stage'],
[374]	820	'test_fit_to_mesh_file failed')
	821
[346]	822	#clean up
	823	os.remove(mesh_file)
[484]	824	os.remove(mesh_output_file)
[346]	825	os.remove(point_file)
	826
[301]	827	#-------------------------------------------------------------
	828	if __name__ == "__main__":
[642]	829	#suite = unittest.makeSuite(TestCase,'test')
	830
[570]	831	suite = unittest.makeSuite(TestCase,'test')
[705]	832	#suite = unittest.makeSuite(TestCase,'test_arbitrary_datapoints')
[475]	833	runner = unittest.TextTestRunner(verbosity=1)
[301]	834	runner.run(suite)
	835
	836
	837
	838
	839
	840

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

Download in other formats: