Changeset 3035

Timestamp:

May 31, 2006, 5:03:59 PM (19 years ago)

Author:

ole

Message:

Working with Suresh at ACFR

Location:

development/stochastic_study

Files:

• create_realisations.py (modified) (1 diff)
• plot_spread.py (modified) (5 diffs)
• project.py (modified) (1 diff)
• read_realisations.py (modified) (5 diffs)
• run_model.py (modified) (1 diff)
• run_simple_model.py (added)
• test_generator.py (added)

development/stochastic_study/create_realisations.py

@@ -72 +72 @@
                         use_cache = False)
 
-        print V.shape
-        print V[:2,:]
         Z = None
 

development/stochastic_study/plot_spread.py

@@ -15 +15 @@
 # Initialise
 # Read in all realisations-timeseries
-study = 'cyclone1'
-study = 'cyclone3'
-#study = 'nautilus3'
+
+#study = 'nautilus1' #~70 realisations, blocks of 100, sequential, stddev= 0.0006
+#study = 'nautilus3' #~100 realisations, blocks of 10, sequential, stddev= 0.0013 
+study = 'cyclone1' #~4000 realisations, blocks of 100, parallel, stddev= 0.0006
+#study = 'cyclone3' #~120 realisations, blocks of 10, parallel, stddev= 0.0013
 
 time, data, filenames = read_realisations(study, 
                                           #max_realisations = 200,
-                                          gauge_number=2,
-                                          use_cache=False)
+                                          gauge_number=0,
+                                          sorting='randomised',
+                                          #sorting='numerical',
+                                          verbose=True,
+                                          use_cache=True)
 number_of_realisations = data.shape[1]
+number_of_timesteps = data.shape[0]
 print 'Read %d realisations' %number_of_realisations
+
 
 # Plot
@@ -34 +41 @@
 
 
-for j in range(20):
+for j in range(1):
     #print j, filenames[j]
     plot(time, data[:,j], 'k-')
@@ -45 +52 @@
           %(study, shortname, j, number_of_realisations-1))
     
-    #title('Study %s: timeseries for %s (realisation %d of %d)'\
-    #      %(study, filenames[j], j, number_of_realisations-1))    
-    raw_input('Next')
+    #raw_input('Next')
+
+
+raw_input('Stats')    
  
 
@@ -53 +61 @@
 hold(False)
 # Plot spread of stage values for each timestep
-for i in range(300,320): #project.number_of_timesteps):
+for i in range(0,project.number_of_timesteps):
     # Plot histogram#
 
@@ -62 +70 @@
     title('Study %s: spread at timestep %d of %d (t=%.2f)'\
           %(study, i, project.number_of_timesteps-1, time[i]))        
-    raw_input('Next')
+raw_input('Next')
 
+show()
+import sys; sys.exit()
 
 hold(False)

development/stochastic_study/project.py

@@ -27 +27 @@
 
 # Stats (Suresh ?)
-number_of_realisations = 1
+number_of_realisations = 2
 
 #number_of_realisations = 16000
 
-#std_dev = 0.0026  #Range is 26.035 cm
-std_dev = 0.0013  #Range is 26.035 cm   (simulation 3) 
-#std_dev = 0.0006  #Range is 26.035 cm  (simulation 1)
+std_dev = 0.0013 #m  #Range is 26.035 cm   (simulation 3) 
+#std_dev = 0.0006  #m  #Range is 26.035 cm  (simulation 1)
 mean = 0.0
 blocksize = 100 #How many realisations to fit at a time

development/stochastic_study/read_realisations.py

@@ -5 +5 @@
 import project
 from Numeric import Float, zeros
+from RandomArray import randint
 from caching import cache
 
 def read_realisations(subdir, max_realisations = None, gauge_number=0,
-                      exclude=None, use_cache=False, verbose=True):
+                      exclude=None,
+                      sorting=None,
+                      use_cache=False, 
+                      verbose=True):
     """Read realisations
     """
@@ -17 +21 @@
                        {'max_realisations': max_realisations,
                         'gauge_number': gauge_number,
-                        'exclude': exclude},
+                        'exclude': exclude,
+                        'sorting': sorting},
                        verbose=verbose)
     else:
@@ -23 +28 @@
                                     max_realisations = max_realisations,
                                     gauge_number=gauge_number,
-                                    exclude=exclude)
+                                    exclude=exclude,
+                                    sorting=sorting)
 
     return result
 
+
 def _read_realisations(subdir, max_realisations = None, gauge_number=0,
-                       exclude=None, use_cache=False, verbose=True):
+                       exclude=None, sorting=None,
+                       use_cache=False, verbose=True):
     """Read realisations
     """
@@ -52 +60 @@
     time = zeros(project.number_of_timesteps, Float)
     data = zeros((project.number_of_timesteps, len(filenames)), Float)
-    
+
+    if sorting is None:
+        pass
+    elif sorting == 'numerical':
+        filenames = sort_numerically(filenames)
+    elif sorting == 'randomised':    
+        filenames = randomise(filenames)
+    else:
+        print 'Invalid value for sorting: %s' %sorting
+        import sys; sys.exit() 
+
+
     # Read data        
 
@@ -68 +87 @@
 
     return time, data, filenames
+
+
+def sort_numerically(filenames):
+    """Sort filenames numerically
+    simulation_realisation_n_ch.... (sort by n)
+    """
+
+    print 'Sorting'
+    sorted_filenames = []
+
+    numbers = []
+    for filename in filenames:
+        fields = filename.split('_')
+        n = int(fields[2])
+        numbers.append(n)
+
+    A = zip(numbers, filenames)
+    A.sort()
+    
+    sorted_filenames = [a[1] for a in A]    
+    
+    return sorted_filenames
+
+
+def randomise(filenames):
+    """Randomise filenames
+    """
+
+    print 'Randomising'
+    randomised_filenames = []
+
+
+    while len(filenames) > 0:
+        n = randint(0, len(filenames))
+        filename = filenames[n]
+        del filenames[n]
+        randomised_filenames.append(filename)
+     
+    return randomised_filenames

development/stochastic_study/run_model.py

@@ -134 +134 @@
                                        finaltime = finaltime):
                     pass
-                    #domain.write_time()
+                    domain.write_time()