Context Navigation

source: inundation/fit_interpolate/test_interpolate.py @ 3231

Last change on this file since 3231 was 3019, checked in by duncan, 19 years ago
interpolate: adding warnings
File size: 47.8 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):
[3019]	159	#Test that data points coinciding with vertices yield a diagonal matrix
	160
[2187]	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):
[3019]	179	#Try datapoints midway on edges -
	180	#each point should affect two matrix entries equally
	181
[2187]	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):
[3019]	199	#Try datapoints on edges -
	200	#each point should affect two matrix entries in proportion
	201
[2187]	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):
[3019]	221	#Try arbitrary datapoints
	222
[2187]	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):
[3019]	240	#Try arbitrary datapoints one outside the triangle.
	241	#That one should be ignored
	242
[2187]	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):
[3019]	437	#Test linear interpolation of known values at vertices to
	438	#new points inside a triangle
	439
[2189]	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	assert allclose(z, answer)
	489
	490
[2879]	491	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	492
	493	#print "z",z
	494	#print "answer",answer
	495	assert allclose(z, answer)
	496
[2189]	497	def test_interpolate_point_outside_of_mesh(self):
[3019]	498	#Test linear interpolation of known values at vertices to
	499	#new points inside a triangle
	500
[2189]	501	a = [0.0, 0.0]
	502	b = [0.0, 5.0]
	503	c = [5.0, 0.0]
	504	d = [5.0, 5.0]
	505
	506	vertices = [a, b, c, d]
	507	triangles = [ [1,0,2], [2,3,1] ] #bac, cdb
	508
	509	#Far away point
	510	d7 = [-1., -1.]
	511
	512	point_coords = [ d7]
	513	interp = Interpolate(vertices, triangles)
	514
	515	#Known values at vertices
	516	#Functions are x+y, x+2y, 2x+y, x-y-5
	517	f = [ [0., 0., 0., -5.], # (0,0)
	518	[5., 10., 5., -10.], # (0,5)
	519	[5., 5., 10.0, 0.], # (5,0)
	520	[10., 15., 15., -5.]] # (5,5)
	521
[3018]	522	z = interp.interpolate(f, point_coords) #, verbose=True)
[2655]	523	answer = array([ [INF, INF, INF, INF]]) # (-1,-1)
[2189]	524
	525	#print "***********"
	526	#print "z",z
	527	#print "answer",answer
	528	#print "***********"
	529
	530	#Should an error message be returned if points are outside
[3018]	531	# of the mesh?
	532	# A warning message is printed, if verbose is on.
[2189]	533
[2655]	534	for i in range(4):
	535	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	536
[2879]	537	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	538
	539	#print "z",z
	540	#print "answer",answer
	541
[2881]	542	for i in range(4):
	543	self.failUnless( z[0,i] == answer[0,i], 'Fail!')
	544
	545
[2189]	546	def test_interpolate_attributes_to_pointsIV(self):
	547	a = [-1.0, 0.0]
	548	b = [3.0, 4.0]
	549	c = [4.0, 1.0]
	550	d = [-3.0, 2.0] #3
	551	e = [-1.0, -2.0]
	552	f = [1.0, -2.0] #5
	553
	554	vertices = [a, b, c, d,e,f]
	555	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	556
	557
	558	point_coords = [[-2.0, 2.0],
	559	[-1.0, 1.0],
	560	[0.0, 2.0],
	561	[1.0, 1.0],
	562	[2.0, 1.0],
	563	[0.0, 0.0],
	564	[1.0, 0.0],
	565	[0.0, -1.0],
	566	[-0.2, -0.5],
	567	[-0.9, -1.5],
	568	[0.5, -1.9],
	569	[3.0, 1.0]]
	570
	571	interp = Interpolate(vertices, triangles)
	572	f = array([linear_function(vertices),2*linear_function(vertices) ])
	573	f = transpose(f)
	574	#print "f",f
	575	z = interp.interpolate(f, point_coords)
	576	answer = [linear_function(point_coords),
	577	2*linear_function(point_coords) ]
	578	answer = transpose(answer)
	579	#print "z",z
	580	#print "answer",answer
	581	assert allclose(z, answer)
	582
[2879]	583	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
[2189]	584
[2879]	585	#print "z",z
	586	#print "answer",answer
	587	assert allclose(z, answer)
	588
[2189]	589	def test_interpolate_blocking(self):
	590	a = [-1.0, 0.0]
	591	b = [3.0, 4.0]
	592	c = [4.0, 1.0]
	593	d = [-3.0, 2.0] #3
	594	e = [-1.0, -2.0]
	595	f = [1.0, -2.0] #5
	596
	597	vertices = [a, b, c, d,e,f]
	598	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	599
	600
	601	point_coords = [[-2.0, 2.0],
	602	[-1.0, 1.0],
	603	[0.0, 2.0],
	604	[1.0, 1.0],
	605	[2.0, 1.0],
	606	[0.0, 0.0],
	607	[1.0, 0.0],
	608	[0.0, -1.0],
	609	[-0.2, -0.5],
	610	[-0.9, -1.5],
	611	[0.5, -1.9],
	612	[3.0, 1.0]]
	613
	614	interp = Interpolate(vertices, triangles)
	615	f = array([linear_function(vertices),2*linear_function(vertices) ])
	616	f = transpose(f)
	617	#print "f",f
	618	for blocking_max in range(len(point_coords)+2):
	619	#if True:
	620	# blocking_max = 5
	621	z = interp.interpolate(f, point_coords,
	622	start_blocking_len=blocking_max)
	623	answer = [linear_function(point_coords),
	624	2*linear_function(point_coords) ]
	625	answer = transpose(answer)
	626	#print "z",z
	627	#print "answer",answer
	628	assert allclose(z, answer)
[2660]	629
	630	f = array([linear_function(vertices),2*linear_function(vertices),
	631	2*linear_function(vertices) - 100 ])
	632	f = transpose(f)
	633	#print "f",f
	634	for blocking_max in range(len(point_coords)+2):
	635	#if True:
	636	# blocking_max = 5
	637	z = interp.interpolate(f, point_coords,
	638	start_blocking_len=blocking_max)
	639	answer = array([linear_function(point_coords),
	640	2*linear_function(point_coords) ,
	641	2*linear_function(point_coords)-100 ])
	642	z = transpose(z)
	643	#print "z",z
	644	#print "answer",answer
	645	assert allclose(z, answer)
[2189]	646
[2787]	647	def test_interpolate_geo_spatial(self):
	648	a = [-1.0, 0.0]
	649	b = [3.0, 4.0]
	650	c = [4.0, 1.0]
	651	d = [-3.0, 2.0] #3
	652	e = [-1.0, -2.0]
	653	f = [1.0, -2.0] #5
	654
	655	vertices = [a, b, c, d,e,f]
	656	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	657
	658
	659	point_coords_absolute = [[-2.0, 2.0],
	660	[-1.0, 1.0],
	661	[0.0, 2.0],
	662	[1.0, 1.0],
	663	[2.0, 1.0],
	664	[0.0, 0.0],
	665	[1.0, 0.0],
	666	[0.0, -1.0],
	667	[-0.2, -0.5],
	668	[-0.9, -1.5],
	669	[0.5, -1.9],
	670	[3.0, 1.0]]
	671
	672	geo = Geo_reference(57,100, 500)
	673	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	674	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	675
	676	interp = Interpolate(vertices, triangles)
	677	f = array([linear_function(vertices),2*linear_function(vertices) ])
	678	f = transpose(f)
	679	#print "f",f
	680	for blocking_max in range(14):
	681	#if True:
	682	# blocking_max = 5
	683	z = interp.interpolate(f, point_coords,
	684	start_blocking_len=blocking_max)
	685	answer = [linear_function(point_coords.get_data_points( \
	686	absolute = True)),
	687	2*linear_function(point_coords.get_data_points( \
	688	absolute = True)) ]
	689	answer = transpose(answer)
	690	#print "z",z
	691	#print "answer",answer
	692	assert allclose(z, answer)
	693
	694	f = array([linear_function(vertices),2*linear_function(vertices),
	695	2*linear_function(vertices) - 100 ])
	696	f = transpose(f)
	697	#print "f",f
	698	for blocking_max in range(14):
	699	#if True:
	700	# blocking_max = 5
	701	z = interp.interpolate(f, point_coords,
	702	start_blocking_len=blocking_max)
	703	answer = array([linear_function(point_coords.get_data_points( \
	704	absolute = True)),
	705	2*linear_function(point_coords.get_data_points( \
	706	absolute = True)) ,
	707	2*linear_function(point_coords.get_data_points( \
	708	absolute = True))-100 ])
	709	z = transpose(z)
	710	#print "z",z
	711	#print "answer",answer
	712	assert allclose(z, answer)
	713
[2879]	714	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	715
	716	#print "z",z
	717	#print "answer",answer
	718	assert allclose(z, answer)
	719
[2787]	720	def test_interpolate_geo_spatial(self):
	721	a = [-1.0, 0.0]
	722	b = [3.0, 4.0]
	723	c = [4.0, 1.0]
	724	d = [-3.0, 2.0] #3
	725	e = [-1.0, -2.0]
	726	f = [1.0, -2.0] #5
	727
	728	vertices = [a, b, c, d,e,f]
	729	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	730
	731
	732	point_coords_absolute = [[-2.0, 2.0],
	733	[-1.0, 1.0],
	734	[0.0, 2.0],
	735	[1.0, 1.0],
	736	[2.0, 1.0],
	737	[0.0, 0.0],
	738	[1.0, 0.0],
	739	[0.0, -1.0],
	740	[-0.2, -0.5],
	741	[-0.9, -1.5],
	742	[0.5, -1.9],
	743	[3.0, 1.0]]
	744
	745	geo = Geo_reference(57,100, 500)
	746	point_coords = geo.change_points_geo_ref(point_coords_absolute)
	747	point_coords = Geospatial_data(point_coords,geo_reference = geo)
	748
	749	interp = Interpolate(vertices, triangles)
	750	f = array([linear_function(vertices),2*linear_function(vertices) ])
	751	f = transpose(f)
	752	#print "f",f
	753	z = interp.interpolate_block(f, point_coords)
	754	answer = [linear_function(point_coords.get_data_points( \
	755	absolute = True)),
	756	2*linear_function(point_coords.get_data_points( \
	757	absolute = True)) ]
	758	answer = transpose(answer)
	759	#print "z",z
	760	#print "answer",answer
	761	assert allclose(z, answer)
	762
[2879]	763	z = interp.interpolate(f, point_coords, start_blocking_len = 2)
	764
	765	#print "z",z
	766	#print "answer",answer
	767	assert allclose(z, answer)
	768
	769
[2189]	770	def test_interpolate_reuse(self):
	771	a = [-1.0, 0.0]
	772	b = [3.0, 4.0]
	773	c = [4.0, 1.0]
	774	d = [-3.0, 2.0] #3
	775	e = [-1.0, -2.0]
	776	f = [1.0, -2.0] #5
	777
	778	vertices = [a, b, c, d,e,f]
	779	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	780
	781
	782	point_coords = [[-2.0, 2.0],
	783	[-1.0, 1.0],
	784	[0.0, 2.0],
	785	[1.0, 1.0],
	786	[2.0, 1.0],
	787	[0.0, 0.0],
	788	[1.0, 0.0],
	789	[0.0, -1.0],
	790	[-0.2, -0.5],
	791	[-0.9, -1.5],
	792	[0.5, -1.9],
	793	[3.0, 1.0]]
	794
	795	interp = Interpolate(vertices, triangles)
	796	f = array([linear_function(vertices),2*linear_function(vertices) ])
	797	f = transpose(f)
	798	z = interp.interpolate(f, point_coords,
	799	start_blocking_len=20)
	800	answer = [linear_function(point_coords),
	801	2*linear_function(point_coords) ]
	802	answer = transpose(answer)
	803	#print "z",z
	804	#print "answer",answer
	805	assert allclose(z, answer)
[2201]	806	assert allclose(interp._A_can_be_reused, True)
[2189]	807
	808	z = interp.interpolate(f)
	809	assert allclose(z, answer)
	810
	811	# This causes blocking to occur.
	812	z = interp.interpolate(f, start_blocking_len=10)
	813	assert allclose(z, answer)
[2201]	814	assert allclose(interp._A_can_be_reused, False)
[2189]	815
	816	#A is recalculated
	817	z = interp.interpolate(f)
	818	assert allclose(z, answer)
[2201]	819	assert allclose(interp._A_can_be_reused, True)
	820
[2189]	821	interp = Interpolate(vertices, triangles)
	822	#Must raise an exception, no points specified
	823	try:
	824	z = interp.interpolate(f)
	825	except:
	826	pass
	827
	828
[2394]	829
	830	def test_interpolation_interface_time_only(self):
	831
[2655]	832	# Test spatio-temporal interpolation
	833	# Test that spatio temporal function performs the correct
	834	# interpolations in both time and space
	835
[2394]	836
[2655]	837
[2394]	838	#Three timesteps
	839	time = [1.0, 5.0, 6.0]
	840
	841
	842	#One quantity
	843	Q = zeros( (3,6), Float )
	844
	845	#Linear in time and space
	846	a = [0.0, 0.0]
	847	b = [0.0, 2.0]
	848	c = [2.0, 0.0]
	849	d = [0.0, 4.0]
	850	e = [2.0, 2.0]
	851	f = [4.0, 0.0]
	852
	853	points = [a, b, c, d, e, f]
	854
	855	for i, t in enumerate(time):
	856	Q[i, :] = t*linear_function(points)
	857
	858
	859	#Check basic interpolation of one quantity using averaging
	860	#(no interpolation points or spatial info)
[2563]	861	I = Interpolation_function(time, [mean(Q[0,:]),
[2394]	862	mean(Q[1,:]),
	863	mean(Q[2,:])])
	864
	865
	866
	867	#Check temporal interpolation
	868	for i in [0,1,2]:
	869	assert allclose(I(time[i]), mean(Q[i,:]))
	870
	871	#Midway
	872	assert allclose(I( (time[0] + time[1])/2 ),
	873	(I(time[0]) + I(time[1]))/2 )
	874
	875	assert allclose(I( (time[1] + time[2])/2 ),
	876	(I(time[1]) + I(time[2]))/2 )
	877
	878	assert allclose(I( (time[0] + time[2])/2 ),
	879	(I(time[0]) + I(time[2]))/2 )
	880
	881	#1/3
	882	assert allclose(I( (time[0] + time[2])/3 ),
	883	(I(time[0]) + I(time[2]))/3 )
	884
	885
	886	#Out of bounds checks
	887	try:
	888	I(time[0]-1)
	889	except:
	890	pass
	891	else:
	892	raise 'Should raise exception'
	893
	894	try:
	895	I(time[-1]+1)
	896	except:
	897	pass
	898	else:
	899	raise 'Should raise exception'
	900
	901
	902
	903
	904	def test_interpolation_interface_spatial_only(self):
[2655]	905	# Test spatio-temporal interpolation with constant time
	906
[2394]	907	#Three timesteps
	908	time = [1.0, 5.0, 6.0]
	909
	910
	911	#Setup mesh used to represent fitted function
	912	a = [0.0, 0.0]
	913	b = [0.0, 2.0]
	914	c = [2.0, 0.0]
	915	d = [0.0, 4.0]
	916	e = [2.0, 2.0]
	917	f = [4.0, 0.0]
	918
	919	points = [a, b, c, d, e, f]
	920	#bac, bce, ecf, dbe
	921	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	922
	923
	924	#New datapoints where interpolated values are sought
	925	interpolation_points = [[ 0.0, 0.0],
	926	[ 0.5, 0.5],
	927	[ 0.7, 0.7],
	928	[ 1.0, 0.5],
	929	[ 2.0, 0.4],
	930	[ 2.8, 1.2]]
	931
	932
	933	#One quantity linear in space
	934	Q = linear_function(points)
	935
	936
	937	#Check interpolation of one quantity using interpolaton points
[2563]	938	I = Interpolation_function(time, Q,
[2394]	939	vertex_coordinates = points,
	940	triangles = triangles,
	941	interpolation_points = interpolation_points,
	942	verbose = False)
	943
	944
	945	answer = linear_function(interpolation_points)
	946
	947	t = time[0]
	948	for j in range(50): #t in [1, 6]
	949	for id in range(len(interpolation_points)):
	950	assert allclose(I(t, id), answer[id])
	951
	952	t += 0.1
	953
	954
	955	try:
	956	I(1)
	957	except:
	958	pass
	959	else:
	960	raise 'Should raise exception'
	961
	962
	963
	964	def test_interpolation_interface(self):
[2655]	965	# Test spatio-temporal interpolation
	966	# Test that spatio temporal function performs the correct
	967	# interpolations in both time and space
	968
[2394]	969	#Three timesteps
[2659]	970	time = [1.0, 5.0, 6.0]
[2394]	971
	972	#Setup mesh used to represent fitted function
	973	a = [0.0, 0.0]
	974	b = [0.0, 2.0]
	975	c = [2.0, 0.0]
	976	d = [0.0, 4.0]
	977	e = [2.0, 2.0]
	978	f = [4.0, 0.0]
	979
	980	points = [a, b, c, d, e, f]
	981	#bac, bce, ecf, dbe
	982	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	983
	984
	985	#New datapoints where interpolated values are sought
	986	interpolation_points = [[ 0.0, 0.0],
	987	[ 0.5, 0.5],
	988	[ 0.7, 0.7],
	989	[ 1.0, 0.5],
	990	[ 2.0, 0.4],
	991	[ 2.8, 1.2]]
	992
	993	#One quantity
	994	Q = zeros( (3,6), Float )
	995
	996	#Linear in time and space
	997	for i, t in enumerate(time):
	998	Q[i, :] = t*linear_function(points)
	999
	1000	#Check interpolation of one quantity using interpolaton points)
[2563]	1001	I = Interpolation_function(time, Q,
[2394]	1002	vertex_coordinates = points,
	1003	triangles = triangles,
	1004	interpolation_points = interpolation_points,
	1005	verbose = False)
	1006
	1007	answer = linear_function(interpolation_points)
	1008
	1009	t = time[0]
	1010	for j in range(50): #t in [1, 6]
	1011	for id in range(len(interpolation_points)):
	1012	assert allclose(I(t, id), t*answer[id])
	1013	t += 0.1
	1014
	1015	try:
	1016	I(1)
	1017	except:
	1018	pass
	1019	else:
	1020	raise 'Should raise exception'
	1021
	1022
[2750]	1023	def test_interpolation_precompute_points(self):
[2684]	1024	# looking at a discrete mesh
	1025	#
	1026
	1027	#Three timesteps
	1028	time = [0.0, 60.0]
	1029
	1030	#Setup mesh used to represent fitted function
	1031	points = [[ 15., -20.],
	1032	[ 15., 10.],
	1033	[ 0., -20.],
	1034	[ 0., 10.],
	1035	[ 0., -20.],
	1036	[ 15., 10.]]
	1037
	1038	triangles = [[0, 1, 2],
	1039	[3, 4, 5]]
	1040
	1041	#New datapoints where interpolated values are sought
	1042	interpolation_points = [[ 1., 0.], [0.,1.]]
	1043
	1044	#One quantity
	1045	Q = zeros( (2,6), Float )
	1046
	1047	#Linear in time and space
	1048	for i, t in enumerate(time):
	1049	Q[i, :] = t*linear_function(points)
	1050	#print "Q", Q
	1051
	1052
	1053
	1054	interp = Interpolate(points, triangles)
	1055	f = array([linear_function(points),2*linear_function(points) ])
	1056	f = transpose(f)
	1057	#print "f",f
	1058	z = interp.interpolate(f, interpolation_points)
	1059	answer = [linear_function(interpolation_points),
	1060	2*linear_function(interpolation_points) ]
	1061	answer = transpose(answer)
[2750]	1062	#print "z",z
	1063	#print "answer",answer
[2684]	1064	assert allclose(z, answer)
	1065
	1066
	1067	#Check interpolation of one quantity using interpolaton points)
	1068	I = Interpolation_function(time, Q,
	1069	vertex_coordinates = points,
	1070	triangles = triangles,
	1071	interpolation_points = interpolation_points,
	1072	verbose = False)
[2755]	1073
[2750]	1074	#print "I.precomputed_values", I.precomputed_values
[2755]	1075
	1076	msg = 'Interpolation failed'
	1077	assert allclose(I.precomputed_values['Attribute'][1], [60, 60]), msg
	1078	#self.failUnless( I.precomputed_values['Attribute'][1] == 60.0,
	1079	# ' failed')
[2684]	1080
[2659]	1081	def test_interpolation_function_outside_point(self):
	1082	# Test spatio-temporal interpolation
	1083	# Test that spatio temporal function performs the correct
	1084	# interpolations in both time and space
	1085
	1086	#Three timesteps
	1087	time = [1.0, 5.0, 6.0]
[2394]	1088
[2659]	1089	#Setup mesh used to represent fitted function
	1090	a = [0.0, 0.0]
	1091	b = [0.0, 2.0]
	1092	c = [2.0, 0.0]
	1093	d = [0.0, 4.0]
	1094	e = [2.0, 2.0]
	1095	f = [4.0, 0.0]
[2394]	1096
[2659]	1097	points = [a, b, c, d, e, f]
	1098	#bac, bce, ecf, dbe
	1099	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
[2577]	1100
	1101
[2659]	1102	#New datapoints where interpolated values are sought
	1103	interpolation_points = [[ 0.0, 0.0],
	1104	[ 0.5, 0.5],
	1105	[ 0.7, 0.7],
	1106	[ 1.0, 0.5],
	1107	[ 2.0, 0.4],
	1108	[ 545354534, 4354354353]] # outside the mesh
	1109
	1110	#One quantity
	1111	Q = zeros( (3,6), Float )
	1112
	1113	#Linear in time and space
	1114	for i, t in enumerate(time):
	1115	Q[i, :] = t*linear_function(points)
	1116
	1117	#Check interpolation of one quantity using interpolaton points)
	1118	I = Interpolation_function(time, Q,
	1119	vertex_coordinates = points,
	1120	triangles = triangles,
	1121	interpolation_points = interpolation_points,
	1122	verbose = False)
	1123
	1124	answer = linear_function(interpolation_points)
	1125
	1126	t = time[0]
	1127	for j in range(50): #t in [1, 6]
	1128	for id in range(len(interpolation_points)-1):
	1129	assert allclose(I(t, id), t*answer[id])
	1130	t += 0.1
	1131
	1132	# Now test the point outside the mesh
	1133	t = time[0]
	1134	for j in range(50): #t in [1, 6]
	1135	self.failUnless(I(t, 5) == INF, 'Fail!')
	1136	t += 0.1
	1137
	1138	try:
	1139	I(1)
	1140	except:
	1141	pass
	1142	else:
	1143	raise 'Should raise exception'
	1144
[2995]	1145
	1146	def test_interpolation_function_time(self):
[3019]	1147	#Test a long time series with an error in it (this did cause an
	1148	#error once)
	1149
[2995]	1150
	1151	time = array(\
	1152	[0.00000000e+00, 5.00000000e-02, 1.00000000e-01, 1.50000000e-01,
	1153	2.00000000e-01, 2.50000000e-01, 3.00000000e-01, 3.50000000e-01,
	1154	4.00000000e-01, 4.50000000e-01, 5.00000000e-01, 5.50000000e-01,
	1155	6.00000000e-01, 6.50000000e-01, 7.00000000e-01, 7.50000000e-01,
	1156	8.00000000e-01, 8.50000000e-01, 9.00000000e-01, 9.50000000e-01,
	1157	1.00000000e-00, 1.05000000e+00, 1.10000000e+00, 1.15000000e+00,
	1158	1.20000000e+00, 1.25000000e+00, 1.30000000e+00, 1.35000000e+00,
	1159	1.40000000e+00, 1.45000000e+00, 1.50000000e+00, 1.55000000e+00,
	1160	1.60000000e+00, 1.65000000e+00, 1.70000000e+00, 1.75000000e+00,
	1161	1.80000000e+00, 1.85000000e+00, 1.90000000e+00, 1.95000000e+00,
	1162	2.00000000e+00, 2.05000000e+00, 2.10000000e+00, 2.15000000e+00,
	1163	2.20000000e+00, 2.25000000e+00, 2.30000000e+00, 2.35000000e+00,
	1164	2.40000000e+00, 2.45000000e+00, 2.50000000e+00, 2.55000000e+00,
	1165	2.60000000e+00, 2.65000000e+00, 2.70000000e+00, 2.75000000e+00,
	1166	2.80000000e+00, 2.85000000e+00, 2.90000000e+00, 2.95000000e+00,
	1167	3.00000000e+00, 3.05000000e+00, 9.96920997e+36, 3.15000000e+00,
	1168	3.20000000e+00, 3.25000000e+00, 3.30000000e+00, 3.35000000e+00,
	1169	3.40000000e+00, 3.45000000e+00, 3.50000000e+00, 3.55000000e+00,
	1170	3.60000000e+00, 3.65000000e+00, 3.70000000e+00, 3.75000000e+00,
	1171	3.80000000e+00, 3.85000000e+00, 3.90000000e+00, 3.95000000e+00,
	1172	4.00000000e+00, 4.05000000e+00, 4.10000000e+00, 4.15000000e+00,
	1173	4.20000000e+00, 4.25000000e+00, 4.30000000e+00, 4.35000000e+00,
	1174	4.40000000e+00, 4.45000000e+00, 4.50000000e+00, 4.55000000e+00,
	1175	4.60000000e+00, 4.65000000e+00, 4.70000000e+00, 4.75000000e+00,
	1176	4.80000000e+00, 4.85000000e+00, 4.90000000e+00, 4.95000000e+00,
	1177	5.00000000e+00, 5.05000000e+00, 5.10000000e+00, 5.15000000e+00,
	1178	5.20000000e+00, 5.25000000e+00, 5.30000000e+00, 5.35000000e+00,
	1179	5.40000000e+00, 5.45000000e+00, 5.50000000e+00, 5.55000000e+00,
	1180	5.60000000e+00, 5.65000000e+00, 5.70000000e+00, 5.75000000e+00,
	1181	5.80000000e+00, 5.85000000e+00, 5.90000000e+00, 5.95000000e+00,
	1182	6.00000000e+00, 6.05000000e+00, 6.10000000e+00, 6.15000000e+00,
	1183	6.20000000e+00, 6.25000000e+00, 6.30000000e+00, 6.35000000e+00,
	1184	6.40000000e+00, 6.45000000e+00, 6.50000000e+00, 6.55000000e+00,
	1185	6.60000000e+00, 6.65000000e+00, 6.70000000e+00, 6.75000000e+00,
	1186	6.80000000e+00, 6.85000000e+00, 6.90000000e+00, 6.95000000e+00,
	1187	7.00000000e+00, 7.05000000e+00, 7.10000000e+00, 7.15000000e+00,
	1188	7.20000000e+00, 7.25000000e+00, 7.30000000e+00, 7.35000000e+00,
	1189	7.40000000e+00, 7.45000000e+00, 7.50000000e+00, 7.55000000e+00,
	1190	7.60000000e+00, 7.65000000e+00, 7.70000000e+00, 7.75000000e+00,
	1191	7.80000000e+00, 7.85000000e+00, 7.90000000e+00, 7.95000000e+00,
	1192	8.00000000e+00, 8.05000000e+00, 8.10000000e+00, 8.15000000e+00,
	1193	8.20000000e+00, 8.25000000e+00, 8.30000000e+00, 8.35000000e+00,
	1194	8.40000000e+00, 8.45000000e+00, 8.50000000e+00, 8.55000000e+00,
	1195	8.60000000e+00, 8.65000000e+00, 8.70000000e+00, 8.75000000e+00,
	1196	8.80000000e+00, 8.85000000e+00, 8.90000000e+00, 8.95000000e+00,
	1197	9.00000000e+00, 9.05000000e+00, 9.10000000e+00, 9.15000000e+00,
	1198	9.20000000e+00, 9.25000000e+00, 9.30000000e+00, 9.35000000e+00,
	1199	9.40000000e+00, 9.45000000e+00, 9.50000000e+00, 9.55000000e+00,
	1200	9.60000000e+00, 9.65000000e+00, 9.70000000e+00, 9.75000000e+00,
	1201	9.80000000e+00, 9.85000000e+00, 9.90000000e+00, 9.95000000e+00,
	1202	1.00000000e+01, 1.00500000e+01, 1.01000000e+01, 1.01500000e+01,
	1203	1.02000000e+01, 1.02500000e+01, 1.03000000e+01, 1.03500000e+01,
	1204	1.04000000e+01, 1.04500000e+01, 1.05000000e+01, 1.05500000e+01,
	1205	1.06000000e+01, 1.06500000e+01, 1.07000000e+01, 1.07500000e+01,
	1206	1.08000000e+01, 1.08500000e+01, 1.09000000e+01, 1.09500000e+01,
	1207	1.10000000e+01, 1.10500000e+01, 1.11000000e+01, 1.11500000e+01,
	1208	1.12000000e+01, 1.12500000e+01, 1.13000000e+01, 1.13500000e+01,
	1209	1.14000000e+01, 1.14500000e+01, 1.15000000e+01, 1.15500000e+01,
	1210	1.16000000e+01, 1.16500000e+01, 1.17000000e+01, 1.17500000e+01,
	1211	1.18000000e+01, 1.18500000e+01, 1.19000000e+01, 1.19500000e+01,
	1212	1.20000000e+01, 1.20500000e+01, 1.21000000e+01, 1.21500000e+01,
	1213	1.22000000e+01, 1.22500000e+01, 1.23000000e+01, 1.23500000e+01,
	1214	1.24000000e+01, 1.24500000e+01, 1.25000000e+01, 1.25500000e+01,
	1215	1.26000000e+01, 1.26500000e+01, 1.27000000e+01, 1.27500000e+01,
	1216	1.28000000e+01, 1.28500000e+01, 1.29000000e+01, 1.29500000e+01,
	1217	1.30000000e+01, 1.30500000e+01, 1.31000000e+01, 1.31500000e+01,
	1218	1.32000000e+01, 1.32500000e+01, 1.33000000e+01, 1.33500000e+01,
	1219	1.34000000e+01, 1.34500000e+01, 1.35000000e+01, 1.35500000e+01,
	1220	1.36000000e+01, 1.36500000e+01, 1.37000000e+01, 1.37500000e+01,
	1221	1.38000000e+01, 1.38500000e+01, 1.39000000e+01, 1.39500000e+01,
	1222	1.40000000e+01, 1.40500000e+01, 1.41000000e+01, 1.41500000e+01,
	1223	1.42000000e+01, 1.42500000e+01, 1.43000000e+01, 1.43500000e+01,
	1224	1.44000000e+01, 1.44500000e+01, 1.45000000e+01, 1.45500000e+01,
	1225	1.46000000e+01, 1.46500000e+01, 1.47000000e+01, 1.47500000e+01,
	1226	1.48000000e+01, 1.48500000e+01, 1.49000000e+01, 1.49500000e+01,
	1227	1.50000000e+01, 1.50500000e+01, 1.51000000e+01, 1.51500000e+01,
	1228	1.52000000e+01, 1.52500000e+01, 1.53000000e+01, 1.53500000e+01,
	1229	1.54000000e+01, 1.54500000e+01, 1.55000000e+01, 1.55500000e+01,
	1230	1.56000000e+01, 1.56500000e+01, 1.57000000e+01, 1.57500000e+01,
	1231	1.58000000e+01, 1.58500000e+01, 1.59000000e+01, 1.59500000e+01,
	1232	1.60000000e+01, 1.60500000e+01, 1.61000000e+01, 1.61500000e+01,
	1233	1.62000000e+01, 1.62500000e+01, 1.63000000e+01, 1.63500000e+01,
	1234	1.64000000e+01, 1.64500000e+01, 1.65000000e+01, 1.65500000e+01,
	1235	1.66000000e+01, 1.66500000e+01, 1.67000000e+01, 1.67500000e+01,
	1236	1.68000000e+01, 1.68500000e+01, 1.69000000e+01, 1.69500000e+01,
	1237	1.70000000e+01, 1.70500000e+01, 1.71000000e+01, 1.71500000e+01,
	1238	1.72000000e+01, 1.72500000e+01, 1.73000000e+01, 1.73500000e+01,
	1239	1.74000000e+01, 1.74500000e+01, 1.75000000e+01, 1.75500000e+01,
	1240	1.76000000e+01, 1.76500000e+01, 1.77000000e+01, 1.77500000e+01,
	1241	1.78000000e+01, 1.78500000e+01, 1.79000000e+01, 1.79500000e+01,
	1242	1.80000000e+01, 1.80500000e+01, 1.81000000e+01, 1.81500000e+01,
	1243	1.82000000e+01, 1.82500000e+01, 1.83000000e+01, 1.83500000e+01,
	1244	1.84000000e+01, 1.84500000e+01, 1.85000000e+01, 1.85500000e+01,
	1245	1.86000000e+01, 1.86500000e+01, 1.87000000e+01, 1.87500000e+01,
	1246	1.88000000e+01, 1.88500000e+01, 1.89000000e+01, 1.89500000e+01,
	1247	1.90000000e+01, 1.90500000e+01, 1.91000000e+01, 1.91500000e+01,
	1248	1.92000000e+01, 1.92500000e+01, 1.93000000e+01, 1.93500000e+01,
	1249	1.94000000e+01, 1.94500000e+01, 1.95000000e+01, 1.95500000e+01,
	1250	1.96000000e+01, 1.96500000e+01, 1.97000000e+01, 1.97500000e+01,
	1251	1.98000000e+01, 1.98500000e+01, 1.99000000e+01, 1.99500000e+01,
	1252	2.00000000e+01, 2.00500000e+01, 2.01000000e+01, 2.01500000e+01,
	1253	2.02000000e+01, 2.02500000e+01, 2.03000000e+01, 2.03500000e+01,
	1254	2.04000000e+01, 2.04500000e+01, 2.05000000e+01, 2.05500000e+01,
	1255	2.06000000e+01, 2.06500000e+01, 2.07000000e+01, 2.07500000e+01,
	1256	2.08000000e+01, 2.08500000e+01, 2.09000000e+01, 2.09500000e+01,
	1257	2.10000000e+01, 2.10500000e+01, 2.11000000e+01, 2.11500000e+01,
	1258	2.12000000e+01, 2.12500000e+01, 2.13000000e+01, 2.13500000e+01,
	1259	2.14000000e+01, 2.14500000e+01, 2.15000000e+01, 2.15500000e+01,
	1260	2.16000000e+01, 2.16500000e+01, 2.17000000e+01, 2.17500000e+01,
	1261	2.18000000e+01, 2.18500000e+01, 2.19000000e+01, 2.19500000e+01,
	1262	2.20000000e+01, 2.20500000e+01, 2.21000000e+01, 2.21500000e+01,
	1263	2.22000000e+01, 2.22500000e+01, 2.23000000e+01, 2.23500000e+01,
	1264	2.24000000e+01, 2.24500000e+01, 2.25000000e+01])
	1265
	1266	#print 'Diff', time[1:] - time[:-1]
	1267
	1268	#Setup mesh used to represent fitted function
	1269	a = [0.0, 0.0]
	1270	b = [0.0, 2.0]
	1271	c = [2.0, 0.0]
	1272	d = [0.0, 4.0]
	1273	e = [2.0, 2.0]
	1274	f = [4.0, 0.0]
	1275
	1276	points = [a, b, c, d, e, f]
	1277	#bac, bce, ecf, dbe
	1278	triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
	1279
	1280
	1281	#New datapoints where interpolated values are sought
	1282	interpolation_points = [[ 0.0, 0.0],
	1283	[ 0.5, 0.5],
	1284	[ 0.7, 0.7],
	1285	[ 1.0, 0.5],
	1286	[ 2.0, 0.4],
	1287	[ 545354534, 4354354353]] # outside the mesh
	1288
	1289	#One quantity
	1290	Q = zeros( (len(time),6), Float )
	1291
	1292	#Linear in time and space
	1293	for i, t in enumerate(time):
	1294	Q[i, :] = t*linear_function(points)
	1295
	1296	#Check interpolation of one quantity using interpolaton points)
	1297	try:
	1298	I = Interpolation_function(time, Q,
	1299	vertex_coordinates = points,
	1300	triangles = triangles,
	1301	interpolation_points = interpolation_points,
	1302	verbose = False)
	1303	except:
	1304	pass
	1305	else:
	1306	raise 'Should raise exception due to time being non-monotoneous'
	1307
	1308
[2655]	1309	def test_points_outside_the_polygon(self):
	1310	a = [-1.0, 0.0]
	1311	b = [3.0, 4.0]
	1312	c = [4.0, 1.0]
	1313	d = [-3.0, 2.0] #3
	1314	e = [-1.0, -2.0]
	1315	f = [1.0, -2.0] #5
	1316
	1317	vertices = [a, b, c, d,e,f]
	1318	triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
	1319
	1320	point_coords = [[-2.0, 2.0],
	1321	[-1.0, 1.0],
	1322	[9999.0, 9999.0], # point Outside poly
	1323	[-9999.0, 1.0], # point Outside poly
	1324	[2.0, 1.0],
	1325	[0.0, 0.0],
	1326	[1.0, 0.0],
	1327	[0.0, -1.0],
	1328	[-0.2, -0.5],
	1329	[-0.9, -1.5],
	1330	[0.5, -1.9],
	1331	[999999, 9999999]] # point Outside poly
	1332	geo_data = Geospatial_data(data_points = point_coords)
	1333
	1334	interp = Interpolate(vertices, triangles)
	1335	f = array([linear_function(vertices),2*linear_function(vertices) ])
	1336	f = transpose(f)
	1337	#print "f",f
	1338	z = interp.interpolate(f, geo_data)
	1339	#z = interp.interpolate(f, point_coords)
	1340	answer = [linear_function(point_coords),
	1341	2*linear_function(point_coords) ]
	1342	answer = transpose(answer)
	1343	answer[2,:] = [INF, INF]
	1344	answer[3,:] = [INF, INF]
	1345	answer[11,:] = [INF, INF]
	1346	#print "z",z
	1347	#print "answer _ fixed",answer
	1348	assert allclose(z[0:1], answer[0:1])
	1349	assert allclose(z[4:10], answer[4:10])
	1350	for i in [2,3,11]:
	1351	self.failUnless( z[i,1] == answer[11,1], 'Fail!')
	1352	self.failUnless( z[i,0] == answer[11,0], 'Fail!')
	1353
[2187]	1354	#-------------------------------------------------------------
	1355	if __name__ == "__main__":
[2655]	1356
[2651]	1357	suite = unittest.makeSuite(Test_Interpolate,'test')
[2684]	1358	#suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_precompute_points')
[2187]	1359	runner = unittest.TextTestRunner(verbosity=1)
	1360	runner.run(suite)
	1361
	1362
	1363
	1364
	1365

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

Download in other formats: