source: anuga_work/development/alpha_validation/find_alpha.py @ 5175

Last change on this file since 5175 was 4667, checked in by nick, 18 years ago

removed old data file

File size: 3.6 KB
Line 
1
2from anuga.abstract_2d_finite_volumes import domain
3from anuga.shallow_water import Domain
4from anuga.geospatial_data.geospatial_data import Geospatial_data
5from anuga.coordinate_transforms import Geo_reference
6from anuga.pmesh.mesh_interface import create_mesh_from_regions
7from anuga.shallow_water import Reflective_boundary
8from anuga.utilities.numerical_tools import cov
9from Numeric import array, resize,shape,Float,zeros,take,argsort,argmin
10from pylab import plot, ion, hold,savefig,semilogx,plotting,loglog
11
12
13poly_topo = [[664000,7752000],[664000,7754000],
14             [666000,7754000],[666000,7752000]]
15mesh_dir_name = 'temp.msh'
16                 
17create_mesh_from_regions(poly_topo,
18                         boundary_tags={'back': [2], 'side': [1,3],
19                                        'ocean': [0]},
20                         maximum_triangle_area=1000,
21                         filename=mesh_dir_name,
22                         use_cache=True,
23                         verbose=True)
24
25
26#topo_dir_name = 'pt_hedland_export.txt'
27topo_dir_name = 'pt_hedland_small.txt'
28#sample = 'pt_hedland_small.txt'
29remainder ='remainder.txt'
30#split topo data
31G = Geospatial_data(file_name = topo_dir_name)
32print 'start split'
33
34
35
36
37
38G_sample,G_loss= G.split(0.01, True)
39
40#G_sample.export_points_file(sample)
41
42G_small, G_other = G_sample.split(0.1,True)
43
44#import sys
45#sys.exit()
46
47alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
48#alphas = [0.001,1,100]
49domains = {}
50#domains = []
51
52print 'Setup computational domain'
53for alpha in alphas:
54   
55    #think about making domain absolute when it is created
56    domain = Domain(mesh_dir_name, use_cache=False, verbose=True)
57    print domain.statistics()
58    domain.set_quantity('elevation', 
59                    geospatial_data = G_other,
60                    use_cache = True,
61                    verbose = True,
62                    alpha = alpha)
63    domains[alpha]=domain
64
65points=G_small.get_data_points()
66data=array([],typecode=Float)
67data=resize(data,(len(points),3+len(alphas)))
68#gets relative point from sample
69
70data[:,0]=points[:,0]
71data[:,1]=points[:,1]
72print'geo',domain.geo_reference,points[0:10],points[0:10,0]
73elevation_sample=G_small.get_attributes(attribute_name='elevation')
74
75data[:,2]=elevation_sample
76
77print'data',data[0:10],len(alphas)
78normal_cov=array(zeros([len(alphas),2]),typecode=Float)
79i=0
80
81for domain in domains:
82#    print'domain',domain
83    points_geo=domains[domain].geo_reference.change_points_geo_ref(points)
84   
85
86    print'poin',points_geo[0:10],elevation_sample[0:10]
87    print'quantities',domains[domain].get_quantity_names()
88    elevation_predicted=domains[domain].quantities['elevation'].get_values(interpolation_points=points_geo)
89    print'elevation_predicted',elevation_predicted[0:10],len(points),len(elevation_predicted)
90   
91    data[:,i+3]=elevation_predicted
92    sample_cov= cov(elevation_sample)
93
94    ele_cov= cov(elevation_sample-elevation_predicted)
95    normal_cov[i,:]= [domain,ele_cov/sample_cov]
96
97    print'cov',alpha,'= ',normal_cov, ele_cov, sample_cov
98    i+=1
99   
100   
101print'data',data[100:200]
102print'normal cov',normal_cov
103#to sort array by column
104normal_cov0=normal_cov[:,0]
105normal_cov_new=take(normal_cov,argsort(normal_cov0))
106semilogx(normal_cov_new[:,0],normal_cov_new[:,1])
107#loglog(normal_cov_new[:,0],normal_cov_new[:,1])
108print 'normal_cov_new',normal_cov_new
109savefig("alphas",dpi=300)
110
111print argmin(normal_cov_new,axis=1)
112print min(normal_cov_new[:,1]) , normal_cov_new[(argmin(normal_cov_new,axis=0))[1],0]
113   
Note: See TracBrowser for help on using the repository browser.