Changeset 4642


Ignore:
Timestamp:
Jul 26, 2007, 9:19:04 AM (10 years ago)
Author:
nick
Message:

first beta version of find_optimal_alpha_value, this will find the best alpha value for the data you provide, using cross validation.

comments and tests to come.

File:
1 edited

Legend:

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

    r4641 r4642  
    15431543    return points
    15441544
    1545 #def file2xya(filename):
    1546    
    1547 #    G = Geospatial_data(filename)
    1548 #    G.export_points_file()
     1545
     1546
     1547def find_optimal_smoothing_parameter(data_file,
     1548                                     smoothing_list=None,
     1549                                     attribute_smoothed='elevation',
     1550                                     mesh_file=None,
     1551                                     mesh_resolution=100000,
     1552                                     north_boundary=None,
     1553                                     south_boundary=None,
     1554                                     east_boundary=None,
     1555                                     west_boundary=None,
     1556                                     plot_name=None
     1557                                     ):
     1558    #alpha list, data file, mesh file or boundary area and resolution
     1559    #plot name and attribute smoothed default elevation
     1560
     1561    from anuga.abstract_2d_finite_volumes import domain
     1562    from anuga.shallow_water import Domain
     1563    from anuga.geospatial_data.geospatial_data import Geospatial_data
     1564    from anuga.coordinate_transforms import Geo_reference
     1565    from anuga.pmesh.mesh_interface import create_mesh_from_regions
     1566    from anuga.shallow_water import Reflective_boundary
     1567   
     1568    from anuga.utilities.numerical_tools import cov
     1569    from Numeric import array, resize,shape,Float,zeros,take,argsort
     1570    from pylab import plot, ion, hold,savefig,semilogx,plotting
     1571
     1572    if north_boundary or south_boundary or east_boundary or west_boundary is None:
     1573        no_boundary=True
     1574       
     1575
     1576    if mesh_file is None:
     1577        mesh_file='temp.msh'
     1578       
     1579        if no_boundary is True:
     1580            msg= 'all boundary information or a mesh file'
     1581            raise msg
     1582
     1583    is_inside_polygon
     1584   
     1585    poly_topo = [[east_boundary,south_boundary],[east_boundary,north_boundary],
     1586             [west_boundary,north_boundary],[west_boundary,south_boundary]]
     1587#    mesh_file = 'temp.msh'
     1588                 
     1589    create_mesh_from_regions(poly_topo,
     1590                         boundary_tags={'back': [2], 'side': [1,3],
     1591                                        'ocean': [0]},
     1592                         maximum_triangle_area=mesh_resolution,
     1593                         filename=mesh_file,
     1594                         use_cache=True,
     1595                         verbose=True)
     1596
     1597
     1598    topo_dir_name = 'pt_hedland_export.txt'
     1599#    sample = 'sample.txt'
     1600#    remainder ='remainder.txt'
     1601    #split topo data
     1602    G = Geospatial_data(file_name = topo_dir_name)
     1603    print 'start split'
     1604    G_small, G_other = G.split(0.1,True)
     1605
     1606
     1607
     1608    alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0,1000.0,100000.0]
     1609    #alphas = [0.001,1,100]
     1610    domains = {}
     1611
     1612    print 'Setup computational domain'
     1613    for alpha in alphas:
     1614   
     1615        #think about making domain absolute when it is created
     1616        domain = Domain(mesh_file, use_cache=False, verbose=True)
     1617        print domain.statistics()
     1618        domain.set_quantity(attribute_smoothed,
     1619                    geospatial_data = G_other,
     1620                    use_cache = True,
     1621                    verbose = True,
     1622                    alpha = alpha)
     1623        domains[alpha]=domain
     1624
     1625    points=G_small.get_data_points()
     1626    data=array([],typecode=Float)
     1627    data=resize(data,(len(points),3+len(alphas)))
     1628    #gets relative point from sample
     1629
     1630    data[:,0]=points[:,0]
     1631    data[:,1]=points[:,1]
     1632    print'geo',domain.geo_reference,points[0:10],points[0:10,0]
     1633    elevation_sample=G_small.get_attributes(attribute_name=attribute_smoothed)
     1634
     1635    data[:,2]=elevation_sample
     1636
     1637    print'data',data[0:10],len(alphas)
     1638    normal_cov=array(zeros([len(alphas),2]),typecode=Float)
     1639    i=0
     1640
     1641    for domain in domains:
     1642    #    print'domain',domain
     1643        points_geo=domains[domain].geo_reference.change_points_geo_ref(points)
     1644   
     1645
     1646        print'poin',points_geo[0:10],elevation_sample[0:10]
     1647        print'quantities',domains[domain].get_quantity_names()
     1648        elevation_predicted=domains[domain].quantities[attribute_smoothed].get_values(interpolation_points=points_geo)
     1649        print'elevation_predicted',elevation_predicted[0:10],len(points),len(elevation_predicted)
     1650   
     1651        data[:,i+3]=elevation_predicted
     1652        sample_cov= cov(elevation_sample)
     1653
     1654        ele_cov= cov(elevation_sample-elevation_predicted)
     1655        normal_cov[i,:]= [domain,ele_cov/sample_cov]
     1656
     1657        print'cov',alpha,'= ',normal_cov, ele_cov, sample_cov
     1658        i+=1
     1659   
     1660   
     1661    print'data',data[100:200]
     1662    print'normal cov',normal_cov
     1663    normal_cov0=normal_cov[:,0]
     1664    normal_cov_new=take(normal_cov,argsort(normal_cov0))
     1665    semilogx(normal_cov_new[:,0],normal_cov_new[:,1])
     1666    savefig("alphas",dpi=300)
     1667   
    15491668
    15501669             
Note: See TracChangeset for help on using the changeset viewer.