source: inundation/pyvolution/test_util.py @ 1794

Last change on this file since 1794 was 1671, checked in by ole, 19 years ago

Moved old file function stuff out and re-tested.
File_function is now all in one function which is based on NetCDF

File size: 34.1 KB
Line 
1#!/usr/bin/env python
2
3
4import unittest
5from Numeric import zeros, array, allclose
6from math import sqrt, pi
7
8from util import *
9from config import epsilon
10from data_manager import timefile2swww
11
12
13def test_function(x, y):
14    return x+y
15
16class Test_Util(unittest.TestCase):
17    def setUp(self):
18        pass
19
20    def tearDown(self):
21        pass
22
23
24    def test_gradient(self):
25        x0 = 0.0; y0 = 0.0; z0 = 0.0
26        x1 = 1.0; y1 = 0.0; z1 = -1.0
27        x2 = 0.0; y2 = 1.0; z2 = 0.0
28
29        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
30
31        assert zx == -1.0
32        assert zy == 0.0
33
34    def test_gradient_more(self):
35        x0 = 2.0/3; y0 = 2.0/3
36        x1=  8.0/3; y1 = 2.0/3
37        x2 = 2.0/3; y2 = 8.0/3
38
39        q0 = 2.0+2.0/3
40        q1 = 8.0+2.0/3
41        q2 = 2.0+8.0/3
42
43        #Gradient of fitted pwl surface
44        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
45
46        assert abs(a - 3.0) < epsilon
47        assert abs(b - 1.0) < epsilon
48
49
50    def test_gradient2(self):
51        """Test two-point gradient
52        """
53       
54        x0 = 5.0; y0 = 5.0; z0 = 10.0
55        x1 = 8.0; y1 = 2.0; z1 = 1.0
56        x2 = 8.0; y2 = 8.0; z2 = 10.0
57
58        #Reference
59        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
60        a, b = gradient2(x0, y0, x1, y1, z0, z1)
61
62        assert zx == a
63        assert zy == b
64
65        z2_computed = z0 + a*(x2-x0) + b*(y2-y0)
66        assert z2_computed == z2
67       
68    def test_gradient2_more(self):
69        """Test two-point gradient more
70        """
71        x0 = 2.0; y0 = 2.0
72        x1 = 8.0; y1 = 3.0
73        x2 = 1.0; y2 = 8.0
74
75        q0 = 2.0
76        q1 = 8.0
77        q2 = q0
78
79        #Gradient of fitted pwl surface
80        a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
81        a, b = gradient2(x0, y0, x1, y1, q0, q1)       
82
83        assert a == a_ref
84        assert b == b_ref
85
86
87    def test_that_C_extension_compiles(self):
88        FN = 'util_ext.c'
89        try:
90            import util_ext
91        except:
92            from compile import compile
93
94            try:
95                compile(FN)
96            except:
97                raise 'Could not compile %s' %FN
98            else:
99                import util_ext
100
101
102    def test_gradient_C_extension(self):
103        from util_ext import gradient as gradient_c
104
105        x0 = 2.0/3; y0 = 2.0/3
106        x1=  8.0/3; y1 = 2.0/3
107        x2 = 2.0/3; y2 = 8.0/3
108
109        q0 = 2.0+2.0/3
110        q1 = 8.0+2.0/3
111        q2 = 2.0+8.0/3
112
113        #Gradient of fitted pwl surface
114        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
115
116        assert abs(a - 3.0) < epsilon
117        assert abs(b - 1.0) < epsilon
118
119
120    def test_gradient_C_extension3(self):
121        from util_ext import gradient as gradient_c
122
123        from RandomArray import uniform, seed
124        seed(17, 53)
125
126        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
127
128        q0 = uniform(0.0, 10.0, 4)
129        q1 = uniform(1.0, 3.0, 4)
130        q2 = uniform(7.0, 20.0, 4)
131
132
133        for i in range(4):
134            #Gradient of fitted pwl surface
135            from util import gradient_python
136            a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2,
137                                    q0[i], q1[i], q2[i])
138
139            #print a_ref, b_ref
140            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
141                              q0[i], q1[i], q2[i])
142
143            #print a, a_ref, b, b_ref
144            assert abs(a - a_ref) < epsilon
145            assert abs(b - b_ref) < epsilon
146
147
148
149    #Geometric
150    #def test_distance(self):
151    #    from util import distance#
152    #
153    #    self.failUnless( distance([4,2],[7,6]) == 5.0,
154    #                     'Distance is wrong!')
155    #    self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475),
156    #                    'distance is wrong!')
157    #    self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591),
158    #                    'distance is wrong!')
159    #
160    #    self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
161    #                    'distance is wrong!')
162
163    def test_angle(self):
164        from util import angle
165
166        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)
167
168
169    def test_anglediff(self):
170        from util import anglediff
171
172        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)
173
174    def test_ensure_numeric(self):
175        from util import ensure_numeric
176        from Numeric import ArrayType, Float, array
177
178        A = [1,2,3,4]
179        B = ensure_numeric(A)
180        assert type(B) == ArrayType
181        assert B.typecode() == 'l'
182        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
183
184
185        A = [1,2,3.14,4]
186        B = ensure_numeric(A)
187        assert type(B) == ArrayType
188        assert B.typecode() == 'd'
189        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
190
191
192        A = [1,2,3,4]
193        B = ensure_numeric(A, Float)
194        assert type(B) == ArrayType
195        assert B.typecode() == 'd'
196        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
197
198
199        A = [1,2,3,4]
200        B = ensure_numeric(A, Float)
201        assert type(B) == ArrayType
202        assert B.typecode() == 'd'
203        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
204
205
206        A = array([1,2,3,4])
207        B = ensure_numeric(A)
208        assert type(B) == ArrayType
209        assert B.typecode() == 'l'       
210        assert A == B   
211        assert A is B   #Same object
212
213
214        A = array([1,2,3,4])
215        B = ensure_numeric(A, Float)
216        assert type(B) == ArrayType
217        assert B.typecode() == 'd'       
218        assert A == B   
219        assert A is not B   #Not the same object       
220
221
222
223    def test_file_function_time1(self):
224        """Test that File function interpolates correctly
225        between given times. No x,y dependency here.
226        """
227
228        #Write file
229        import os, time
230        from config import time_format
231        from math import sin, pi
232
233        #Typical ASCII file
234        finaltime = 1200
235        filename = 'test_file_function'
236        fid = open(filename + '.txt', 'w')
237        start = time.mktime(time.strptime('2000', '%Y'))
238        dt = 60  #One minute intervals
239        t = 0.0
240        while t <= finaltime:
241            t_string = time.strftime(time_format, time.gmtime(t+start))
242            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
243            t += dt
244
245        fid.close()
246
247        #Convert ASCII file to NetCDF (Which is what we really like!)
248        timefile2swww(filename)
249
250
251        #Create file function
252        F = file_function(filename + '.sww', quantities = ['Attribute0',
253                                                           'Attribute1',
254                                                           'Attribute2'])
255       
256        #Now try interpolation
257        for i in range(20):
258            t = i*10
259            q = F(t)
260
261            #Exact linear intpolation
262            assert allclose(q[0], 2*t)
263            if i%6 == 0:
264                assert allclose(q[1], t**2)
265                assert allclose(q[2], sin(t*pi/600))
266
267        #Check non-exact
268
269        t = 90 #Halfway between 60 and 120
270        q = F(t)
271        assert allclose( (120**2 + 60**2)/2, q[1] )
272        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
273
274
275        t = 100 #Two thirds of the way between between 60 and 120
276        q = F(t)
277        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
278        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
279
280        os.remove(filename + '.txt')
281        os.remove(filename + '.sww')       
282
283
284       
285    def test_spatio_temporal_file_function(self):
286        """Test that spatio temporal file function performs the correct
287        interpolations in both time and space
288        NetCDF version (x,y,t dependency)       
289        """
290        import time
291
292        #Create sww file of simple propagation from left to right
293        #through rectangular domain
294        from shallow_water import Domain, Dirichlet_boundary
295        from mesh_factory import rectangular
296        from Numeric import take, concatenate, reshape
297
298        #Create basic mesh and shallow water domain
299        points, vertices, boundary = rectangular(3, 3)
300        domain1 = Domain(points, vertices, boundary)
301
302        from util import mean
303        domain1.reduction = mean
304        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
305                              # only one value.
306
307        domain1.default_order = 2
308        domain1.store = True
309        domain1.set_datadir('.')
310        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
311        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
312
313        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
314        domain1.set_quantity('elevation', 0)
315        domain1.set_quantity('friction', 0)
316        domain1.set_quantity('stage', 0)
317
318        # Boundary conditions
319        B0 = Dirichlet_boundary([0,0,0])
320        B6 = Dirichlet_boundary([0.6,0,0])
321        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
322        domain1.check_integrity()
323
324        finaltime = 8
325        #Evolution
326        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
327            pass
328            #domain1.write_time()
329
330
331        #Now read data from sww and check
332        from Scientific.IO.NetCDF import NetCDFFile
333        filename = domain1.get_name() + '.' + domain1.format
334        fid = NetCDFFile(filename)
335
336        x = fid.variables['x'][:]
337        y = fid.variables['y'][:]
338        stage = fid.variables['stage'][:]
339        xmomentum = fid.variables['xmomentum'][:]
340        ymomentum = fid.variables['ymomentum'][:]
341        time = fid.variables['time'][:]
342
343        #Take stage vertex values at last timestep on diagonal
344        #Diagonal is identified by vertices: 0, 5, 10, 15
345
346        timestep = len(time)-1 #Last timestep
347        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
348        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
349        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
350        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
351
352        #Reference interpolated values at midpoints on diagonal at
353        #this timestep are
354        r0 = (D[0] + D[1])/2
355        r1 = (D[1] + D[2])/2
356        r2 = (D[2] + D[3])/2
357
358        #And the midpoints are found now
359        Dx = take(reshape(x, (16,1)), [0,5,10,15])
360        Dy = take(reshape(y, (16,1)), [0,5,10,15])
361
362        diag = concatenate( (Dx, Dy), axis=1)
363        d_midpoints = (diag[1:] + diag[:-1])/2
364
365        #Let us see if the file function can find the correct
366        #values at the midpoints at the last timestep:
367        f = file_function(filename, domain1,
368                          interpolation_points = d_midpoints)
369
370        q = f(timestep/10., point_id=0); assert allclose(r0, q)
371        q = f(timestep/10., point_id=1); assert allclose(r1, q)
372        q = f(timestep/10., point_id=2); assert allclose(r2, q)
373
374
375        ##################
376        #Now do the same for the first timestep
377
378        timestep = 0 #First timestep
379        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
380        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
381        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
382        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
383
384        #Reference interpolated values at midpoints on diagonal at
385        #this timestep are
386        r0 = (D[0] + D[1])/2
387        r1 = (D[1] + D[2])/2
388        r2 = (D[2] + D[3])/2
389
390        #Let us see if the file function can find the correct
391        #values
392        q = f(0, point_id=0); assert allclose(r0, q)
393        q = f(0, point_id=1); assert allclose(r1, q)
394        q = f(0, point_id=2); assert allclose(r2, q)
395
396
397        ##################
398        #Now do it again for a timestep in the middle
399
400        timestep = 33
401        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
402        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
403        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
404        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
405
406        #Reference interpolated values at midpoints on diagonal at
407        #this timestep are
408        r0 = (D[0] + D[1])/2
409        r1 = (D[1] + D[2])/2
410        r2 = (D[2] + D[3])/2
411
412        q = f(timestep/10., point_id=0); assert allclose(r0, q)
413        q = f(timestep/10., point_id=1); assert allclose(r1, q)
414        q = f(timestep/10., point_id=2); assert allclose(r2, q)
415
416
417        ##################
418        #Now check temporal interpolation
419        #Halfway between timestep 15 and 16
420
421        timestep = 15
422        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
423        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
424        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
425        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
426
427        #Reference interpolated values at midpoints on diagonal at
428        #this timestep are
429        r0_0 = (D[0] + D[1])/2
430        r1_0 = (D[1] + D[2])/2
431        r2_0 = (D[2] + D[3])/2
432
433        #
434        timestep = 16
435        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
436        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
437        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
438        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
439
440        #Reference interpolated values at midpoints on diagonal at
441        #this timestep are
442        r0_1 = (D[0] + D[1])/2
443        r1_1 = (D[1] + D[2])/2
444        r2_1 = (D[2] + D[3])/2
445
446        # The reference values are
447        r0 = (r0_0 + r0_1)/2
448        r1 = (r1_0 + r1_1)/2
449        r2 = (r2_0 + r2_1)/2
450
451        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
452        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
453        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
454
455        ##################
456        #Finally check interpolation 2 thirds of the way
457        #between timestep 15 and 16
458
459        # The reference values are
460        r0 = (r0_0 + 2*r0_1)/3
461        r1 = (r1_0 + 2*r1_1)/3
462        r2 = (r2_0 + 2*r2_1)/3
463
464        #And the file function gives
465        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
466        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
467        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
468
469        fid.close()
470        import os
471        os.remove(filename)
472
473       
474
475    def test_spatio_temporal_file_function_time(self):
476        """Test that File function interpolates correctly
477        between given times.
478        NetCDF version (x,y,t dependency)
479        """
480
481        #Create NetCDF (sww) file to be read
482        # x: 0, 5, 10, 15
483        # y: -20, -10, 0, 10
484        # t: 0, 60, 120, ...., 1200
485        #
486        # test quantities (arbitrary but non-trivial expressions):
487        #
488        #   stage     = 3*x - y**2 + 2*t
489        #   xmomentum = exp( -((x-7)**2 + (y+5)**2)/20 ) * t**2
490        #   ymomentum = x**2 + y**2 * sin(t*pi/600)
491
492        #NOTE: Nice test that may render some of the others redundant.
493
494        import os, time
495        from config import time_format
496        from Numeric import sin, pi, exp
497        from mesh_factory import rectangular
498        from shallow_water import Domain
499        import data_manager
500
501        finaltime = 1200
502        filename = 'test_file_function'
503
504        #Create a domain to hold test grid
505        #(0:15, -20:10)
506        points, vertices, boundary =\
507                rectangular(4, 4, 15, 30, origin = (0, -20))
508
509
510        #print 'Number of elements', len(vertices)
511        domain = Domain(points, vertices, boundary)
512        domain.smooth = False
513        domain.default_order = 2
514        domain.set_datadir('.')
515        domain.set_name(filename)
516        domain.store = True
517        domain.format = 'sww'   #Native netcdf visualisation format
518
519        #print 'E', domain.get_extent()
520        #print points
521        start = time.mktime(time.strptime('2000', '%Y'))
522        domain.starttime = start
523
524
525        #Store structure
526        domain.initialise_storage()
527
528        #Compute artificial time steps and store
529        dt = 60  #One minute intervals
530        t = 0.0
531        while t <= finaltime:
532            #Compute quantities
533            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
534            domain.set_quantity('stage', f1)
535
536            f2 = lambda x,y: x+y+t**2
537            domain.set_quantity('xmomentum', f2)
538
539            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
540            domain.set_quantity('ymomentum', f3)
541
542            #Store and advance time
543            domain.time = t
544            domain.store_timestep(domain.conserved_quantities)
545            t += dt
546
547
548        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]]
549
550
551
552        #Deliberately set domain.starttime to too early
553        domain.starttime = start - 1
554
555        #Create file function
556        F = file_function(filename + '.sww', domain,
557                          quantities = domain.conserved_quantities,
558                          interpolation_points = interpolation_points)
559
560        #Check that FF updates fixes domain starttime
561        assert allclose(domain.starttime, start)
562
563        #Check that domain.starttime isn't updated if later
564        domain.starttime = start + 1
565        F = file_function(filename + '.sww', domain,
566                          quantities = domain.conserved_quantities,
567                          interpolation_points = interpolation_points)
568        assert allclose(domain.starttime, start+1)
569        domain.starttime = start
570
571
572        #Check linear interpolation in time
573        F = file_function(filename + '.sww', domain,
574                          quantities = domain.conserved_quantities,
575                          interpolation_points = interpolation_points)               
576        for id in range(len(interpolation_points)):
577            x = interpolation_points[id][0]
578            y = interpolation_points[id][1]
579
580            for i in range(20):
581                t = i*10
582                k = i%6
583
584                if k == 0:
585                    q0 = F(t, point_id=id)
586                    q1 = F(t+60, point_id=id)
587
588
589                q = F(t, point_id=id)
590                #print i, k, t, q
591                #print ' ', q0
592                #print ' ', q1
593                #print 's', (k*q1 + (6-k)*q0)/6
594                #print
595                assert allclose(q, (k*q1 + (6-k)*q0)/6)
596
597
598        #Another check of linear interpolation in time
599        for id in range(len(interpolation_points)):
600            q60 = F(60, point_id=id)
601            q120 = F(120, point_id=id)
602
603            t = 90 #Halfway between 60 and 120
604            q = F(t, point_id=id)
605            assert allclose( (q120+q60)/2, q )
606
607            t = 100 #Two thirds of the way between between 60 and 120
608            q = F(t, point_id=id)
609            assert allclose(q60/3 + 2*q120/3, q)
610
611
612
613        #Check that domain.starttime isn't updated if later than file starttime but earlier
614        #than file end time
615        delta = 23
616        domain.starttime = start + delta
617        F = file_function(filename + '.sww', domain,
618                          quantities = domain.conserved_quantities,
619                          interpolation_points = interpolation_points)
620        assert allclose(domain.starttime, start+delta)
621
622
623
624
625        #Now try interpolation with delta offset
626        for id in range(len(interpolation_points)):           
627            x = interpolation_points[id][0]
628            y = interpolation_points[id][1]
629
630            for i in range(20):
631                t = i*10
632                k = i%6
633
634                if k == 0:
635                    q0 = F(t-delta, point_id=id)
636                    q1 = F(t+60-delta, point_id=id)
637
638                q = F(t-delta, point_id=id)
639                assert allclose(q, (k*q1 + (6-k)*q0)/6)
640
641
642        os.remove(filename + '.sww')
643
644
645
646    def test_file_function_time_with_domain(self):
647        """Test that File function interpolates correctly
648        between given times. No x,y dependency here.
649        Use domain with starttime
650        """
651
652        #Write file
653        import os, time, calendar
654        from config import time_format
655        from math import sin, pi
656        from domain import Domain
657
658        finaltime = 1200
659        filename = 'test_file_function'
660        fid = open(filename + '.txt', 'w')
661        start = time.mktime(time.strptime('2000', '%Y'))
662        dt = 60  #One minute intervals
663        t = 0.0
664        while t <= finaltime:
665            t_string = time.strftime(time_format, time.gmtime(t+start))
666            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
667            t += dt
668
669        fid.close()
670
671
672        #Convert ASCII file to NetCDF (Which is what we really like!)
673        timefile2swww(filename)
674
675
676
677        a = [0.0, 0.0]
678        b = [4.0, 0.0]
679        c = [0.0, 3.0]
680
681        points = [a, b, c]
682        vertices = [[0,1,2]]
683        domain = Domain(points, vertices)
684
685        #Check that domain.starttime is updated if non-existing
686        F = file_function(filename + '.sww', domain)
687       
688        assert allclose(domain.starttime, start)
689
690        #Check that domain.starttime is updated if too early
691        domain.starttime = start - 1
692        F = file_function(filename + '.sww', domain)
693        assert allclose(domain.starttime, start)
694
695        #Check that domain.starttime isn't updated if later
696        domain.starttime = start + 1
697        F = file_function(filename + '.sww', domain)       
698        assert allclose(domain.starttime, start+1)
699
700        domain.starttime = start
701        F = file_function(filename + '.sww', domain,
702                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])
703       
704
705        #print F.T
706        #print F.precomputed_values
707        #print 'F(60)', F(60)
708       
709        #Now try interpolation
710        for i in range(20):
711            t = i*10
712            q = F(t)
713
714            #Exact linear intpolation
715            assert allclose(q[0], 2*t)
716            if i%6 == 0:
717                assert allclose(q[1], t**2)
718                assert allclose(q[2], sin(t*pi/600))
719
720        #Check non-exact
721
722        t = 90 #Halfway between 60 and 120
723        q = F(t)
724        assert allclose( (120**2 + 60**2)/2, q[1] )
725        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
726
727
728        t = 100 #Two thirds of the way between between 60 and 120
729        q = F(t)
730        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
731        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
732
733        os.remove(filename + '.sww')
734        os.remove(filename + '.txt')       
735
736    def test_file_function_time_with_domain_different_start(self):
737        """Test that File function interpolates correctly
738        between given times. No x,y dependency here.
739        Use domain with a starttime later than that of file
740
741        ASCII version
742        """
743
744        #Write file
745        import os, time, calendar
746        from config import time_format
747        from math import sin, pi
748        from domain import Domain
749
750        finaltime = 1200
751        filename = 'test_file_function'
752        fid = open(filename + '.txt', 'w')
753        start = time.mktime(time.strptime('2000', '%Y'))
754        dt = 60  #One minute intervals
755        t = 0.0
756        while t <= finaltime:
757            t_string = time.strftime(time_format, time.gmtime(t+start))
758            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
759            t += dt
760
761        fid.close()
762
763        #Convert ASCII file to NetCDF (Which is what we really like!)
764        timefile2swww(filename)       
765
766        a = [0.0, 0.0]
767        b = [4.0, 0.0]
768        c = [0.0, 3.0]
769
770        points = [a, b, c]
771        vertices = [[0,1,2]]
772        domain = Domain(points, vertices)
773
774        #Check that domain.starttime isn't updated if later than file starttime but earlier
775        #than file end time
776        delta = 23
777        domain.starttime = start + delta
778        #F = file_function(filename + '.sww', domain)
779        F = file_function(filename + '.sww', domain,
780                          quantities = ['Attribute0', 'Attribute1', 'Attribute2'])       
781        assert allclose(domain.starttime, start+delta)
782
783
784
785
786        #Now try interpolation with delta offset
787        for i in range(20):
788            t = i*10
789            q = F(t-delta)
790
791            #Exact linear intpolation
792            assert allclose(q[0], 2*t)
793            if i%6 == 0:
794                assert allclose(q[1], t**2)
795                assert allclose(q[2], sin(t*pi/600))
796
797        #Check non-exact
798
799        t = 90 #Halfway between 60 and 120
800        q = F(t-delta)
801        assert allclose( (120**2 + 60**2)/2, q[1] )
802        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )
803
804
805        t = 100 #Two thirds of the way between between 60 and 120
806        q = F(t-delta)
807        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
808        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
809
810
811        os.remove(filename + '.sww')
812        os.remove(filename + '.txt')               
813
814
815
816    def test_xya_ascii(self):
817        import time, os
818        FN = 'xyatest' + str(time.time()) + '.xya'
819        fid = open(FN, 'w')
820        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )
821        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
822        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
823        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )
824        fid.close()
825
826        points, attributes = read_xya(FN, format = 'asc')
827
828        assert allclose(points, [ [0,1], [1,0], [1,1] ])
829        assert allclose(attributes['a1'], [10,30,40.2])
830        assert allclose(attributes['a2'], [20,20,40.3])
831        assert allclose(attributes['a3'], [30,10,40.4])
832
833        os.remove(FN)
834
835    def test_xya_ascii_w_names(self):
836        import time, os
837        FN = 'xyatest' + str(time.time()) + '.xya'
838        fid = open(FN, 'w')
839        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )
840        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
841        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
842        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )
843        fid.close()
844
845        points, attributes = read_xya(FN, format = 'asc')
846
847        assert allclose(points, [ [0,1], [1,0], [1,1] ])
848
849        assert allclose(attributes['a1'], [10,30,40.2])
850        assert allclose(attributes['a2'], [20,20,40.3])
851        assert allclose(attributes['a3'], [30,10,40.4])
852
853
854        os.remove(FN)
855
856
857
858
859    #Polygon stuff
860    def test_polygon_function_constants(self):
861        p1 = [[0,0], [10,0], [10,10], [0,10]]
862        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
863
864        f = Polygon_function( [(p1, 1.0)] )
865        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
866        assert allclose(z, [1,1,0,0])
867
868
869        f = Polygon_function( [(p2, 2.0)] )
870        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
871        assert allclose(z, [2,0,0,2])
872
873
874        #Combined
875        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
876        z = f([5, 5, 27, 35], [5, 9, 8, -5])
877        assert allclose(z, [2,1,0,2])
878
879
880    def test_polygon_function_callable(self):
881        """Check that values passed into Polygon_function can be callable
882        themselves.
883        """
884        p1 = [[0,0], [10,0], [10,10], [0,10]]
885        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
886
887        f = Polygon_function( [(p1, test_function)] )
888        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
889        assert allclose(z, [10,14,0,0])
890
891        #Combined
892        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
893        z = f([5, 5, 27, 35], [5, 9, 8, -5])
894        assert allclose(z, [2,14,0,2])
895
896
897        #Combined w default
898        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
899        z = f([5, 5, 27, 35], [5, 9, 8, -5])
900        assert allclose(z, [2,14,3.14,2])
901
902
903        #Combined w default func
904        f = Polygon_function( [(p1, test_function), (p2, 2.0)],
905                              default = test_function)
906        z = f([5, 5, 27, 35], [5, 9, 8, -5])
907        assert allclose(z, [2,14,35,2])
908
909
910    def test_point_on_line(self):
911
912        #Endpoints first
913        assert point_on_line( 0, 0, 0,0, 1,0 )
914        assert point_on_line( 1, 0, 0,0, 1,0 )
915
916        #Then points on line
917        assert point_on_line( 0.5, 0, 0,0, 1,0 )
918        assert point_on_line( 0, 0.5, 0,1, 0,0 )
919        assert point_on_line( 1, 0.5, 1,1, 1,0 )
920        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )
921
922        #Then points not on line
923        assert not point_on_line( 0.5, 0, 0,1, 1,1 )
924        assert not point_on_line( 0, 0.5, 0,0, 1,1 )
925
926        #From real example that failed
927        assert not point_on_line( 40,50, 40,20, 40,40 )
928
929
930        #From real example that failed
931        assert not point_on_line( 40,19, 40,20, 40,40 )
932
933
934
935
936    def test_inside_polygon_main(self):
937
938
939        #Simplest case: Polygon is the unit square
940        polygon = [[0,0], [1,0], [1,1], [0,1]]
941
942        assert inside_polygon( (0.5, 0.5), polygon )
943        assert not inside_polygon( (0.5, 1.5), polygon )
944        assert not inside_polygon( (0.5, -0.5), polygon )
945        assert not inside_polygon( (-0.5, 0.5), polygon )
946        assert not inside_polygon( (1.5, 0.5), polygon )
947
948        #Try point on borders
949        assert inside_polygon( (1., 0.5), polygon, closed=True)
950        assert inside_polygon( (0.5, 1), polygon, closed=True)
951        assert inside_polygon( (0., 0.5), polygon, closed=True)
952        assert inside_polygon( (0.5, 0.), polygon, closed=True)
953
954        assert not inside_polygon( (0.5, 1), polygon, closed=False)
955        assert not inside_polygon( (0., 0.5), polygon, closed=False)
956        assert not inside_polygon( (0.5, 0.), polygon, closed=False)
957        assert not inside_polygon( (1., 0.5), polygon, closed=False)
958
959
960
961        #From real example (that failed)
962        polygon = [[20,20], [40,20], [40,40], [20,40]]
963        points = [ [40, 50] ]
964        res = inside_polygon(points, polygon)
965        assert len(res) == 0
966
967        polygon = [[20,20], [40,20], [40,40], [20,40]]
968        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
969        res = inside_polygon(points, polygon)
970        assert len(res) == 2
971        assert allclose(res, [0,1])
972
973
974
975        #More convoluted and non convex polygon
976        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
977        assert inside_polygon( (0.5, 0.5), polygon )
978        assert inside_polygon( (1, -0.5), polygon )
979        assert inside_polygon( (1.5, 0), polygon )
980
981        assert not inside_polygon( (0.5, 1.5), polygon )
982        assert not inside_polygon( (0.5, -0.5), polygon )
983
984
985        #Very convoluted polygon
986        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
987        assert inside_polygon( (5, 5), polygon )
988        assert inside_polygon( (17, 7), polygon )
989        assert inside_polygon( (27, 2), polygon )
990        assert inside_polygon( (35, -5), polygon )
991        assert not inside_polygon( (15, 7), polygon )
992        assert not inside_polygon( (24, 3), polygon )
993        assert not inside_polygon( (25, -10), polygon )
994
995
996
997        #Another combination (that failed)
998        polygon = [[0,0], [10,0], [10,10], [0,10]]
999        assert inside_polygon( (5, 5), polygon )
1000        assert inside_polygon( (7, 7), polygon )
1001        assert not inside_polygon( (-17, 7), polygon )
1002        assert not inside_polygon( (7, 17), polygon )
1003        assert not inside_polygon( (17, 7), polygon )
1004        assert not inside_polygon( (27, 8), polygon )
1005        assert not inside_polygon( (35, -5), polygon )
1006
1007
1008
1009
1010        #Multiple polygons
1011
1012        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
1013                   [10,10], [11,10], [11,11], [10,11], [10,10]]
1014        assert inside_polygon( (0.5, 0.5), polygon )
1015        assert inside_polygon( (10.5, 10.5), polygon )
1016
1017        #FIXME: Fails if point is 5.5, 5.5
1018        assert not inside_polygon( (0, 5.5), polygon )
1019
1020        #Polygon with a hole
1021        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
1022                   [0,0], [1,0], [1,1], [0,1], [0,0]]
1023
1024        assert inside_polygon( (0, -0.5), polygon )
1025        assert not inside_polygon( (0.5, 0.5), polygon )
1026
1027    def test_inside_polygon_vector_version(self):
1028        #Now try the vector formulation returning indices
1029        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1030        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1031        res = inside_polygon( points, polygon, verbose=False )
1032
1033        assert allclose( res, [0,1,2] )
1034
1035    def test_outside_polygon(self):
1036        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1037
1038        assert not outside_polygon( [0.5, 0.5], U )
1039        #evaluate to False as the point 0.5, 0.5 is inside the unit square
1040       
1041        assert outside_polygon( [1.5, 0.5], U )
1042        #evaluate to True as the point 1.5, 0.5 is outside the unit square
1043       
1044        indices = outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
1045        assert allclose( indices, [1] )
1046       
1047        #One more test of vector formulation returning indices
1048        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1049        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1050        res = outside_polygon( points, polygon )
1051
1052        assert allclose( res, [3, 4] )
1053
1054
1055
1056        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1057        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1058        res = outside_polygon( points, polygon )
1059
1060        assert allclose( res, [0, 4, 5] )       
1061     
1062    def test_outside_polygon2(self):
1063        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1064   
1065        assert not outside_polygon( [0.5, 1.0], U, closed = True )
1066        #evaluate to False as the point 0.5, 1.0 is inside the unit square
1067       
1068        assert outside_polygon( [0.5, 1.0], U, closed = False )
1069        #evaluate to True as the point 0.5, 1.0 is outside the unit square
1070
1071    def test_separate_points_by_polygon(self):
1072        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square   
1073
1074        indices, count = separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U )
1075        assert allclose( indices, [0,2,1] )
1076        assert count == 2
1077       
1078        #One more test of vector formulation returning indices
1079        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1080        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1081        res, count = separate_points_by_polygon( points, polygon )
1082
1083        assert allclose( res, [0,1,2,4,3] )
1084        assert count == 3
1085
1086
1087        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
1088        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
1089        res, count = separate_points_by_polygon( points, polygon )
1090
1091        assert allclose( res, [1,2,3,5,4,0] )       
1092        assert count == 3
1093       
1094
1095    def test_populate_polygon(self):
1096
1097        polygon = [[0,0], [1,0], [1,1], [0,1]]
1098        points = populate_polygon(polygon, 5)
1099
1100        assert len(points) == 5
1101        for point in points:
1102            assert inside_polygon(point, polygon)
1103
1104
1105        #Very convoluted polygon
1106        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
1107
1108        points = populate_polygon(polygon, 5)
1109
1110        assert len(points) == 5
1111        for point in points:
1112            assert inside_polygon(point, polygon)
1113
1114
1115
1116#-------------------------------------------------------------
1117if __name__ == "__main__":
1118    suite = unittest.makeSuite(Test_Util,'test')
1119    #suite = unittest.makeSuite(Test_Util,'test_file_function_time')
1120    runner = unittest.TextTestRunner()
1121    runner.run(suite)
1122
1123
1124
1125
Note: See TracBrowser for help on using the repository browser.