Changeset 4816


Ignore:
Timestamp:
Nov 15, 2007, 10:24:02 AM (17 years ago)
Author:
nick
Message:

changed find_optimal_smoothing_parameter to use a passed mesh_file and added some checks.
plus updated split to be more robust
plus added two more test

Location:
anuga_core/source/anuga/geospatial_data
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r4776 r4816  
    1111                        size, shape
    1212#from Array import tolist
    13 from RandomArray import randint, seed
     13from RandomArray import randint, seed, get_seed
    1414from copy import deepcopy
    1515
     
    670670       
    671671        # Set seed if provided, mainly important for unit test!
     672        # plus recalcule seed when no seed provided.
    672673        if seed_num != None:
    673674            seed(seed_num,seed_num)
     675        else:
     676            seed()
     677        if verbose: print "seed:", get_seed()
    674678       
    675679        random_num = randint(0,self_size-1,(int(new_size),))
     
    13931397    alpha_list: the alpha values to test in a single list
    13941398   
    1395     mesh_file: name of the created mesh file
     1399    mesh_file: name of the created mesh file or if passed in will read it.
     1400            NOTE, if there is a mesh file mesh_resolution,
     1401            north_boundary, south... etc will be ignored.
    13961402   
    13971403    mesh_resolution: the maximum area size for a triangle
     
    14041410   
    14051411    USAGE:
    1406         find_optimal_smoothing_parameter(data_file=fileName,
     1412        value, alpha = find_optimal_smoothing_parameter(data_file=fileName,
    14071413                                             alpha_list=[0.0001, 0.01, 1],
    14081414                                             mesh_file=None,
     
    14301436    attribute_smoothed='elevation'
    14311437
    1432     if north_boundary is None or south_boundary is None or \
    1433        east_boundary is None or west_boundary is None:
    1434         no_boundary=True
    1435     else:
    1436         no_boundary=False
    1437    
    1438     if no_boundary is True:
    1439         msg= 'All boundaries must be defined'
    1440         raise msg
    1441 
    14421438    if mesh_file is None:
    14431439        mesh_file='temp.msh'
    14441440
    1445        
    1446        
    1447     poly_topo = [[east_boundary,south_boundary],
    1448                  [east_boundary,north_boundary],
    1449                  [west_boundary,north_boundary],
    1450                  [west_boundary,south_boundary]]
    1451                  
    1452     create_mesh_from_regions(poly_topo,
    1453                              boundary_tags={'back': [2],
    1454                                             'side': [1,3],
    1455                                             'ocean': [0]},
    1456                          maximum_triangle_area=mesh_resolution,
    1457                          filename=mesh_file,
    1458                          use_cache=cache,
    1459                          verbose=verbose)
     1441        if north_boundary is None or south_boundary is None or \
     1442           east_boundary is None or west_boundary is None:
     1443            no_boundary=True
     1444        else:
     1445            no_boundary=False
     1446       
     1447        if no_boundary is True:
     1448            msg= 'All boundaries must be defined'
     1449            raise msg
     1450
     1451        poly_topo = [[east_boundary,south_boundary],
     1452                     [east_boundary,north_boundary],
     1453                     [west_boundary,north_boundary],
     1454                     [west_boundary,south_boundary]]
     1455                     
     1456        create_mesh_from_regions(poly_topo,
     1457                                 boundary_tags={'back': [2],
     1458                                                'side': [1,3],
     1459                                                'ocean': [0]},
     1460                             maximum_triangle_area=mesh_resolution,
     1461                             filename=mesh_file,
     1462                             use_cache=cache,
     1463                             verbose=verbose)
     1464
     1465    else: # if mesh file provided
     1466        #test mesh file exists?
     1467        if access(mesh_file,F_OK) == 0:
     1468            msg="file %s doesn't exist!" %mesh_file
     1469            raise IOError, msg
    14601470
    14611471    #split topo data
     
    15301540        loglog(normal_cov_new[:,0],normal_cov_new[:,1])
    15311541        savefig(plot_name,dpi=300)
    1532    
    1533     remove(mesh_file)
     1542    if mesh_file == 'temp.msh':
     1543        remove(mesh_file)
    15341544   
    15351545    return min(normal_cov_new[:,1]) , normal_cov_new[(argmin(normal_cov_new,axis=0))[1],0]
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r4745 r4816  
    18211821        assert allclose(A,[24, 18, 17, 11, 8])
    18221822 
    1823        
    1824 
    1825 #    def test_find_optimal_smoothing_parameter1(self):
    1826 #        from anuga.pmesh.mesh_interface import create_mesh_from_regions
    1827 #        from anuga.shallow_water import Domain
    1828 #        from anuga.shallow_water import Reflective_boundary
    1829 #        from cmath import cos
    1830 #       
    1831 #        "make sphere domain"
    1832 #
    1833 #        fileName = tempfile.mktemp(".csv")
    1834 #        file = open(fileName,"w")
    1835 #        file.write("x,y,elevation \n")
    1836 #
    1837 #        for i in range(-5,6):
    1838 #            for j in range(-5,6):
    1839 #                z = abs(cos(((i*i) + (j*j))*.1)*2)
    1840 ##                print 'x,y,f',i,j,z
    1841 #                file.write("%s, %s, %s\n" %(i, j, z))
    1842 #               
    1843 #        file.close()
    1844 #
    1845 #
    1846 #
    1847 #        poly = [[5,-5], [5,5], [-5,5], [-5,-5]]
    1848 #        mesh_file='temp.msh'
    1849 #       
    1850 #        create_mesh_from_regions(poly,
    1851 #                        boundary_tags={'back': [2],
    1852 #                                       'side': [1,3],
    1853 #                                       'ocean': [0]},
    1854 #                        maximum_triangle_area=3,
    1855 #                        filename=mesh_file,
    1856 #                        use_cache=False,
    1857 #                        verbose=False)
    1858 #
    1859 #        domain = Domain(mesh_file, use_cache=False, verbose=True)
    1860 #        domain.set_name('test')
    1861 #        domain.set_quantity('elevation', filename = fileName ,verbose=True, alpha=0.1)
    1862 #       
    1863 #
    1864 #        Br = Reflective_boundary(domain)
    1865 #
    1866 #        domain.set_boundary({'back': Br,
    1867 #                         'side': Br,
    1868 #                         'ocean': Br})
    1869 #
    1870 #        for t in domain.evolve(yieldstep = 1, finaltime = 1):
    1871 #            domain.write_time()
    1872 
    18731823    def test_find_optimal_smoothing_parameter(self):
    18741824        """
     
    18801830        from cmath import cos
    18811831       
    1882         fileName = tempfile.mktemp(".csv")
    1883         file = open(fileName,"w")
     1832        filename = tempfile.mktemp(".csv")
     1833        file = open(filename,"w")
    18841834        file.write("x,y,elevation \n")
    18851835
     
    18931843        file.close()
    18941844 
    1895         value, alpha = find_optimal_smoothing_parameter(data_file=fileName,
     1845        value, alpha = find_optimal_smoothing_parameter(data_file=filename,
    18961846                                             alpha_list=[0.0001, 0.01, 1],
    18971847                                             mesh_file=None,
     
    19051855                                             verbose=False)
    19061856
    1907 
    1908         results = Geospatial_data(fileName)
    1909         os.remove(fileName)
    1910        
    1911  #       print value, alpha
     1857        os.remove(filename)
     1858       
     1859#        print value, alpha
    19121860        assert (alpha==0.01)
    19131861
     1862    def test_find_optimal_smoothing_parameter1(self):
     1863        """
     1864        Creates a elevation file represting hill (sort of) and
     1865        Then creates a mesh file and passes the mesh file and the elevation
     1866        file to find_optimal_smoothing_parameter for 3 different alphas,
     1867       
     1868        NOTE the random number seed is provided to control the results
     1869        """
     1870        from cmath import cos
     1871        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     1872       
     1873        filename = tempfile.mktemp(".csv")
     1874        file = open(filename,"w")
     1875        file.write("x,y,elevation \n")
     1876
     1877        for i in range(-5,6):
     1878            for j in range(-5,6):
     1879                #this equation made surface like a circle ripple
     1880                z = abs(cos(((i*i) + (j*j))*.1)*2)
     1881#                print 'x,y,f',i,j,z
     1882                file.write("%s, %s, %s\n" %(i, j, z))
     1883               
     1884        file.close()
     1885        poly=[[5,5],[5,-5],[-5,-5],[-5,5]]
     1886        internal_poly=[[[[1,1],[1,-1],[-1,-1],[-1,1]],.5]]
     1887        mesh_filename= tempfile.mktemp(".msh")
     1888       
     1889        create_mesh_from_regions(poly,
     1890                             boundary_tags={'back': [2],
     1891                                            'side': [1,3],
     1892                                            'ocean': [0]},
     1893                         maximum_triangle_area=3,
     1894                         interior_regions=internal_poly,
     1895                         filename=mesh_filename,
     1896                         use_cache=False,
     1897                         verbose=False)
     1898 
     1899        value, alpha = find_optimal_smoothing_parameter(data_file=filename,
     1900                                             alpha_list=[0.0001, 0.01, 1],
     1901                                             mesh_file=mesh_filename,
     1902                                             plot_name=None,
     1903                                             seed_num=174,
     1904                                             verbose=False)
     1905
     1906        os.remove(filename)
     1907        os.remove(mesh_filename)
     1908       
     1909#        print value, alpha
     1910        assert (alpha==0.01)
     1911
     1912    def test_find_optimal_smoothing_parameter2(self):
     1913        """
     1914        Tests requirement that mesh file must exist or IOError is thrown
     1915       
     1916        NOTE the random number seed is provided to control the results
     1917        """
     1918        from cmath import cos
     1919        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     1920       
     1921        filename = tempfile.mktemp(".csv")
     1922        mesh_filename= tempfile.mktemp(".msh")
     1923       
     1924        try:
     1925            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
     1926                                             alpha_list=[0.0001, 0.01, 1],
     1927                                             mesh_file=mesh_filename,
     1928                                             plot_name=None,
     1929                                             seed_num=174,
     1930                                             verbose=False)
     1931        except IOError:
     1932            pass
     1933        else:
     1934            self.failUnless(0 ==1,  'Error not thrown error!')
    19141935       
    19151936         
     
    19171938
    19181939    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long')
    1919     suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
     1940#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
    19201941#    suite = unittest.makeSuite(Test_Geospatial_data, 'test_split')
    1921 #    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
     1942    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
    19221943    runner = unittest.TextTestRunner() #verbosity=2)
    19231944    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.