Context Navigation

source: inundation/fit_interpolate/test_interpolate.py @ 2891

Last change on this file since 2891 was 2881, checked in by duncan, 19 years ago
add blocking to tests
File size: 37.6 KB

Rev	Line
[2187]	1	#!/usr/bin/env python
	2
	3	#TEST
[2394]	4
	5	#import time, os
	6
	7
[2187]	8	import sys
[2394]	9	import os
[2187]	10	import unittest
	11	from math import sqrt
[2394]	12	import tempfile
[2187]	13
[2394]	14	from Scientific.IO.NetCDF import NetCDFFile
	15	from Numeric import allclose, array, transpose, zeros, Float
[2187]	16
[2394]	17
	18	# ANUGA code imports
[2187]	19	from interpolate import *
	20	from coordinate_transforms.geo_reference import Geo_reference
[2651]	21	from shallow_water import Domain, Transmissive_boundary
[2655]	22	from utilities.numerical_tools import mean, INF
[2394]	23	from data_manager import get_dataobject
[2577]	24	from geospatial_data.geospatial_data import Geospatial_data
[2187]	25
	26	def distance(x, y):
	27	return sqrt( sum( (array(x)-array(y))**2 ))
	28
	29	def linear_function(point):
	30	point = array(point)
	31	return point[:,0]+point[:,1]
	32
	33
	34	class Test_Interpolate(unittest.TestCase):
	35
	36	def setUp(self):
	37
[2394]	38	import time
	39	from mesh_factory import rectangular
	40
	41
	42	#Create basic mesh
	43	points, vertices, boundary = rectangular(2, 2)
	44
	45	#Create shallow water domain
	46	domain = Domain(points, vertices, boundary)
	47	domain.default_order=2
	48	domain.beta_h = 0
	49
	50
	51	#Set some field values
	52	domain.set_quantity('elevation', lambda x,y: -x)
	53	domain.set_quantity('friction', 0.03)
	54
	55
	56	######################
	57	# Boundary conditions
	58	B = Transmissive_boundary(domain)
	59	domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
	60
	61
	62	######################
	63	#Initial condition - with jumps
	64
	65	bed = domain.quantities['elevation'].vertex_values
	66	stage = zeros(bed.shape, Float)
	67
	68	h = 0.3
	69	for i in range(stage.shape[0]):
	70	if i % 2 == 0:
	71	stage[i,:] = bed[i,:] + h
	72	else:
	73	stage[i,:] = bed[i,:]
	74
	75	domain.set_quantity('stage', stage)
	76
	77	domain.distribute_to_vertices_and_edges()
	78
	79
	80	self.domain = domain
	81
	82	C = domain.get_vertex_coordinates()
	83	self.X = C[:,0:6:2].copy()
	84	self.Y = C[:,1:6:2].copy()
	85
	86	self.F = bed
	87
	88
	89
[2187]	90	def tearDown(self):
	91	pass
	92
	93	def test_datapoint_at_centroid(self):
	94	a = [0.0, 0.0]
	95	b = [0.0, 2.0]
	96	c = [2.0,0.0]
	97	points = [a, b, c]
	98	vertices = [ [1,0,2] ] #bac
	99
	100	data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
	101
	102	interp = Interpolate(points, vertices)
	103	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	104	[[1./3, 1./3, 1./3]])
	105
	106	def test_quad_tree(self):
	107	p0 = [-10.0, -10.0]
	108	p1 = [20.0, -10.0]
	109	p2 = [-10.0, 20.0]
	110	p3 = [10.0, 50.0]
	111	p4 = [30.0, 30.0]
	112	p5 = [50.0, 10.0]
	113	p6 = [40.0, 60.0]
	114	p7 = [60.0, 40.0]
	115	p8 = [-66.0, 20.0]
	116	p9 = [10.0, -66.0]
	117
	118	points = [p0, p1, p2, p3, p4, p5, p6, p7, p8, p9]
	119	triangles = [ [0, 1, 2],
	120	[3, 2, 4],
	121	[4, 2, 1],
	122	[4, 1, 5],
	123	[3, 4, 6],
	124	[6, 4, 7],
	125	[7, 4, 5],
	126	[8, 0, 2],
	127	[0, 9, 1]]
	128
	129	data = [ [4,4] ]
	130	interp = Interpolate(points, triangles,
[2201]	131	max_vertices_per_cell = 4)
[2187]	132	#print "PDSG - interp.get_A()", interp.get_A()
	133	answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0.,
	134	0., 0. , 0., 0., 0., 0.]]
	135
	136
	137	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	138	answer)
	139	#interp.set_point_coordinates([[-30, -30]]) #point outside of mesh
	140	#print "PDSG - interp.get_A()", interp.get_A()
	141	data = [[-30, -30]]
	142	answer = [ [ 0.0, 0.0, 0.0, 0.,
	143	0., 0. , 0., 0., 0., 0.]]
	144
	145	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	146	answer)
	147
	148
	149	#point outside of quad tree root cell
	150	#interp.set_point_coordinates([[-70, -70]])
	151	#print "PDSG - interp.get_A()", interp.get_A()
	152	data = [[-70, -70]]
	153	answer = [ [ 0.0, 0.0, 0.0, 0.,
	154	0., 0. , 0., 0., 0., 0.]]
	155	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	156	answer)
	157
	158	def test_datapoints_at_vertices(self):
	159	"""Test that data points coinciding with vertices yield a diagonal matrix
	160	"""
	161
	162	a = [0.0, 0.0]
	163	b = [0.0, 2.0]
	164	c = [2.0,0.0]
	165	points = [a, b, c]
	166	vertices = [ [1,0,2] ] #bac
	167
	168	data = points #Use data at vertices
	169
	170	interp = Interpolate(points, vertices)
	171	answer = [[1., 0., 0.],
	172	[0., 1., 0.],
	173	[0., 0., 1.]]
	174	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	175	answer)
	176
	177
	178	def test_datapoints_on_edge_midpoints(self):
	179	"""Try datapoints midway on edges -
	180	each point should affect two matrix entries equally
	181	"""
	182
	183	a = [0.0, 0.0]
	184	b = [0.0, 2.0]
	185	c = [2.0,0.0]
	186	points = [a, b, c]
	187	vertices = [ [1,0,2] ] #bac
	188
	189	data = [ [0., 1.], [1., 0.], [1., 1.] ]
	190	answer = [[0.5, 0.5, 0.0], #Affects vertex 1 and 0
	191	[0.5, 0.0, 0.5], #Affects vertex 0 and 2
	192	[0.0, 0.5, 0.5]]
[2690]	193	interp = Interpolate(points, vertices)
[2187]	194
	195	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	196	answer)
	197
	198	def test_datapoints_on_edges(self):
	199	"""Try datapoints on edges -
	200	each point should affect two matrix entries in proportion
	201	"""
	202
	203	a = [0.0, 0.0]
	204	b = [0.0, 2.0]
	205	c = [2.0,0.0]
	206	points = [a, b, c]
	207	vertices = [ [1,0,2] ] #bac
	208
	209	data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
	210	answer = [[0.25, 0.75, 0.0], #Affects vertex 1 and 0
	211	[0.25, 0.0, 0.75], #Affects vertex 0 and 2
	212	[0.0, 0.25, 0.75]]
	213
[2690]	214	interp = Interpolate(points, vertices)
[2187]	215
	216	assert allclose(interp._build_interpolation_matrix_A(data).todense(),
	217	answer)
	218
	219
	220	def test_arbitrary_datapoints(self):
	221	"""Try arbitrary datapoints
	222	"""
	223
	224	from Numeric import sum
	225
	226	a = [0.0, 0.0]
	227	b = [0.0, 2.0]
	228	c = [2.0,0.0]
	229	points = [a, b, c]
	230	vertices = [ [1,0,2] ] #bac
	231
	232	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
	233
[2690]	234	interp = Interpolate(points, vertices)
[2187]	235	#print "interp.get_A()", interp.get_A()
	236	results = interp._build_interpolation_matrix_A(data).todense()
	237	assert allclose(sum(results, axis=1), 1.0)
	238
[2655]	239	def test_arbitrary_datapoints_some_outside(self):
[2187]	240	"""Try arbitrary datapoints one outside the triangle.
	241	That one should be ignored
	242	"""
	243
	244	from Numeric import sum
	245
	246	a = [0.0, 0.0]
	247	b = [0.0, 2.0]
	248	c = [2.0,0.0]
	249	points = [a, b, c]
	250	vertices = [ [1,0,2] ] #bac
	251
	252	data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
	253
[2690]	254	interp = Interpolate(points, vertices)
[2187]	255	results = interp._build_interpolation_matrix_A(data).todense()
	256	assert allclose(sum(results, axis=1), [1,1,1,0])
	257
	258
	259
	260	# this causes a memory error in scipy.sparse
	261	def test_more_triangles(self):
	262
	263	a = [-1.0, 0.0]
	264	b = [3.0, 4.0]
	265	c = [4.0,1.0]
	266	d = [-3.0, 2.0] #3
	267	e = [-1.0,-2.0]
	268	f = [1.0, -2.0] #5
	269
	270	points = [a, b, c, d,e,f]
	271	triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
	272
	273	#Data points
	274	data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
	275	interp = Interpolate(points, triangles)
	276
	277	answer = [[0.0, 0.0, 0.0, 1.0, 0.0, 0.0], #Affects point d
	278	[0.5, 0.0, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	279	[0.75, 0.25, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	280	[0.0, 0.5, 0.0, 0.5, 0.0, 0.0], #Affects points a and d
	281	[0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
	282	[1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a, e and f
	283
	284
	285	A = interp._build_interpolation_matrix_A(data).todense()
	286	for i in range(A.shape[0]):
	287	for j in range(A.shape[1]):
	288	if not allclose(A[i,j], answer[i][j]):
	289	print i,j,':',A[i,j], answer[i][j]
	290
	291
[2190]	292	#results = interp._build_interpolation_matrix_A(data).todense()
[2187]	293
[2190]	294	assert allclose(A, answer)
[2690]	295
	296	def test_geo_ref(self):
	297	v0 = [0.0, 0.0]
	298	v1 = [0.0, 5.0]
	299	v2 = [5.0, 0.0]
[2187]	300
[2690]	301	vertices_absolute = [v0, v1, v2]
	302	triangles = [ [1,0,2] ] #bac
	303
	304	geo = Geo_reference(57,100, 500)
	305
	306	vertices = geo.change_points_geo_ref(vertices_absolute)
	307	#print "vertices",vertices
	308
	309	d0 = [1.0, 1.0]
	310	d1 = [1.0, 2.0]
	311	d2 = [3.0, 1.0]
	312	point_coords = [ d0, d1, d2]
	313
	314	interp = Interpolate(vertices, triangles, mesh_origin=geo)
	315	f = linear_function(vertices_absolute)
	316	z = interp.interpolate(f, point_coords)
	317	answer = linear_function(point_coords)
	318
	319	#print "z",z
	320	#print "answer",answer
	321	assert allclose(z, answer)
[2879]	322
[2690]	323
[2879]	324	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	325	answer = linear_function(point_coords)
	326
	327	#print "z",z
	328	#print "answer",answer
	329	assert allclose(z, answer)
	330
	331
[2690]	332	def test_Geospatial_verts(self):
	333	v0 = [0.0, 0.0]
	334	v1 = [0.0, 5.0]
	335	v2 = [5.0, 0.0]
	336
	337	vertices_absolute = [v0, v1, v2]
	338	triangles = [ [1,0,2] ] #bac
	339
	340	geo = Geo_reference(57,100, 500)
	341	vertices = geo.change_points_geo_ref(vertices_absolute)
	342	geopoints = Geospatial_data(vertices,geo_reference = geo)
	343	#print "vertices",vertices
	344
	345	d0 = [1.0, 1.0]
	346	d1 = [1.0, 2.0]
	347	d2 = [3.0, 1.0]
	348	point_coords = [ d0, d1, d2]
	349
	350	interp = Interpolate(geopoints, triangles)
	351	f = linear_function(vertices_absolute)
	352	z = interp.interpolate(f, point_coords)
	353	answer = linear_function(point_coords)
	354
	355	#print "z",z
	356	#print "answer",answer
	357	assert allclose(z, answer)
	358
[2879]	359	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	360	answer = linear_function(point_coords)
	361
	362	#print "z",z
	363	#print "answer",answer
	364	assert allclose(z, answer)
	365
[2189]	366	def test_interpolate_attributes_to_points(self):
	367	v0 = [0.0, 0.0]
	368	v1 = [0.0, 5.0]
	369	v2 = [5.0, 0.0]
[2187]	370
[2189]	371	vertices = [v0, v1, v2]
	372	triangles = [ [1,0,2] ] #bac
[2187]	373
[2189]	374	d0 = [1.0, 1.0]
	375	d1 = [1.0, 2.0]
	376	d2 = [3.0, 1.0]
	377	point_coords = [ d0, d1, d2]
	378
[2690]	379	interp = Interpolate(vertices, triangles)
[2189]	380	f = linear_function(vertices)
	381	z = interp.interpolate(f, point_coords)
	382	answer = linear_function(point_coords)
	383
[2690]	384	#print "z",z
	385	#print "answer",answer
[2189]	386	assert allclose(z, answer)
	387
	388
[2879]	389	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	390	answer = linear_function(point_coords)
[2189]	391
[2879]	392	#print "z",z
	393	#print "answer",answer
	394	assert allclose(z, answer)
	395
[2189]	396	def test_interpolate_attributes_to_pointsII(self):
	397	a = [-1.0, 0.0]
	398	b = [3.0, 4.0]
	399	c = [4.0, 1.0]
	400	d = [-3.0, 2.0] #3
	401	e = [-1.0, -2.0]
	402	f = [1.0, -2.0] #5
	403
	404	vertices = [a, b, c, d,e,f]
	405	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	406
	407
	408	point_coords = [[-2.0, 2.0],
	409	[-1.0, 1.0],
	410	[0.0, 2.0],
	411	[1.0, 1.0],
	412	[2.0, 1.0],
	413	[0.0, 0.0],
	414	[1.0, 0.0],
	415	[0.0, -1.0],
	416	[-0.2, -0.5],
	417	[-0.9, -1.5],
	418	[0.5, -1.9],
	419	[3.0, 1.0]]
	420
	421	interp = Interpolate(vertices, triangles)
	422	f = linear_function(vertices)
	423	z = interp.interpolate(f, point_coords)
	424	answer = linear_function(point_coords)
	425	#print "z",z
	426	#print "answer",answer
	427	assert allclose(z, answer)
	428
[2879]	429	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	430	answer = linear_function(point_coords)
	431
	432	#print "z",z
	433	#print "answer",answer
	434	assert allclose(z, answer)
	435
[2189]	436	def test_interpolate_attributes_to_pointsIII(self):
	437	"""Test linear interpolation of known values at vertices to
	438	new points inside a triangle
	439	"""
	440	a = [0.0, 0.0]
	441	b = [0.0, 5.0]
	442	c = [5.0, 0.0]
	443	d = [5.0, 5.0]
	444
	445	vertices = [a, b, c, d]
	446	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	447
	448	#Points within triangle 1
	449	d0 = [1.0, 1.0]
	450	d1 = [1.0, 2.0]
	451	d2 = [3.0, 1.0]
	452
	453	#Point within triangle 2
	454	d3 = [4.0, 3.0]
	455
	456	#Points on common edge
	457	d4 = [2.5, 2.5]
	458	d5 = [4.0, 1.0]
	459
	460	#Point on common vertex
	461	d6 = [0., 5.]
	462
	463	point_coords = [d0, d1, d2, d3, d4, d5, d6]
	464
	465	interp = Interpolate(vertices, triangles)
	466
	467	#Known values at vertices
	468	#Functions are x+y, x+2y, 2x+y, x-y-5
	469	f = [ [0., 0., 0., -5.], # (0,0)
	470	[5., 10., 5., -10.], # (0,5)
	471	[5., 5., 10.0, 0.], # (5,0)
	472	[10., 15., 15., -5.]] # (5,5)
	473
	474	z = interp.interpolate(f, point_coords)
	475	answer = [ [2., 3., 3., -5.], # (1,1)
	476	[3., 5., 4., -6.], # (1,2)
	477	[4., 5., 7., -3.], # (3,1)
	478	[7., 10., 11., -4.], # (4,3)
	479	[5., 7.5, 7.5, -5.], # (2.5, 2.5)
	480	[5., 6., 9., -2.], # (4,1)
	481	[5., 10., 5., -10.]] # (0,5)
	482
	483	#print "***********"
	484	#print "z",z
	485	#print "answer",answer
	486	#print "***********"
	487
	488	#Should an error message be returned if points are outside
	489	# of the mesh? Not currently.
	490
	491	assert allclose(z, answer)
	492
	493
[2879]	494	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	495
	496	#print "z",z
	497	#print "answer",answer
	498	assert allclose(z, answer)
	499
[2189]	500	def test_interpolate_point_outside_of_mesh(self):
	501	"""Test linear interpolation of known values at vertices to
	502	new points inside a triangle
	503	"""
	504	a = [0.0, 0.0]
	505	b = [0.0, 5.0]
	506	c = [5.0, 0.0]
	507	d = [5.0, 5.0]
	508
	509	vertices = [a, b, c, d]
	510	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	511
	512	#Far away point
	513	d7 = [-1., -1.]
	514
	515	point_coords = [ d7]
	516	interp = Interpolate(vertices, triangles)
	517
	518	#Known values at vertices
	519	#Functions are x+y, x+2y, 2x+y, x-y-5
	520	f = [ [0., 0., 0., -5.], # (0,0)
	521	[5., 10., 5., -10.], # (0,5)
	522	[5., 5., 10.0, 0.], # (5,0)
	523	[10., 15., 15., -5.]] # (5,5)
	524
	525	z = interp.interpolate(f, point_coords)
[2655]	526	answer = array([ [INF, INF, INF, INF]]) # (-1,-1)
[2189]	527
	528	#print "***********"
	529	#print "z",z
	530	#print "answer",answer
	531	#print "***********"
	532
	533	#Should an error message be returned if points are outside
	534	# of the mesh? Not currently.
	535
[2655]	536	for i in range(4):
	537	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	538
[2879]	539	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	540
	541	#print "z",z
	542	#print "answer",answer
	543
[2881]	544	for i in range(4):
	545	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	546
	547
[2189]	548	def test_interpolate_attributes_to_pointsIV(self):
	549	a = [-1.0, 0.0]
	550	b = [3.0, 4.0]
	551	c = [4.0, 1.0]
	552	d = [-3.0, 2.0] #3
	553	e = [-1.0, -2.0]
	554	f = [1.0, -2.0] #5
	555
	556	vertices = [a, b, c, d,e,f]
	557	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	558
	559
	560	point_coords = [[-2.0, 2.0],
	561	[-1.0, 1.0],
	562	[0.0, 2.0],
	563	[1.0, 1.0],
	564	[2.0, 1.0],
	565	[0.0, 0.0],
	566	[1.0, 0.0],
	567	[0.0, -1.0],
	568	[-0.2, -0.5],
	569	[-0.9, -1.5],
	570	[0.5, -1.9],
	571	[3.0, 1.0]]
	572
	573	interp = Interpolate(vertices, triangles)
	574	f = array([linear_function(vertices),2*linear_function(vertices) ])
	575	f = transpose(f)
	576	#print "f",f
	577	z = interp.interpolate(f, point_coords)
	578	answer = [linear_function(point_coords),
	579	2*linear_function(point_coords) ]
	580	answer = transpose(answer)
	581	#print "z",z
	582	#print "answer",answer
	583	assert allclose(z, answer)
	584
[2879]	585	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
[2189]	586
[2879]	587	#print "z",z
	588	#print "answer",answer
	589	assert allclose(z, answer)
	590
[2189]	591	def test_interpolate_blocking(self):
	592	a = [-1.0, 0.0]
	593	b = [3.0, 4.0]
	594	c = [4.0, 1.0]
	595	d = [-3.0, 2.0] #3
	596	e = [-1.0, -2.0]
	597	f = [1.0, -2.0] #5
	598
	599	vertices = [a, b, c, d,e,f]
	600	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	601
	602
	603	point_coords = [[-2.0, 2.0],
	604	[-1.0, 1.0],
	605	[0.0, 2.0],
	606	[1.0, 1.0],
	607	[2.0, 1.0],
	608	[0.0, 0.0],
	609	[1.0, 0.0],
	610	[0.0, -1.0],
	611	[-0.2, -0.5],
	612	[-0.9, -1.5],
	613	[0.5, -1.9],
	614	[3.0, 1.0]]
	615
	616	interp = Interpolate(vertices, triangles)
	617	f = array([linear_function(vertices),2*linear_function(vertices) ])
	618	f = transpose(f)
	619	#print "f",f
	620	for blocking_max in range(len(point_coords)+2):
	621	#if True:
	622	# blocking_max = 5
	623	z = interp.interpolate(f, point_coords,
	624	start_blocking_len=blocking_max)
	625	answer = [linear_function(point_coords),
	626	2*linear_function(point_coords) ]
	627	answer = transpose(answer)
	628	#print "z",z
	629	#print "answer",answer
	630	assert allclose(z, answer)
[2660]	631
	632	f = array([linear_function(vertices),2*linear_function(vertices),
	633	2*linear_function(vertices) - 100 ])
	634	f = transpose(f)
	635	#print "f",f
	636	for blocking_max in range(len(point_coords)+2):
	637	#if True:
	638	# blocking_max = 5
	639	z = interp.interpolate(f, point_coords,
	640	start_blocking_len=blocking_max)
	641	answer = array([linear_function(point_coords),
	642	2*linear_function(point_coords) ,
	643	2*linear_function(point_coords)-100 ])
	644	z = transpose(z)
	645	#print "z",z
	646	#print "answer",answer
	647	assert allclose(z, answer)
[2189]	648
[2787]	649	def test_interpolate_geo_spatial(self):
	650	a = [-1.0, 0.0]
	651	b = [3.0, 4.0]
	652	c = [4.0, 1.0]
	653	d = [-3.0, 2.0] #3
	654	e = [-1.0, -2.0]
	655	f = [1.0, -2.0] #5
	656
	657	vertices = [a, b, c, d,e,f]
	658	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	659
	660
	661	point_coords_absolute = [[-2.0, 2.0],
	662	[-1.0, 1.0],
	663	[0.0, 2.0],
	664	[1.0, 1.0],
	665	[2.0, 1.0],
	666	[0.0, 0.0],
	667	[1.0, 0.0],
	668	[0.0, -1.0],
	669	[-0.2, -0.5],
	670	[-0.9, -1.5],
	671	[0.5, -1.9],
	672	[3.0, 1.0]]
	673
	674	geo = Geo_reference(57,100, 500)
	675	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	676	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	677
	678	interp = Interpolate(vertices, triangles)
	679	f = array([linear_function(vertices),2*linear_function(vertices) ])
	680	f = transpose(f)
	681	#print "f",f
	682	for blocking_max in range(14):
	683	#if True:
	684	# blocking_max = 5
	685	z = interp.interpolate(f, point_coords,
	686	start_blocking_len=blocking_max)
	687	answer = [linear_function(point_coords.get_data_points( \
	688	absolute = True)),
	689	2*linear_function(point_coords.get_data_points( \
	690	absolute = True)) ]
	691	answer = transpose(answer)
	692	#print "z",z
	693	#print "answer",answer
	694	assert allclose(z, answer)
	695
	696	f = array([linear_function(vertices),2*linear_function(vertices),
	697	2*linear_function(vertices) - 100 ])
	698	f = transpose(f)
	699	#print "f",f
	700	for blocking_max in range(14):
	701	#if True:
	702	# blocking_max = 5
	703	z = interp.interpolate(f, point_coords,
	704	start_blocking_len=blocking_max)
	705	answer = array([linear_function(point_coords.get_data_points( \
	706	absolute = True)),
	707	2*linear_function(point_coords.get_data_points( \
	708	absolute = True)) ,
	709	2*linear_function(point_coords.get_data_points( \
	710	absolute = True))-100 ])
	711	z = transpose(z)
	712	#print "z",z
	713	#print "answer",answer
	714	assert allclose(z, answer)
	715
[2879]	716	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	717
	718	#print "z",z
	719	#print "answer",answer
	720	assert allclose(z, answer)
	721
[2787]	722	def test_interpolate_geo_spatial(self):
	723	a = [-1.0, 0.0]
	724	b = [3.0, 4.0]
	725	c = [4.0, 1.0]
	726	d = [-3.0, 2.0] #3
	727	e = [-1.0, -2.0]
	728	f = [1.0, -2.0] #5
	729
	730	vertices = [a, b, c, d,e,f]
	731	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	732
	733
	734	point_coords_absolute = [[-2.0, 2.0],
	735	[-1.0, 1.0],
	736	[0.0, 2.0],
	737	[1.0, 1.0],
	738	[2.0, 1.0],
	739	[0.0, 0.0],
	740	[1.0, 0.0],
	741	[0.0, -1.0],
	742	[-0.2, -0.5],
	743	[-0.9, -1.5],
	744	[0.5, -1.9],
	745	[3.0, 1.0]]
	746
	747	geo = Geo_reference(57,100, 500)
	748	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	749	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	750
	751	interp = Interpolate(vertices, triangles)
	752	f = array([linear_function(vertices),2*linear_function(vertices) ])
	753	f = transpose(f)
	754	#print "f",f
	755	z = interp.interpolate_block(f, point_coords)
	756	answer = [linear_function(point_coords.get_data_points( \
	757	absolute = True)),
	758	2*linear_function(point_coords.get_data_points( \
	759	absolute = True)) ]
	760	answer = transpose(answer)
	761	#print "z",z
	762	#print "answer",answer
	763	assert allclose(z, answer)
	764
[2879]	765	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	766
	767	#print "z",z
	768	#print "answer",answer
	769	assert allclose(z, answer)
	770
	771
[2189]	772	def test_interpolate_reuse(self):
	773	a = [-1.0, 0.0]
	774	b = [3.0, 4.0]
	775	c = [4.0, 1.0]
	776	d = [-3.0, 2.0] #3
	777	e = [-1.0, -2.0]
	778	f = [1.0, -2.0] #5
	779
	780	vertices = [a, b, c, d,e,f]
	781	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	782
	783
	784	point_coords = [[-2.0, 2.0],
	785	[-1.0, 1.0],
	786	[0.0, 2.0],
	787	[1.0, 1.0],
	788	[2.0, 1.0],
	789	[0.0, 0.0],
	790	[1.0, 0.0],
	791	[0.0, -1.0],
	792	[-0.2, -0.5],
	793	[-0.9, -1.5],
	794	[0.5, -1.9],
	795	[3.0, 1.0]]
	796
	797	interp = Interpolate(vertices, triangles)
	798	f = array([linear_function(vertices),2*linear_function(vertices) ])
	799	f = transpose(f)
	800	z = interp.interpolate(f, point_coords,
	801	start_blocking_len=20)
	802	answer = [linear_function(point_coords),
	803	2*linear_function(point_coords) ]
	804	answer = transpose(answer)
	805	#print "z",z
	806	#print "answer",answer
	807	assert allclose(z, answer)
[2201]	808	assert allclose(interp._A_can_be_reused, True)
[2189]	809
	810	z = interp.interpolate(f)
	811	assert allclose(z, answer)
	812
	813	# This causes blocking to occur.
	814	z = interp.interpolate(f, start_blocking_len=10)
	815	assert allclose(z, answer)
[2201]	816	assert allclose(interp._A_can_be_reused, False)
[2189]	817
	818	#A is recalculated
	819	z = interp.interpolate(f)
	820	assert allclose(z, answer)
[2201]	821	assert allclose(interp._A_can_be_reused, True)
	822
[2189]	823	interp = Interpolate(vertices, triangles)
	824	#Must raise an exception, no points specified
	825	try:
	826	z = interp.interpolate(f)
	827	except:
	828	pass
	829
	830
[2394]	831
	832	def test_interpolation_interface_time_only(self):
	833
[2655]	834	# Test spatio-temporal interpolation
	835	# Test that spatio temporal function performs the correct
	836	# interpolations in both time and space
	837
[2394]	838
[2655]	839
[2394]	840	#Three timesteps
	841	time = [1.0, 5.0, 6.0]
	842
	843
	844	#One quantity
	845	Q = zeros( (3,6), Float )
	846
	847	#Linear in time and space
	848	a = [0.0, 0.0]
	849	b = [0.0, 2.0]
	850	c = [2.0, 0.0]
	851	d = [0.0, 4.0]
	852	e = [2.0, 2.0]
	853	f = [4.0, 0.0]
	854
	855	points = [a, b, c, d, e, f]
	856
	857	for i, t in enumerate(time):
	858	Q[i, :] = t*linear_function(points)
	859
	860
	861	#Check basic interpolation of one quantity using averaging
	862	#(no interpolation points or spatial info)
[2563]	863	I = Interpolation_function(time, [mean(Q[0,:]),
[2394]	864	mean(Q[1,:]),
	865	mean(Q[2,:])])
	866
	867
	868
	869	#Check temporal interpolation
	870	for i in [0,1,2]:
	871	assert allclose(I(time[i]), mean(Q[i,:]))
	872
	873	#Midway
	874	assert allclose(I( (time[0] + time[1])/2 ),
	875	(I(time[0]) + I(time[1]))/2 )
	876
	877	assert allclose(I( (time[1] + time[2])/2 ),
	878	(I(time[1]) + I(time[2]))/2 )
	879
	880	assert allclose(I( (time[0] + time[2])/2 ),
	881	(I(time[0]) + I(time[2]))/2 )
	882
	883	#1/3
	884	assert allclose(I( (time[0] + time[2])/3 ),
	885	(I(time[0]) + I(time[2]))/3 )
	886
	887
	888	#Out of bounds checks
	889	try:
	890	I(time[0]-1)
	891	except:
	892	pass
	893	else:
	894	raise 'Should raise exception'
	895
	896	try:
	897	I(time[-1]+1)
	898	except:
	899	pass
	900	else:
	901	raise 'Should raise exception'
	902
	903
	904
	905
	906	def test_interpolation_interface_spatial_only(self):
[2655]	907	# Test spatio-temporal interpolation with constant time
	908
[2394]	909	#Three timesteps
	910	time = [1.0, 5.0, 6.0]
	911
	912
	913	#Setup mesh used to represent fitted function
	914	a = [0.0, 0.0]
	915	b = [0.0, 2.0]
	916	c = [2.0, 0.0]
	917	d = [0.0, 4.0]
	918	e = [2.0, 2.0]
	919	f = [4.0, 0.0]
	920
	921	points = [a, b, c, d, e, f]
	922	#bac, bce, ecf, dbe
	923	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	924
	925
	926	#New datapoints where interpolated values are sought
	927	interpolation_points = [[ 0.0, 0.0],
	928	[ 0.5, 0.5],
	929	[ 0.7, 0.7],
	930	[ 1.0, 0.5],
	931	[ 2.0, 0.4],
	932	[ 2.8, 1.2]]
	933
	934
	935	#One quantity linear in space
	936	Q = linear_function(points)
	937
	938
	939	#Check interpolation of one quantity using interpolaton points
[2563]	940	I = Interpolation_function(time, Q,
[2394]	941	vertex_coordinates = points,
	942	triangles = triangles,
	943	interpolation_points = interpolation_points,
	944	verbose = False)
	945
	946
	947	answer = linear_function(interpolation_points)
	948
	949	t = time[0]
	950	for j in range(50): #t in [1, 6]
	951	for id in range(len(interpolation_points)):
	952	assert allclose(I(t, id), answer[id])
	953
	954	t += 0.1
	955
	956
	957	try:
	958	I(1)
	959	except:
	960	pass
	961	else:
	962	raise 'Should raise exception'
	963
	964
	965
	966	def test_interpolation_interface(self):
[2655]	967	# Test spatio-temporal interpolation
	968	# Test that spatio temporal function performs the correct
	969	# interpolations in both time and space
	970
[2394]	971	#Three timesteps
[2659]	972	time = [1.0, 5.0, 6.0]
[2394]	973
	974	#Setup mesh used to represent fitted function
	975	a = [0.0, 0.0]
	976	b = [0.0, 2.0]
	977	c = [2.0, 0.0]
	978	d = [0.0, 4.0]
	979	e = [2.0, 2.0]
	980	f = [4.0, 0.0]
	981
	982	points = [a, b, c, d, e, f]
	983	#bac, bce, ecf, dbe
	984	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	985
	986
	987	#New datapoints where interpolated values are sought
	988	interpolation_points = [[ 0.0, 0.0],
	989	[ 0.5, 0.5],
	990	[ 0.7, 0.7],
	991	[ 1.0, 0.5],
	992	[ 2.0, 0.4],
	993	[ 2.8, 1.2]]
	994
	995	#One quantity
	996	Q = zeros( (3,6), Float )
	997
	998	#Linear in time and space
	999	for i, t in enumerate(time):
	1000	Q[i, :] = t*linear_function(points)
	1001
	1002	#Check interpolation of one quantity using interpolaton points)
[2563]	1003	I = Interpolation_function(time, Q,
[2394]	1004	vertex_coordinates = points,
	1005	triangles = triangles,
	1006	interpolation_points = interpolation_points,
	1007	verbose = False)
	1008
	1009	answer = linear_function(interpolation_points)
	1010
	1011	t = time[0]
	1012	for j in range(50): #t in [1, 6]
	1013	for id in range(len(interpolation_points)):
	1014	assert allclose(I(t, id), t*answer[id])
	1015	t += 0.1
	1016
	1017	try:
	1018	I(1)
	1019	except:
	1020	pass
	1021	else:
	1022	raise 'Should raise exception'
	1023
	1024
[2750]	1025	def test_interpolation_precompute_points(self):
[2684]	1026	# looking at a discrete mesh
	1027	#
	1028
	1029	#Three timesteps
	1030	time = [0.0, 60.0]
	1031
	1032	#Setup mesh used to represent fitted function
	1033	points = [[ 15., -20.],
	1034	[ 15., 10.],
	1035	[ 0., -20.],
	1036	[ 0., 10.],
	1037	[ 0., -20.],
	1038	[ 15., 10.]]
	1039
	1040	triangles = [[0, 1, 2],
	1041	[3, 4, 5]]
	1042
	1043	#New datapoints where interpolated values are sought
	1044	interpolation_points = [[ 1., 0.], [0.,1.]]
	1045
	1046	#One quantity
	1047	Q = zeros( (2,6), Float )
	1048
	1049	#Linear in time and space
	1050	for i, t in enumerate(time):
	1051	Q[i, :] = t*linear_function(points)
	1052	#print "Q", Q
	1053
	1054
	1055
	1056	interp = Interpolate(points, triangles)
	1057	f = array([linear_function(points),2*linear_function(points) ])
	1058	f = transpose(f)
	1059	#print "f",f
	1060	z = interp.interpolate(f, interpolation_points)
	1061	answer = [linear_function(interpolation_points),
	1062	2*linear_function(interpolation_points) ]
	1063	answer = transpose(answer)
[2750]	1064	#print "z",z
	1065	#print "answer",answer
[2684]	1066	assert allclose(z, answer)
	1067
	1068
	1069	#Check interpolation of one quantity using interpolaton points)
	1070	I = Interpolation_function(time, Q,
	1071	vertex_coordinates = points,
	1072	triangles = triangles,
	1073	interpolation_points = interpolation_points,
	1074	verbose = False)
[2755]	1075
[2750]	1076	#print "I.precomputed_values", I.precomputed_values
[2755]	1077
	1078	msg = 'Interpolation failed'
	1079	assert allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg
	1080	#self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
	1081	# ' failed')
[2684]	1082
[2659]	1083	def test_interpolation_function_outside_point(self):
	1084	# Test spatio-temporal interpolation
	1085	# Test that spatio temporal function performs the correct
	1086	# interpolations in both time and space
	1087
	1088	#Three timesteps
	1089	time = [1.0, 5.0, 6.0]
[2394]	1090
[2659]	1091	#Setup mesh used to represent fitted function
	1092	a = [0.0, 0.0]
	1093	b = [0.0, 2.0]
	1094	c = [2.0, 0.0]
	1095	d = [0.0, 4.0]
	1096	e = [2.0, 2.0]
	1097	f = [4.0, 0.0]
[2394]	1098
[2659]	1099	points = [a, b, c, d, e, f]
	1100	#bac, bce, ecf, dbe
	1101	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
[2577]	1102
	1103
[2659]	1104	#New datapoints where interpolated values are sought
	1105	interpolation_points = [[ 0.0, 0.0],
	1106	[ 0.5, 0.5],
	1107	[ 0.7, 0.7],
	1108	[ 1.0, 0.5],
	1109	[ 2.0, 0.4],
	1110	[ 545354534, 4354354353]] # outside the mesh
	1111
	1112	#One quantity
	1113	Q = zeros( (3,6), Float )
	1114
	1115	#Linear in time and space
	1116	for i, t in enumerate(time):
	1117	Q[i, :] = t*linear_function(points)
	1118
	1119	#Check interpolation of one quantity using interpolaton points)
	1120	I = Interpolation_function(time, Q,
	1121	vertex_coordinates = points,
	1122	triangles = triangles,
	1123	interpolation_points = interpolation_points,
	1124	verbose = False)
	1125
	1126	answer = linear_function(interpolation_points)
	1127
	1128	t = time[0]
	1129	for j in range(50): #t in [1, 6]
	1130	for id in range(len(interpolation_points)-1):
	1131	assert allclose(I(t, id), t*answer[id])
	1132	t += 0.1
	1133
	1134	# Now test the point outside the mesh
	1135	t = time[0]
	1136	for j in range(50): #t in [1, 6]
	1137	self.failUnless(I(t, 5) == INF, 'Fail!')
	1138	t += 0.1
	1139
	1140	try:
	1141	I(1)
	1142	except:
	1143	pass
	1144	else:
	1145	raise 'Should raise exception'
	1146
[2655]	1147	def test_points_outside_the_polygon(self):
	1148	a = [-1.0, 0.0]
	1149	b = [3.0, 4.0]
	1150	c = [4.0, 1.0]
	1151	d = [-3.0, 2.0] #3
	1152	e = [-1.0, -2.0]
	1153	f = [1.0, -2.0] #5
	1154
	1155	vertices = [a, b, c, d,e,f]
	1156	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	1157
	1158	point_coords = [[-2.0, 2.0],
	1159	[-1.0, 1.0],
	1160	[9999.0, 9999.0], # point Outside poly
	1161	[-9999.0, 1.0], # point Outside poly
	1162	[2.0, 1.0],
	1163	[0.0, 0.0],
	1164	[1.0, 0.0],
	1165	[0.0, -1.0],
	1166	[-0.2, -0.5],
	1167	[-0.9, -1.5],
	1168	[0.5, -1.9],
	1169	[999999, 9999999]] # point Outside poly
	1170	geo_data = Geospatial_data(data_points = point_coords)
	1171
	1172	interp = Interpolate(vertices, triangles)
	1173	f = array([linear_function(vertices),2*linear_function(vertices) ])
	1174	f = transpose(f)
	1175	#print "f",f
	1176	z = interp.interpolate(f, geo_data)
	1177	#z = interp.interpolate(f, point_coords)
	1178	answer = [linear_function(point_coords),
	1179	2*linear_function(point_coords) ]
	1180	answer = transpose(answer)
	1181	answer[2,:] = [INF, INF]
	1182	answer[3,:] = [INF, INF]
	1183	answer[11,:] = [INF, INF]
	1184	#print "z",z
	1185	#print "answer _ fixed",answer
	1186	assert allclose(z[0:1], answer[0:1])
	1187	assert allclose(z[4:10], answer[4:10])
	1188	for i in [2,3,11]:
	1189	self.failUnless( z[i,1] == answer[11,1], 'Fail!')
	1190	self.failUnless( z[i,0] == answer[11,0], 'Fail!')
	1191
[2187]	1192	#-------------------------------------------------------------
	1193	if __name__ == "__main__":
[2655]	1194
[2651]	1195	suite = unittest.makeSuite(Test_Interpolate,'test')
[2684]	1196	#suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_precompute_points')
[2187]	1197	runner = unittest.TextTestRunner(verbosity=1)
	1198	runner.run(suite)
	1199
	1200
	1201
	1202
	1203

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

Download in other formats: