Index: development/stochastic_study/create_realisations.py
===================================================================
--- development/stochastic_study/create_realisations.py	(revision 3011)
+++ development/stochastic_study/create_realisations.py	(revision 3035)
@@ -72,6 +72,4 @@
                         use_cache = False)
 
-        print V.shape
-        print V[:2,:]
         Z = None
 
Index: development/stochastic_study/plot_spread.py
===================================================================
--- development/stochastic_study/plot_spread.py	(revision 3011)
+++ development/stochastic_study/plot_spread.py	(revision 3035)
@@ -15,14 +15,21 @@
 # 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,5 +41,5 @@
 
 
-for j in range(20):
+for j in range(1):
     #print j, filenames[j]
     plot(time, data[:,j], 'k-')
@@ -45,7 +52,8 @@
           %(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,5 +61,5 @@
 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,6 +70,8 @@
     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)
Index: development/stochastic_study/project.py
===================================================================
--- development/stochastic_study/project.py	(revision 3011)
+++ development/stochastic_study/project.py	(revision 3035)
@@ -27,11 +27,10 @@
 
 # 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
Index: development/stochastic_study/read_realisations.py
===================================================================
--- development/stochastic_study/read_realisations.py	(revision 3011)
+++ development/stochastic_study/read_realisations.py	(revision 3035)
@@ -5,8 +5,12 @@
 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,5 +21,6 @@
                        {'max_realisations': max_realisations,
                         'gauge_number': gauge_number,
-                        'exclude': exclude},
+                        'exclude': exclude,
+                        'sorting': sorting},
                        verbose=verbose)
     else:
@@ -23,10 +28,13 @@
                                     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,5 +60,16 @@
     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,2 +87,41 @@
 
     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
Index: development/stochastic_study/run_model.py
===================================================================
--- development/stochastic_study/run_model.py	(revision 3011)
+++ development/stochastic_study/run_model.py	(revision 3035)
@@ -134,5 +134,5 @@
                                        finaltime = finaltime):
                     pass
-                    #domain.write_time()
+                    domain.write_time()
                     
         
Index: development/stochastic_study/run_simple_model.py
===================================================================
--- development/stochastic_study/run_simple_model.py	(revision 3035)
+++ development/stochastic_study/run_simple_model.py	(revision 3035)
@@ -0,0 +1,175 @@
+"""Stochastic study of the ANUGA implementation of the
+shallow water wave equation.
+
+Very simple flat bed with one spike in the middle where different values are sampled.
+Three sides are reflective and one Dirichlet condition
+
+Output is measured at gauges around the spike.
+
+The system will evolve dynamically towards a steady state.
+
+Suresh Kumar and Ole Nielsen 2006
+"""
+
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+
+# Related major packages
+from Numeric import allclose, zeros, Float
+from RandomArray import normal, seed
+
+# Application specific imports
+from pyvolution.mesh_factory import rectangular
+from pyvolution.shallow_water import Domain
+from pyvolution.shallow_water import Reflective_boundary
+from pyvolution.shallow_water import Dirichlet_boundary
+from pyvolution.util import file_function
+from caching.caching import cache
+import project                 # Definition of file names and polygons
+
+
+#-----------------------------------------------------------------------------
+# Read in processor information
+#-----------------------------------------------------------------------------
+
+try:
+    import pypar
+except:
+    print 'Could not import pypar'
+    myid = 0
+    numprocs = 1
+    processor_name = 'local host'
+else:    
+    myid = pypar.rank()
+    numprocs = pypar.size()
+    processor_name = pypar.Get_processor_name()
+
+print 'I am process %d of %d running on %s' %(myid, numprocs, processor_name)
+
+
+#-----------------------------------------------------------------------------
+# Setup computational domain
+#----------------------------------------------------------------------------- 
+
+# Create basic triangular mesh
+points, vertices, boundary = rectangular(10, 10, 10, 10)
+
+#Create shallow water domain
+domain = Domain(points, vertices, boundary)
+domain.set_name(project.basename)
+domain.set_datadir(project.working_dir) 
+print domain.statistics()
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+
+Br = Reflective_boundary(domain) #Wall
+Bd = Dirichlet_boundary([1, 0, 0])
+
+# Bind boundary objects to tags
+domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('friction', 0.0)
+domain.set_quantity('stage', 0.0)
+
+# Initial conditions
+def f(x,y):
+    z = zeros( x.shape, Float )
+    for i in range(len(x)):
+        if allclose([x[i], y[i]], 5):
+            print i
+            z[i] = 1
+    return z
+
+domain.set_quantity('elevation', f)
+
+
+z = domain.get_quantity('elevation').get_values()
+print z
+#print z != 0
+
+import sys; sys.exit() 
+
+# Run model for different realisations of spike 
+print 'Create %d different realisations' %project.number_of_realisations
+
+# Array of bed elevation values at vertices
+
+#print z
+N = z.shape[0]
+
+timestep = 0.1
+finaltime = 1
+
+seed(13,17)
+for realisation in range(project.number_of_realisations):
+
+    print 'Creating realisation #%d' %realisation
+    #z[N/2-1,:] = normal(project.mean, project.std_dev, 1)
+
+    # Distribute work in round-robin fashion
+    if realisation%numprocs == myid:
+        name = project.basename + '_P%d' %myid    
+        domain.set_name(name)                # Output name
+        #domain.set_quantity('elevation', z)  # Assign bathymetry
+        domain.set_time(0.0)                 # Reset time
+
+        #---------------------------------------------------
+        # Evolve system through time
+        #---------------------------------------------------
+        print 'P%d: Running realisation %d of %d'\
+              %(myid, realisation, project.number_of_realisations)
+
+        t0 = time.time()
+        for t in domain.evolve(yieldstep = timestep,
+                               finaltime = finaltime):
+            pass
+            domain.write_time()
+                    
+        
+        print 'P%d: Simulation of realisation %d took %.2f seconds'\
+              %(myid, realisation, time.time()-t0)
+
+
+
+        
+        #---------------------------------------------------
+        # Now extract the 3 timeseries (Ch 5-7-9) and store them
+        # in three files for this realisation
+
+        print 'P%d: Extracting time series for realisation %d from file %s'\
+              %(myid, realisation, project.working_dir + domain.filename + '.sww')
+        f = file_function(project.working_dir + domain.filename + '.sww',
+                          quantities='stage',
+                          interpolation_points=project.gauges,
+                          verbose=False)
+
+
+        simulation_name = project.working_dir + \
+                          project.basename + '_realisation_%d' %realisation
+
+        print 'P%d: Writing to file %s'\
+              %(myid, simulation_name + '_' + name + '.txt')
+        
+        for k, name in enumerate(project.gauge_names):
+            fid = open(simulation_name + '_' + name + '.txt', 'w')
+            for t in f.get_time():
+                #For all precomputed timesteps
+                val = f(t, point_id = k)[0]
+                fid.write('%f %f\n' %(t, val))
+                
+            fid.close()
+
+            
+pypar.finalize()
Index: development/stochastic_study/test_generator.py
===================================================================
--- development/stochastic_study/test_generator.py	(revision 3035)
+++ development/stochastic_study/test_generator.py	(revision 3035)
@@ -0,0 +1,13 @@
+
+from RandomArray import normal
+
+
+x = normal(50,.1,10000)
+
+# Plot
+from pylab import ion, hold, plot, title, xlabel, ylabel, legend, savefig, show
+
+ion()
+plot(x, 'k.')
+    
+show()
