Changeset 5880

Timestamp:

Oct 31, 2008, 1:27:11 PM (16 years ago)

Author:

ole

Message:

Added test for platform problem with read_mux.
This one passes on Linux.

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (4 diffs)
• test_data_manager.py (modified) (6 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4907 +4907 @@
 
     numSrc=len(filenames)
-
+    
     file_params=-1*ones(3,Float) #[nsta,dt,nt]
     
@@ -4919 +4919 @@
         permutation = ensure_numeric([], Float)    
 
+    #print 'filenames', filenames
+    #print 'weights', weights
+                
     # Call underlying C implementation urs2sts_ext.c    
     data = read_mux2(numSrc, filenames, weights, file_params, permutation, verbose)
@@ -4972 +4975 @@
     quantity=zeros((number_of_selected_stations, parameters_index), Float)
     
-    
+
+    #print 'number_of_selected_stations', number_of_selected_stations
     starttime=1e16
     for i in range(number_of_selected_stations):
         quantity[i][:]=data[i][:parameters_index]
+        
+        #print i, parameters_index
+        #print quantity[i][:]
     
         latitudes[i]=data[i][parameters_index]
@@ -5150 +5157 @@
     mux={}
     for i, quantity in enumerate(quantities):
+    
+        #print
+        #print quantity
     
         # For each quantity read the associated list of source mux2 file with 

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -6265 +6265 @@
         weights=ones(1, Float)
         #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
+        times, latitudes, longitudes, elevation, stage, starttime=read_mux2_py([files[0]],weights)
+        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity, starttime_ua=read_mux2_py([files[1]],weights)
         msg='ha and ua have different gauge meta data'
         assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
@@ -6303 +6303 @@
         msg='incorrect gauge va time series returned'
         assert allclose(yvelocity,va)
+        
+        
+    def test_read_mux_platform_problem(self):
+        """test_read_mux_platform_problem
+        
+        This is to test a situation where read_mux returned 
+        wrong values Win32
+        """
+        
+        from Numeric import sin, cos
+        from urs_ext import read_mux2 
+        
+        verbose = False
+                
+        tide = 1.5
+        time_step_count = 10
+        time_step = 0.2
+        
+        times_ref = arange(0, time_step_count*time_step, time_step)
+        #print 'time vector', times_ref
+        
+        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
+        n = len(lat_long_points)
+        
+        
+        # Create different timeseries starting and ending at different times 
+        first_tstep=ones(n,Int)
+        first_tstep[0]+=2   # Point 0 starts at 2
+        first_tstep[1]+=4   # Point 1 starts at 4        
+        first_tstep[2]+=3   # Point 2 starts at 3
+        
+        last_tstep=(time_step_count)*ones(n,Int)
+        last_tstep[0]-=1    # Point 0 ends 1 step early
+        last_tstep[1]-=2    # Point 1 ends 2 steps early                
+        last_tstep[4]-=3    # Point 4 ends 3 steps early        
+        
+        
+        
+        # Create varying elevation data (positive values for seafloor)
+        gauge_depth=20*ones(n,Float)
+        for i in range(n):
+            gauge_depth[i] += i**2
+            
+        #print 'gauge_depth', gauge_depth
+        
+        # Create data to be written to first mux file        
+        ha0=2*ones((n,time_step_count),Float)
+        ha0[0]=arange(0,time_step_count)
+        ha0[1]=arange(time_step_count,2*time_step_count)
+        ha0[2]=arange(2*time_step_count,3*time_step_count)
+        ha0[3]=arange(3*time_step_count,4*time_step_count)
+        ua0=5*ones((n,time_step_count),Float)
+        va0=-10*ones((n,time_step_count),Float)
+
+        # Ensure data used to write mux file to be zero when gauges are
+        # not recording
+        for i in range(n):
+             # For each point
+             
+             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
+                 # For timesteps before and after recording range
+                 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0                                  
+
+
+        
+        # Write first mux file to be combined by urs2sts
+        base_nameI, filesI = self.write_mux2(lat_long_points,
+                                             time_step_count, time_step,
+                                             first_tstep, last_tstep,
+                                             depth=gauge_depth,
+                                             ha=ha0,
+                                             ua=ua0,
+                                             va=va0)
+
+                                             
+                                             
+        # Create data to be written to second mux file        
+        ha1=ones((n,time_step_count),Float)
+        ha1[0]=sin(times_ref)
+        ha1[1]=2*sin(times_ref - 3)
+        ha1[2]=5*sin(4*times_ref)
+        ha1[3]=sin(times_ref)
+        ha1[4]=sin(2*times_ref-0.7)
+                
+        ua1=zeros((n,time_step_count),Float)
+        ua1[0]=3*cos(times_ref)        
+        ua1[1]=2*sin(times_ref-0.7)   
+        ua1[2]=arange(3*time_step_count,4*time_step_count)
+        ua1[4]=2*ones(time_step_count)
+        
+        va1=zeros((n,time_step_count),Float)
+        va1[0]=2*cos(times_ref-0.87)        
+        va1[1]=3*ones(time_step_count)
+        va1[3]=2*sin(times_ref-0.71)        
+        
+        
+        # Ensure data used to write mux file to be zero when gauges are
+        # not recording
+        for i in range(n):
+             # For each point
+             
+             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
+                 # For timesteps before and after recording range
+                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                  
+
+
+                                             
+        # Write second mux file to be combined by urs2sts                                             
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha1,
+                                               ua=ua1,
+                                               va=va1)
+
+                                               
+        # Create ordering file
+        permutation = ensure_numeric([4,0,2])
+
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][1], 
+                                             lat_long_points[index][0]))
+        order_fid.close()
+        
+        
+
+        # -------------------------------------
+        # Now read files back and check values
+        weights = ensure_numeric([1.0])
+
+        # For each quantity read the associated list of source mux2 file with 
+        # extention associated with that quantity
+        file_params=-1*ones(3,Float) #[nsta,dt,nt]
+        OFFSET = 5
+
+        # FILE I
+        for j, file in enumerate(filesI):
+            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
+
+            nsta=int(file_params[0])
+            dt=file_params[1]        
+        
+            number_of_selected_stations = data.shape[0]
+
+            # Index where data ends and parameters begin
+            parameters_index = data.shape[1]-OFFSET          
+                 
+            times=dt*arange(parameters_index)    
+            latitudes=zeros(number_of_selected_stations, Float)
+            longitudes=zeros(number_of_selected_stations, Float)
+            elevation=zeros(number_of_selected_stations, Float)
+            quantity=zeros((number_of_selected_stations, parameters_index), Float)
+            
+            starttime=1e16
+            
+            for i in range(number_of_selected_stations):
+                quantity[i][:]=data[i][:parameters_index]
+        
+                #print i, j, parameters_index
+                #print quantity[i][:]
+
+                if j == 0:
+                    # HA
+                    if i == 0:
+                        assert allclose(quantity[i][:],
+                                        [2., 2., 2., 2., 2., 2., 2., 0., 0., 0.])
+                    if i == 1:
+                        assert allclose(quantity[i][:],                                        
+                                        [0., 0., 2., 3., 4., 5., 6., 7., 8., 0.])
+                    if i == 2:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 0., 23., 24., 25., 26., 27., 28., 29.])
+                                        
+                if j == 1:
+                    # UA
+                    if i == 0:
+                        assert allclose(quantity[i][:],
+                                        [5., 5., 5., 5., 5., 5., 5., 0., 0., 0.])
+                    if i == 1:
+                        assert allclose(quantity[i][:],                                        
+                                        [0., 0., 5., 5., 5., 5., 5., 5., 5., 0.])
+                    if i == 2:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 0., 5., 5., 5., 5., 5., 5., 5.])
+                if j == 2:
+                    # VA
+                    if i == 0:
+                        assert allclose(quantity[i][:],                                        
+                                        [-10., -10., -10., -10., -10., -10., -10., 0., 0., 0.])
+                    if i == 1:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., -10., -10., -10., -10., -10., -10., -10., 0.])
+                    if i == 2:
+                        assert allclose(quantity[i][:],                        
+                                        [0., 0., 0., -10., -10., -10., -10., -10., -10., -10.])
+                                                                                
+        
+        # FILE II
+        for j, file in enumerate(filesII):
+            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
+
+            nsta=int(file_params[0])
+            dt=file_params[1]        
+        
+            number_of_selected_stations = data.shape[0]
+
+            # Index where data ends and parameters begin
+            parameters_index = data.shape[1]-OFFSET          
+                 
+            times=dt*arange(parameters_index)    
+            latitudes=zeros(number_of_selected_stations, Float)
+            longitudes=zeros(number_of_selected_stations, Float)
+            elevation=zeros(number_of_selected_stations, Float)
+            quantity=zeros((number_of_selected_stations, parameters_index), Float)
+            
+            starttime=1e16
+            
+            for i in range(number_of_selected_stations):
+                quantity[i][:]=data[i][:parameters_index]
+        
+                #print i, parameters_index
+                #print quantity[i][:]
+                
+                
+                if j == 0:
+                    # HA
+                    if i == 0:
+                        assert allclose(quantity[i][:],
+                                        [-0.64421767, -0.29552022,  0.09983341,  0.47942555,
+                                          0.78332692,  0.9635582, 0.99166483,  0., 0., 0.])
+                    if i == 1:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 0.38941833, 0.56464249, 0.71735609, 0.84147096,
+                                         0.93203908, 0.98544973, 0.99957359, 0.])
+                    if i == 2:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 0., 3.377316, -0.29187071, -3.78401256,
+                                         -4.98082304, -3.15633321, 0.58274603, 3.9683392])                        
+
+                                        
+                if j == 1:
+                    # UA
+                    if i == 0:
+                        assert allclose(quantity[i][:],
+                                        [2., 2., 2., 2., 2., 2., 2., 0., 0., 0.])                        
+                    if i == 1:
+                        assert allclose(quantity[i][:],                                        
+                                        [0., 0., 2.76318288, 2.47600675, 2.09012008, 1.62090695,
+                                         1.08707321, 0.5099014, -0.08759857, 0.])
+                    if i == 2:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 0., 33., 34., 35., 36., 37., 38., 39.])                        
+
+                if j == 2:
+                    # VA
+                    if i == 0:
+                        assert allclose(quantity[i][:],                                       
+                                        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])                         
+                    if i == 1:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 1.78313661, 1.92754185, 1.99510205, 1.98312378,
+                                         1.89208472, 1.72561419, 1.49034882, 0.])
+                    if i == 2:
+                        assert allclose(quantity[i][:],
+                                        [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
+            
+        
+        
 
     def test_urs2sts0(self):
@@ -7447 +7720 @@
         sts_file = base_nameI + '.sts'
         fid = NetCDFFile(sts_file)
+        
 
         # Check that original indices have been stored
@@ -7504 +7778 @@
         va_ref = take(va0, permutation)                
 
-        gauge_depth_ref = take(gauge_depth, permutation)                         
-
-
-        #print 
-        #print stage0
-        #print transpose(ha_ref)+tide - stage0
-
+        gauge_depth_ref = take(gauge_depth, permutation)                      
         
         assert allclose(transpose(ha_ref)+tide, stage0)  # Meters
+        
+        
         
         #Check the momentums - ua
@@ -7561 +7831 @@
         sts_file = base_nameI + '.sts'
         fid = NetCDFFile(sts_file)
+        
         
         # Check that original indices have been stored
@@ -10661 +10932 @@
     #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_combined_sources')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts')       
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_ordering_different_sources')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_prob')