Changeset 2972

Timestamp:

May 25, 2006, 3:42:31 PM (19 years ago)

Author:

ole

Message:

Added create_bins and refactored pyvolution/mesh.py

Files:

• development/stochastic_study/plot_spread.py (modified) (3 diffs)
• development/stochastic_study/project.py (modified) (2 diffs)
• inundation/pyvolution/mesh.py (modified) (2 diffs)
• inundation/utilities/numerical_tools.py (modified) (2 diffs)

development/stochastic_study/plot_spread.py

@@ -2 +2 @@
 """
 
+# Standard modules
 import sys, os
+
+# Related major packages
+from Numeric import zeros, Float, allclose, arange
 from caching import cache
+
+# Application specific imports
+from utilities.numerical_tools import histogram, create_bins
 import project
 
+# Initialise
 gauge_name = project.gauge_names[0]
-from Numeric import zeros, Float, allclose
-
-
 number_of_realisations = project.number_of_realisations
 
 time = zeros(project.number_of_timesteps, Float)
 data = zeros((project.number_of_timesteps, number_of_realisations), Float)
+
+
+# Read in all realisations-timeseries
 
 j = 0 # Count realisations
@@ -37 +45 @@
 
 ion()
-hold(True)
+#hold(True)
 hold(False)
 
-#for i in range(project.number_of_timesteps):
-#    print i, data[i,:]
-#    plot(data[i,:], 'k.')
-#    raw_input('Next')    
-
+# Simple plot of timeseries for different realisations 
 for j in range(number_of_realisations):
-    print j, data[:,j]
+    #print j, data[:,j]
     plot(data[:,j], 'k-')
     
     raw_input('Next')
+    xlabel('time(s)')
+    ylabel('stage (m)')
+    title('Realisation %d of %d' %(j, number_of_realisations))    
+    
+
+
+# Plot histogram of stage values for each timestep
+for i in range(project.number_of_timesteps):
+    # Plot histogram
+
+    w = data[i,:]
+     
+    bins = create_bins(w, project.number_of_bins)
+    print 'bins', bins
+    hist = histogram(w, bins)
+    print 'hist', hist
+    
+    #plot(w, 'k.')
+    plot(hist, 'k.')
+    xlabel('stage (m)')
+    ylabel('count')
+    title('Timestep %d of %d (t=%.2f)'\
+          %(i, project.number_of_timesteps, time[i]))        
+    raw_input('Next')
+
+
 #title('Gauge %s' %name)
 #xlabel('time(s)')
@@ -55 +85 @@
 #legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
 #savefig(name, dpi = 300)
+
+
+
+
 
 

development/stochastic_study/project.py

@@ -19 +19 @@
 
 # Stats (Suresh ?)
-number_of_realisations = 4
+number_of_realisations = 1
 #std_dev = 0.0026  #Range is 26.035 cm
 std_dev = 0.0013  #Range is 26.035 cm
@@ -25 +25 @@
 blocksize = 100 #How many realisations to fit at a time
 
-
+number_of_bins = 10
 
 

inundation/pyvolution/mesh.py

@@ -689 +689 @@
 
         from Numeric import arange
-        from utilities.numerical_tools import histogram
+        from utilities.numerical_tools import histogram, create_bins
 
         vertex_coordinates = self.vertex_coordinates
@@ -698 +698 @@
 
         #Setup 10 bins for area histogram
-        m = max(areas)
-        bins = arange(0., m, m/10)
+        bins = create_bins(areas, 10)
+        #m = max(areas)
+        #bins = arange(0., m, m/10)
         hist = histogram(areas, bins)
 

inundation/utilities/numerical_tools.py

@@ -8 +8 @@
 #(this should move to somewhere central)
 try:
-    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float    
+    from scipy import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange    
 except:
     #print 'Could not find scipy - using Numeric'
-    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float    
+    from Numeric import ArrayType, array, sum, innerproduct, ravel, sqrt, searchsorted, sort, concatenate, Float, arange    
 
 # Getting an infinite number to use when using Numeric
@@ -212 +212 @@
     return n[1:]-n[:-1]
 
+def create_bins(data, number_of_bins = None):
+    """Safely create bins for use with histogram
+    If data contains only one point or is constant, one bin will be created.
+    If number_of_bins in omitted 10 bins will be created
+    """
+
+    mx = max(data)
+    mn = min(data)
+
+    if mx == mn:
+        bins = array([mn])
+    else:
+        if number_of_bins is None:
+            number_of_bins = 10
+            
+        bins = arange(mn, mx, (mx-mn)/number_of_bins)
+
+    return bins
+