source: inundation/ga/storm_surge/pyvolution/test_util.py @ 947

#!/usr/bin/env python


import unittest
from Numeric import zeros, array, allclose
from math import sqrt, pi

from util import *
from config import epsilon


def test_function(x, y):
    return x+y

class TestCase(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass


    def test_gradient(self):
        x0 = 0.0; y0 = 0.0; z0 = 0.0
        x1 = 1.0; y1 = 0.0; z1 = -1.0
        x2 = 0.0; y2 = 1.0; z2 = 0.0
        
        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)

        assert zx == -1.0
        assert zy == 0.0




    def test_gradient2(self):

        from util import gradient
        x0 = 2.0/3; y0 = 2.0/3
        x1=  8.0/3; y1 = 2.0/3
        x2 = 2.0/3; y2 = 8.0/3

        q0 = 2.0+2.0/3
        q1 = 8.0+2.0/3
        q2 = 2.0+8.0/3
        
        #Gradient of fitted pwl surface    
        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
        
        assert abs(a - 3.0) < epsilon
        assert abs(b - 1.0) < epsilon
    


    def test_that_C_extension_compiles(self):
        FN = 'util_ext.c'
        try:
            import util_ext
        except:
            from compile import compile

            try:
                compile(FN)
            except:
                raise 'Could not compile %s' %FN
            else:
                import util_ext

            
    def test_gradient_C_extension(self):
        from util_ext import gradient as gradient_c

        x0 = 2.0/3; y0 = 2.0/3
        x1=  8.0/3; y1 = 2.0/3
        x2 = 2.0/3; y2 = 8.0/3

        q0 = 2.0+2.0/3
        q1 = 8.0+2.0/3
        q2 = 2.0+8.0/3
        
        #Gradient of fitted pwl surface    
        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
        
        assert abs(a - 3.0) < epsilon
        assert abs(b - 1.0) < epsilon

        
    def test_gradient_C_extension3(self):
        from util_ext import gradient as gradient_c
        
        from RandomArray import uniform, seed
        seed(17, 53)

        x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)

        q0 = uniform(0.0, 10.0, 4)
        q1 = uniform(1.0, 3.0, 4)
        q2 = uniform(7.0, 20.0, 4)        
        

        for i in range(4):
            #Gradient of fitted pwl surface
            from util import gradient_python            
            a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2,
                                    q0[i], q1[i], q2[i])

            #print a_ref, b_ref
            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
                              q0[i], q1[i], q2[i])

            #print a, a_ref, b, b_ref
            assert abs(a - a_ref) < epsilon
            assert abs(b - b_ref) < epsilon
        


    #Geometric        
    #def test_distance(self):
    #    from util import distance#
    #
    #    self.failUnless( distance([4,2],[7,6]) == 5.0,
    #                     'Distance is wrong!')
    #    self.failUnless( allclose(distance([7,6],[9,8]), 2.82842712475),
    #                    'distance is wrong!')
    #    self.failUnless( allclose(distance([9,8],[4,2]), 7.81024967591),
    #                    'distance is wrong!')
    #
    #    self.failUnless( distance([9,8],[4,2]) == distance([4,2],[9,8]),
    #                    'distance is wrong!')        

    def test_angle(self):
        from util import angle

        assert allclose(angle([1.0, 1.0])/pi*180, 45.0)


    def test_anglediff(self):
        from util import anglediff        

        assert allclose(anglediff([0.0, 1.], [1.0, 1.0])/pi*180, 45.0)




    def test_file_function_time(self):
        """Test that File function interpolates correctly
        between given times. No x,y dependency here.
        """

        #Write file
        import os, time
        from config import time_format
        from math import sin, pi

        finaltime = 1200
        filename = 'test_file_function.txt'
        fid = open(filename, 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt
    
        fid.close()

        F = file_function(filename)

        #Now try interpolation
        for i in range(20):
            t = i*10
            q = F(t)

            #Exact linear intpolation
            assert allclose(q[0], 2*t) 
            if i%6 == 0:
                assert allclose(q[1], t**2)
                assert allclose(q[2], sin(t*pi/600))                

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        q = F(t) 
        assert allclose( (120**2 + 60**2)/2, q[1] )
        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )


        t = 100 #Two thirds of the way between between 60 and 120
        q = F(t) 
        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
        
        os.remove(filename)
        

    def test_file_function_time_with_domain(self):
        """Test that File function interpolates correctly
        between given times. No x,y dependency here.
        Use domain with starttime
        """

        #Write file
        import os, time, calendar
        from config import time_format
        from math import sin, pi
        from domain import Domain

        finaltime = 1200
        filename = 'test_file_function.txt'
        fid = open(filename, 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt
   
        fid.close()

        a = [0.0, 0.0]
        b = [4.0, 0.0]
        c = [0.0, 3.0]
       
        points = [a, b, c]
        vertices = [[0,1,2]]
        domain = Domain(points, vertices) 

        #Check that domain.starttime is updated if non-existing
        F = file_function(filename, domain)
        assert allclose(domain.starttime, start)        

        #Check that domain.starttime is updated if too early
        domain.starttime = start - 1                
        F = file_function(filename, domain)
        assert allclose(domain.starttime, start)

        #Check that domain.starttime isn't updated if later
        #domain.starttime = start + 1                
        #F = file_function(filename, domain)
        #assert allclose(domain.starttime, start+1)                

        
       
        #Now try interpolation
        for i in range(20):
            t = i*10
            q = F(t)

            #Exact linear intpolation
            assert allclose(q[0], 2*t) 
            if i%6 == 0:
                assert allclose(q[1], t**2)
                assert allclose(q[2], sin(t*pi/600))                

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        q = F(t) 
        assert allclose( (120**2 + 60**2)/2, q[1] )
        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )


        t = 100 #Two thirds of the way between between 60 and 120
        q = F(t) 
        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
        
        os.remove(filename)
        

    def test_file_function_time_with_domain_different_start(self):
        """Test that File function interpolates correctly
        between given times. No x,y dependency here.
        Use domain with a starttime later than that of file
        """

        #Write file
        import os, time, calendar
        from config import time_format
        from math import sin, pi
        from domain import Domain

        finaltime = 1200
        filename = 'test_file_function.txt'
        fid = open(filename, 'w')
        start = time.mktime(time.strptime('2000', '%Y'))
        dt = 60  #One minute intervals
        t = 0.0
        while t <= finaltime:
            t_string = time.strftime(time_format, time.gmtime(t+start))
            fid.write('%s, %f %f %f\n' %(t_string, 2*t, t**2, sin(t*pi/600)))
            t += dt
   
        fid.close()

        a = [0.0, 0.0]
        b = [4.0, 0.0]
        c = [0.0, 3.0]
       
        points = [a, b, c]
        vertices = [[0,1,2]]
        domain = Domain(points, vertices) 

        #Check that domain.starttime isn't updated if later
        delta = 23
        domain.starttime = start + delta                
        F = file_function(filename, domain)
        assert allclose(domain.starttime, start+delta)                

        
        
       
        #Now try interpolation with delta offset 
        for i in range(20):
            t = i*10
            q = F(t-delta)

            #Exact linear intpolation
            assert allclose(q[0], 2*t) 
            if i%6 == 0:
                assert allclose(q[1], t**2)
                assert allclose(q[2], sin(t*pi/600))                

        #Check non-exact

        t = 90 #Halfway between 60 and 120
        q = F(t-delta) 
        assert allclose( (120**2 + 60**2)/2, q[1] )
        assert allclose( (sin(120*pi/600) + sin(60*pi/600))/2, q[2] )


        t = 100 #Two thirds of the way between between 60 and 120
        q = F(t-delta) 
        assert allclose( 2*120**2/3 + 60**2/3, q[1] )
        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
        
        os.remove(filename)


    def test_spatio_temporal_file_function(self):
        """Test that spatio temporal file function performs the correct
        interpolations in both time and space
        """
        import time
       
        #Create sww file of simple propagation from left to right 
        #through rectangular domain
        from shallow_water import Domain, Dirichlet_boundary
        from mesh_factory import rectangular
        from Numeric import take, concatenate, reshape        

        #Create basic mesh and shallow water domain
        points, vertices, boundary = rectangular(3, 3)
        domain1 = Domain(points, vertices, boundary)
        
        from util import mean
        domain1.reduction = mean 
        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
                              # only one value.

        domain1.default_order = 2
        domain1.store = True    
        domain1.set_datadir('.')
        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
                
        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
        domain1.set_quantity('elevation', 0)
        domain1.set_quantity('friction', 0)
        domain1.set_quantity('stage', 0)        

        # Boundary conditions
        B0 = Dirichlet_boundary([0,0,0])                
        B6 = Dirichlet_boundary([0.6,0,0])      
        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
        domain1.check_integrity()

        finaltime = 8
        #Evolution
        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
            pass
            #domain1.write_time()


        #Now read data from sww and check
        from Scientific.IO.NetCDF import NetCDFFile
        filename = domain1.get_name() + '.' + domain1.format
        fid = NetCDFFile(filename)  

        x = fid.variables['x'][:]
        y = fid.variables['y'][:]        
        stage = fid.variables['stage'][:]
        xmomentum = fid.variables['xmomentum'][:]
        ymomentum = fid.variables['ymomentum'][:]
        time = fid.variables['time'][:]

        #Take stage vertex values at last timestep on diagonal
        #Diagonal is identified by vertices: 0, 5, 10, 15

        timestep = len(time)-1 #Last timestep
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/2
        r1 = (D[1] + D[2])/2         
        r2 = (D[2] + D[3])/2         

        #And the midpoints are found now
        Dx = take(reshape(x, (16,1)), [0,5,10,15])
        Dy = take(reshape(y, (16,1)), [0,5,10,15])                

        diag = concatenate( (Dx, Dy), axis=1)
        d_midpoints = (diag[1:] + diag[:-1])/2

        #Let us see if the file function can find the correct
        #values at the midpoints at the last timestep:
        f = file_function(filename, domain1,
                          interpolation_points = d_midpoints)

        q = f(timestep/10., point_id=0); assert allclose(r0, q)
        q = f(timestep/10., point_id=1); assert allclose(r1, q)        
        q = f(timestep/10., point_id=2); assert allclose(r2, q)


        ##################
        #Now do the same for the first timestep

        timestep = 0 #First timestep
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/2
        r1 = (D[1] + D[2])/2         
        r2 = (D[2] + D[3])/2         

        #Let us see if the file function can find the correct
        #values 
        q = f(0, point_id=0); assert allclose(r0, q)
        q = f(0, point_id=1); assert allclose(r1, q)
        q = f(0, point_id=2); assert allclose(r2, q)


        ##################
        #Now do it again for a timestep in the middle

        timestep = 33 
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0 = (D[0] + D[1])/2
        r1 = (D[1] + D[2])/2         
        r2 = (D[2] + D[3])/2         

        q = f(timestep/10., point_id=0); assert allclose(r0, q)
        q = f(timestep/10., point_id=1); assert allclose(r1, q)        
        q = f(timestep/10., point_id=2); assert allclose(r2, q)


        ##################
        #Now check temporal interpolation
        #Halfway between timestep 15 and 16

        timestep = 15
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0_0 = (D[0] + D[1])/2
        r1_0 = (D[1] + D[2])/2         
        r2_0 = (D[2] + D[3])/2         

        #
        timestep = 16
        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
        D = concatenate( (d_stage, d_uh, d_vh), axis=1)

        #Reference interpolated values at midpoints on diagonal at
        #this timestep are
        r0_1 = (D[0] + D[1])/2
        r1_1 = (D[1] + D[2])/2         
        r2_1 = (D[2] + D[3])/2         

        # The reference values are
        r0 = (r0_0 + r0_1)/2
        r1 = (r1_0 + r1_1)/2
        r2 = (r2_0 + r2_1)/2        

        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)        
        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)

        ##################
        #Finally check interpolation 2 thirds of the way
        #between timestep 15 and 16

        # The reference values are
        r0 = (r0_0 + 2*r0_1)/3
        r1 = (r1_0 + 2*r1_1)/3
        r2 = (r2_0 + 2*r2_1)/3        

        #And the file function gives
        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)   
        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
        
        fid.close()

    def test_xya_ascii(self):
        import time, os
        FN = 'xyatest' + str(time.time()) + '.xya'
        fid = open(FN, 'w')
        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )        
        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )            
        fid.close()
        
        points, attributes = read_xya(FN, format = 'asc')

        assert allclose(points, [ [0,1], [1,0], [1,1] ])
        assert allclose(attributes['a1'], [10,30,40.2])
        assert allclose(attributes['a2'], [20,20,40.3])
        assert allclose(attributes['a3'], [30,10,40.4])        
        
        os.remove(FN)
        
    def test_xya_ascii_w_names(self):
        import time, os
        FN = 'xyatest' + str(time.time()) + '.xya'
        fid = open(FN, 'w')
        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )        
        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )            
        fid.close()
        
        points, attributes = read_xya(FN, format = 'asc')

        assert allclose(points, [ [0,1], [1,0], [1,1] ])

        assert allclose(attributes['a1'], [10,30,40.2])
        assert allclose(attributes['a2'], [20,20,40.3])
        assert allclose(attributes['a3'], [30,10,40.4])        

        
        os.remove(FN)




    #Polygon stuff      
    def test_polygon_function_constants(self):
        p1 = [[0,0], [10,0], [10,10], [0,10]]
        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
        
        f = Polygon_function( [(p1, 1.0)] )     
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
        assert allclose(z, [1,1,0,0])
                
        
        f = Polygon_function( [(p2, 2.0)] )
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
        assert allclose(z, [2,0,0,2])   
                
        
        #Combined       
        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
        assert allclose(z, [2,1,0,2])   
        
        
    def test_polygon_function_callable(self):
        """Check that values passed into Polygon_function can be callable 
        themselves.
        """
        p1 = [[0,0], [10,0], [10,10], [0,10]]
        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
        
        f = Polygon_function( [(p1, test_function)] )   
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
        assert allclose(z, [10,14,0,0])
                
        #Combined       
        f = Polygon_function( [(p1, test_function), (p2, 2.0)] )
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
        assert allclose(z, [2,14,0,2])  
        
        
        #Combined w default     
        f = Polygon_function( [(p1, test_function), (p2, 2.0)], default = 3.14)
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
        assert allclose(z, [2,14,3.14,2])               
        
        
        #Combined w default func        
        f = Polygon_function( [(p1, test_function), (p2, 2.0)], 
                              default = test_function)
        z = f([5, 5, 27, 35], [5, 9, 8, -5]) 
        assert allclose(z, [2,14,35,2])                 
        
        
    def test_point_on_line(self):
                
        #Endpoints first        
        assert point_on_line( 0, 0, 0,0, 1,0 )  
        assert point_on_line( 1, 0, 0,0, 1,0 )          

        #Then points on line    
        assert point_on_line( 0.5, 0, 0,0, 1,0 ) 
        assert point_on_line( 0, 0.5, 0,1, 0,0 )                
        assert point_on_line( 1, 0.5, 1,1, 1,0 )        
        assert point_on_line( 0.5, 0.5, 0,0, 1,1 )              
        
        #Then points not on line                
        assert not point_on_line( 0.5, 0, 0,1, 1,1 )    
        assert not point_on_line( 0, 0.5, 0,0, 1,1 ) 
        
        #From real example that failed
        assert not point_on_line( 40,50, 40,20, 40,40 )
        
        
        #From real example that failed
        assert not point_on_line( 40,19, 40,20, 40,40 ) 
                        
        
        
        
    def test_inside_polygon(self):

        
        #Simplest case: Polygon is the unit square
        polygon = [[0,0], [1,0], [1,1], [0,1]]
 
        assert inside_polygon( (0.5, 0.5), polygon )                 
        assert not inside_polygon( (0.5, 1.5), polygon )   
        assert not inside_polygon( (0.5, -0.5), polygon )       
        assert not inside_polygon( (-0.5, 0.5), polygon )               
        assert not inside_polygon( (1.5, 0.5), polygon )        
        
        #Try point on borders   
        assert inside_polygon( (1., 0.5), polygon, closed=True) 
        assert inside_polygon( (0.5, 1), polygon, closed=True)  
        assert inside_polygon( (0., 0.5), polygon, closed=True) 
        assert inside_polygon( (0.5, 0.), polygon, closed=True) 
        
        assert not inside_polygon( (0.5, 1), polygon, closed=False)     
        assert not inside_polygon( (0., 0.5), polygon, closed=False)    
        assert not inside_polygon( (0.5, 0.), polygon, closed=False)    
        assert not inside_polygon( (1., 0.5), polygon, closed=False)    


    
        #From real example (that failed)
        polygon = [[20,20], [40,20], [40,40], [20,40]]         
        points = [ [40, 50] ]
        res = inside_polygon(points, polygon)
        assert len(res) == 0
        
        polygon = [[20,20], [40,20], [40,40], [20,40]]         
        points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ]
        res = inside_polygon(points, polygon)
        assert len(res) == 2
        assert allclose(res, [0,1])
        
                
              
        #More convoluted and non convex polygon
        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
        assert inside_polygon( (0.5, 0.5), polygon )
        assert inside_polygon( (1, -0.5), polygon ) 
        assert inside_polygon( (1.5, 0), polygon )
        
        assert not inside_polygon( (0.5, 1.5), polygon )        
        assert not inside_polygon( (0.5, -0.5), polygon )               
        
        #Now try the vector formulation returning indices
        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
        res = inside_polygon( points, polygon )
        
        assert allclose( res, [0,1,2] )
        
        #Very convoluted polygon
        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
        assert inside_polygon( (5, 5), polygon )        
        assert inside_polygon( (17, 7), polygon )               
        assert inside_polygon( (27, 2), polygon )                       
        assert inside_polygon( (35, -5), polygon )                      
        assert not inside_polygon( (15, 7), polygon )
        assert not inside_polygon( (24, 3), polygon )
        assert not inside_polygon( (25, -10), polygon ) 
        
        
        
        #Another combination (that failed)
        polygon = [[0,0], [10,0], [10,10], [0,10]]
        assert inside_polygon( (5, 5), polygon )        
        assert inside_polygon( (7, 7), polygon )        
        assert not inside_polygon( (-17, 7), polygon )                  
        assert not inside_polygon( (7, 17), polygon )                   
        assert not inside_polygon( (17, 7), polygon )                   
        assert not inside_polygon( (27, 8), polygon )                   
        assert not inside_polygon( (35, -5), polygon )                  
        
        
                        
        
        #Multiple polygons      
        
        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
                   [10,10], [11,10], [11,11], [10,11], [10,10]]
        assert inside_polygon( (0.5, 0.5), polygon )               
        assert inside_polygon( (10.5, 10.5), polygon )                  
        
        #FIXME: Fails if point is 5.5, 5.5
        assert not inside_polygon( (0, 5.5), polygon )                          

        #Polygon with a hole    
        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
                   [0,0], [1,0], [1,1], [0,1], [0,0]]
                        
        assert inside_polygon( (0, -0.5), polygon )
        assert not inside_polygon( (0.5, 0.5), polygon )        
        

    def test_populate_polygon(self):

        polygon = [[0,0], [1,0], [1,1], [0,1]]
        points = populate_polygon(polygon, 5)

        assert len(points) == 5
        for point in points:
            assert inside_polygon(point, polygon)


        #Very convoluted polygon
        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]

        points = populate_polygon(polygon, 5)

        assert len(points) == 5
        for point in points:
            assert inside_polygon(point, polygon)
            

                      
#-------------------------------------------------------------
if __name__ == "__main__":
    suite = unittest.makeSuite(TestCase,'test')
    runner = unittest.TextTestRunner()
    runner.run(suite)