diff --git a/VERSION b/VERSION
index 6c6aa7c..a3df0a6 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-0.1.0
\ No newline at end of file
+0.8.0
diff --git a/molnsutil/molnsutil.py b/molnsutil/molnsutil.py
index 5fed0d3..621c153 100755
--- a/molnsutil/molnsutil.py
+++ b/molnsutil/molnsutil.py
@@ -31,18 +31,21 @@ def run(seed, number_of_trajectories):
 from os import environ
 import logging
 from boto.s3.connection import S3Connection
-logging.basicConfig(filename="boto.log", level=logging.DEBUG)
+#logging.basicConfig(filename="boto.log", level=logging.DEBUG)
 from boto.s3.key import Key
 import uuid
 import math
 import molns_cloudpickle as cloudpickle
 import random
 import copy
+import inspect
 
 import swiftclient.client
 import IPython.parallel
 import uuid
 from IPython.display import HTML, Javascript, display
+import os
+import sys
 
 import itertools
 
@@ -68,7 +71,8 @@ class MolnsUtilStorageException(Exception):
 #   s3.json needs to be created and put in .molns/s3.json in the root of the home directory.
 
 import os
-def get_s3config():
+def get_persisistent_storage_config():
+    """ Return the configuration for the persistent storage. """
     try:
         with open(os.environ['HOME']+'/.molns/s3.json','r') as fh:
             s3config = json.loads(fh.read())
@@ -82,6 +86,7 @@ def get_s3config():
 class LocalStorage():
     """ This class provides an abstraction for storing and reading objects on/from
         the ephemeral storage. """
+    
 
     def __init__(self, folder_name="/home/ubuntu/localarea"):
         self.folder_name = folder_name
@@ -127,7 +132,7 @@ def delete(self,filename):
 
 class S3Provider():
     def __init__(self, bucket_name):
-        s3config = get_s3config()
+        s3config = get_persisistent_storage_config()
         self.connection = S3Connection(is_secure=False,
                                  calling_format=boto.s3.connection.OrdinaryCallingFormat(),
                                  **s3config['credentials']
@@ -189,7 +194,7 @@ def list(self):
 
 class SwiftProvider():
     def __init__(self, bucket_name):
-        s3config = get_s3config()
+        s3config = get_persisistent_storage_config()
         self.connection = swiftclient.client.Connection(auth_version=2.0,**s3config['credentials'])
         self.set_bucket(bucket_name)
 
@@ -253,7 +258,7 @@ class PersistentStorage():
     """
 
     def __init__(self, bucket_name=None):
-        s3config = get_s3config()
+        s3config = get_persisistent_storage_config()
         if bucket_name is None:
             # try reading it from the config file
             try:
@@ -397,11 +402,13 @@ def run_ensemble_map_and_aggregate(model_class, parameters, param_set_id, seed_b
     """ Generate an ensemble, then run the mappers are aggreator.  This will not store the results. """
     import sys
     import uuid
+    import molnsutil.molns_cloudpickle as cp
+
     if aggregator is None:
         aggregator = builtin_aggregator_list_append
     # Create the model
     try:
-        model_class_cls = cloudpickle.loads(model_class)
+        model_class_cls = cp.loads(model_class)
         if parameters is not None:
             model = model_class_cls(**parameters)
         else:
@@ -459,6 +466,8 @@ def run_ensemble(model_class, parameters, param_set_id, seed_base, number_of_tra
     import sys
     import uuid
     from molnsutil import PersistentStorage, LocalStorage, SharedStorage
+    import molnsutil.molns_cloudpickle as cp
+
 
     if storage_mode=="Shared":
         storage  = SharedStorage()
@@ -468,7 +477,7 @@ def run_ensemble(model_class, parameters, param_set_id, seed_base, number_of_tra
         raise MolnsUtilException("Unknown storage type '{0}'".format(storage_mode))
     # Create the model
     try:
-        model_class_cls = cloudpickle.loads(model_class)
+        model_class_cls = cp.loads(model_class)
         if parameters is not None:
             model = model_class_cls(**parameters)
         else:
@@ -561,34 +570,76 @@ def map_and_aggregate(results, param_set_id, mapper, aggregator=None, cache_resu
 
     #return res
 
+############################################################################
 class DistributedEnsemble():
     """ A class to provide an API for execution of a distributed ensemble. """
 
-    def __init__(self, model_class=None, parameters=None, client=None, num_engines=None):
+    my_class_name = 'DistributedEnsemble'
+
+    #-----------------------------------------------------------------------------------
+    @classmethod
+    def delete(cls, name):
+        """ Static method to remove the state of a distributed comptutation from the system."""
+        # delete realization
+        try:
+            with open('.molnsutil/{1}-{0}'.format(name, cls.my_class_name)) as fd:
+                state = pickle.load(fd)
+                
+                if state['storage_mode'] is not None:
+                    if state['storage_mode'] == "Shared":
+                        ss = SharedStorage()
+                    elif state['storage_mode'] == "Persistent":
+                        ss = PersistentStorage()
+                    for param_set_id in state['result_list']:
+                        for filename in state['result_list'][param_set_id]:
+                            try:
+                                ss.delete(filename)
+                            except OSError as e:
+                                pass
+            # delete database file
+            os.remove('.molnsutil/{1}-{0}'.format(name, cls.my_class_name))
+        except Exception as e:
+            sys.stderr.write('delete(): {0}'.format(e))
+
+    #-----------------------------------------------------------------------------------
+    def __init__(self, name=None, model_class=None, parameters=None, client=None, num_engines=None, ignore_model_mismatch=False):
         """ Constructor """
-        self.my_class_name = 'DistributedEnsemble'
+        if not isinstance(name, str):
+            raise MolnsUtilException("name not specified")
+        self.name = name
+        if not inspect.isclass(model_class):
+            raise MolnsUtilException("model_class not a class")
         self.model_class = cloudpickle.dumps(model_class)
-        self.parameters = [parameters]
-        self.number_of_trajectories = 0
-        self.seed_base = self.generate_seed_base()
-        self.storage_mode = None
-        # A chunk list
-        self.result_list = {}
         # Set the Ipython.parallel client
         self.num_engines = num_engines
         self._update_client(client)
-
+        
+        if not self.load_state(ignore_model_mismatch=ignore_model_mismatch):
+            # Set defaults if state not found
+            self.parameters = [parameters]
+            self.number_of_trajectories = 0
+            self.seed_base = self.generate_seed_base()
+            self.storage_mode = None
+            self.result_list = {}
+            self.running_MapReduceTask = None
+            self.running_SimulationTask = None
+            self.reduced_results = None
+            self.mapped_results = None
+            self.number_of_results = 0
+            self.step1_complete = False
+            self.step2_complete = False
+            self.step3_complete = False
+            self.mapper_fn = None
+            self.aggregator_fn  = None
+            self.reducer_fn = None
+
+    #-----------------------------------------------------------------------------------
     def generate_seed_base(self):
         """ Create a random number and truncate to 64 bits. """
-        x = int(uuid.uuid4())
-        if x.bit_length() >= 64:
-            x = x & ((1<<64)-1)
-            if x > (1 << 63) -1:
-                x -= 1 << 64
-        return x
+        return abs(int(random.getrandbits(31)))
 
-    #--------------------------
-    def save_state(self, name):
+    #-----------------------------------------------------------------------------------
+    def save_state(self):
         """ Serialize the state of the ensemble, for persistence beyond memory."""
         state = {}
         state['model_class'] = self.model_class
@@ -597,99 +648,238 @@ def save_state(self, name):
         state['seed_base'] = self.seed_base
         state['result_list'] = self.result_list
         state['storage_mode'] = self.storage_mode
+        #
+        state['reduced_results'] = self.reduced_results
+        state['mapped_results'] = self.mapped_results
+        state['number_of_results'] = self.number_of_results
+        state['step1_complete'] = self.step1_complete
+        state['step2_complete'] = self.step2_complete
+        state['step3_complete'] = self.step3_complete
+        state['mapper_fn'] = self.mapper_fn
+        state['aggregator_fn'] = self.aggregator_fn
+        state['reducer_fn'] = self.reducer_fn
+        #
+        if self.running_MapReduceTask is None:
+            state['running_MapReduceTask'] = None
+        else:
+            state['running_MapReduceTask'] = self.running_MapReduceTask.msg_ids
+        if self.running_SimulationTask is None:
+            state['running_SimulationTask'] = None
+        else:
+            state['running_SimulationTask'] = self.running_SimulationTask.msg_ids
         if not os.path.isdir('.molnsutil'):
             os.makedirs('.molnsutil')
-        with open('.molnsutil/{1}-{0}'.format(name, self.my_class_name)) as fd:
+        with open('.molnsutil/{1}-{0}'.format(self.name, self.my_class_name),'w+') as fd:
             pickle.dump(state, fd)
 
-    def load_state(self, name):
+    #-----------------------------------------------------------------------------------
+    def load_state(self, ignore_model_mismatch=False):
         """ Recover the state of an ensemble from a previous save. """
-        with open('.molnsutil/{1}-{0}'.format(name, self.my_class_name)) as fd:
-            state = pickle.load(fd)
-        if state['model_class'] is not self.model_class:
-            raise MolnsUtilException("Can only load state of a class that is identical to the original class")
-        self.parameters = state['parameters']
-        self.number_of_trajectories = state['number_of_trajectories']
-        self.seed_base = state['seed_base']
-        self.result_list = state['result_list']
-        self.storage_mode = state['storage_mode']
+        try:
+            with open('.molnsutil/{1}-{0}'.format(self.name, self.my_class_name)) as fd:
+                state = pickle.load(fd)
+            if state['model_class'] != self.model_class and not ignore_model_mismatch:
+                #sys.stderr.write("Error loading saved state\n\n")
+                #sys.stderr.write("state['model_class']={0}\n\n".format(state['model_class']))
+                #sys.stderr.write("self.model_class={0}\n\n".format(self.model_class))
+                #sys.stderr.write("state['model_class'] != self.model_class {0}\n\n".format(state['model_class'] != self.model_class))
+                #TODO: Minor differences show up in the pickled string, but the classes are still identical.  Find a way around this.
+                raise MolnsUtilException("Can only load state of a class that is identical to the original class. Use '{0}.delete(name=\"{1}\")' to remove previous state.  Use the argument 'ignore_model_mismatch=True' to override".format(self.my_class_name, self.name))
+                #TODO: Check to be sure the state is sane.  Both tasks can not be running, result list and number_of_trajectories should match up, (others?).
+            
+            self.parameters = state['parameters']
+            self.number_of_trajectories = state['number_of_trajectories']
+            self.seed_base = state['seed_base']
+            self.result_list = state['result_list']
+            self.storage_mode = state['storage_mode']
+            #
+            self.reduced_results = state['reduced_results']
+            self.mapped_results = state['mapped_results']
+            self.number_of_results = state['number_of_results']
+            self.step1_complete = state['step1_complete']
+            self.step2_complete = state['step2_complete']
+            self.step3_complete = state['step3_complete']
+            self.mapper_fn = state['mapper_fn']
+            self.aggregator_fn  = state['aggregator_fn']
+            self.reducer_fn = state['reducer_fn']
+            #
+            if state['running_MapReduceTask'] is None:
+                self.running_MapReduceTask = None
+            else:
+                self.running_MapReduceTask = self.c.get_result(state['running_MapReduceTask'])
+            if state['running_SimulationTask'] is None:
+                self.running_SimulationTask = None
+            else:
+                self.running_SimulationTask = self.c.get_result(state['running_SimulationTask'])
+            return True
+        except IOError as e:
+            return False
 
-    #--------------------------
+    #-----------------------------------------------------------------------------------
     # MAIN FUNCTION
-    #--------------------------
-    def run(self, mapper, aggregator=None, reducer=None, number_of_trajectories=None, chunk_size=None, verbose=True, progress_bar=True, store_realizations=True, storage_mode="Shared", cache_results=False):
-        """ Main entry point """
-        if store_realizations:
-            if self.storage_mode is None:
-                if storage_mode != "Persistent" and storage_mode != "Shared":
-                    raise MolnsUtilException("Acceptable values for 'storage_mode' are 'Persistent' or 'Shared'")
-                self.storage_mode = storage_mode
-            elif self.storage_mode != storage_mode:
-                raise MolnsUtilException("Storage mode already set to {0}, can not mix storage modes".format(self.storage_mode))
-            # Do we have enough trajectores yet?
-            if number_of_trajectories is None and self.number_of_trajectories == 0:
-                raise MolnsUtilException("number_of_trajectories is zero")
-            # Run simulations
-            if self.number_of_trajectories < number_of_trajectories:
-                self.add_realizations( number_of_trajectories - self.number_of_trajectories, chunk_size=chunk_size, verbose=verbose, storage_mode=storage_mode)
+    #-----------------------------------------------------------------------------------
+    def run(self, mapper=None, aggregator=None, reducer=None, number_of_trajectories=None, chunk_size=None, verbose=True, progress_bar=True, store_realizations=True, storage_mode="Shared", cache_results=False):
+        """ Main entry point for executing parallel MOLNs computations.
+        
+        Arguments:
+            mapper                      [required] Python function that takes as input the simulation result. This function is applied to each simulation trajectory.
+            aggregator                  [optional] Python function that aggregates the output of the mapper function on each worker engine.
+            reducer                     [optional] Python function that takes as input the output of all the aggregator functions on each worker.  One insteance of this function is run for each parameter point.
+            number_of_trajectories      [required] Integer, number of simulation trajectories to execute for each parameter point.
+            chunk_size                  [optional] Integer, group a number of trajectories into a block for efficicnet execution. One aggregator will be run for each chunk.
+            verbose                     [optional, default True] Print the status of the computation.
+            progress_bar                [optional, default True] Display a javascript progress bar to indicate the progress of the computation.
+            store_realizations          [optional, default True] If set to False, the intermediary results will be deleted as soon as the computation is complete.
+            storage_mode                [required, default 'Shared'] Either 'Shared' or 'Persistent'. Store the intermediary results in the ephemeral shared filesystem ('Shared'), or the cloud object store ('Persistent' e.g. Amazon S3 or OpenStack Swift).
+            cache_results               [optional, experimental, default False] Store the intermediary results in the ephemeral storage on each compute node.
+        Returns:
+            The output of the computation is returned.  For 'DistributedEnsemble' class, this is the output of the 'reducer' function.  For the 'ParameterSweep' class, this is a 'ParameterSweepResultList' object.
+        Raises:
+            MolnsUtilException on error.
+        """
+        #####
+        # 0. Validate input.
+        if mapper is None or not hasattr(mapper, '__call__'):
+            raise MolnsUtilException("mapper function not specified")
+        if self.storage_mode is None:
+            if storage_mode != "Persistent" and storage_mode != "Shared":
+                raise MolnsUtilException("Acceptable values for 'storage_mode' are 'Persistent' or 'Shared'")
+            self.storage_mode = storage_mode
+        elif self.storage_mode != storage_mode:
+            raise MolnsUtilException("Storage mode already set to {0}, can not mix storage modes".format(self.storage_mode))
+        if number_of_trajectories is None or number_of_trajectories == 0:
+            raise MolnsUtilException("'number_of_trajectories' is zero.")
+        elif not store_realizations and self.number_of_trajectories > 0 and self.number_of_trajectories != number_of_trajectories:
+            raise MolnsUtilException("'number_of_trajectories' changed since first call.  Value can only be changed if store_realizations is True")
+            #TODO: Fix, if store_realizations is false, number_of_trajectories=0 after success
+        elif store_realizations and self.number_of_trajectories > 0 and self.number_of_trajectories > number_of_trajectories:
+            raise MolnsUtilException("'number_of_trajectories' less than first call. Can not reduce number of realizations.")
+        # Check if the mapper, aggregator, and reducer functions are the same as previously run.  If not throw error.  Use 'clear_results()' to reset
+        mapper_fn_pkl = cloudpickle.dumps(mapper)
+        if self.mapper_fn is not None and self.mapper_fn != mapper_fn_pkl:
+            raise MolnsUtilException("'mapper' function has changed since results have been computed.  Use 'clear_results()' to reset.")
+        else:
+            self.mapper_fn = mapper_fn_pkl
+        aggregator_fn_pkl = cloudpickle.dumps(aggregator)
+        if self.aggregator_fn is not None and self.aggregator_fn != aggregator_fn_pkl:
+            raise MolnsUtilException("'aggregator' function has changed since results have been computed.  Use 'clear_results()' to reset.")
+        else:
+            self.aggregator_fn = aggregator_fn_pkl
+        reducer_fn_pkl = cloudpickle.dumps(reducer)
+        if self.reducer_fn is not None and self.reducer_fn != reducer_fn_pkl:
+            raise MolnsUtilException("'reducer' function has changed since results have been computed.  Use 'clear_results()' to reset.")
+        else:
+            self.reducer_fn = reducer_fn_pkl
+        if self.number_of_trajectories > 0 and not store_realizations and self.number_of_trajectories != number_of_trajectories:
+            raise MolnsUtilException("'number_of_trajectories' is not the same as the original call. It can only be changed if 'store_realizations' is True.")
 
-            if chunk_size is None:
-                chunk_size = self._determine_chunk_size(self.number_of_trajectories)
+        #####
+        # Shortcut, check if computation is complete
+        if self.step3_complete:
             if verbose:
-                print "Running mapper & aggregator on the result objects (number of results={0}, chunk size={1})".format(self.number_of_trajectories*len(self.parameters), chunk_size)
-            else:
-                progress_bar=False
+                print "Using previously computed results."
+            return self.reduced_results
 
-            # chunks per parameter
-            num_chunks = int(math.ceil(self.number_of_trajectories/float(chunk_size)))
-            chunks = [chunk_size]*(num_chunks-1)
-            chunks.append(self.number_of_trajectories-chunk_size*(num_chunks-1))
-            # total chunks
-            pchunks = chunks*len(self.parameters)
-            num_pchunks = num_chunks*len(self.parameters)
-            pparams = []
-            param_set_ids = []
-            presult_list = []
-            for id, param in enumerate(self.parameters):
-                param_set_ids.extend( [id]*num_chunks )
-                pparams.extend( [param]*num_chunks )
-                for i in range(num_chunks):
-                    presult_list.append( self.result_list[id][i*chunk_size:(i+1)*chunk_size] )
 
-            results = self.lv.map_async(map_and_aggregate, presult_list, param_set_ids, [mapper]*num_pchunks,[aggregator]*num_pchunks,[cache_results]*num_pchunks)
+        ######
+        # 1. Run simulations
+        #sys.stderr.write("[1] self.number_of_trajectories < number_of_trajectories: {0} {1} {2}\n".format(self.number_of_trajectories,number_of_trajectories,self.number_of_trajectories < number_of_trajectories))
+        #sys.stderr.write("[1] (not self.step1_complete): {0} {1}\n".format(self.step1_complete, (not store_realizations)))
+        if (not self.step1_complete) or (store_realizations == True and self.number_of_trajectories < number_of_trajectories):
+            if verbose:
+                print "Step 1: Computing simulation trajectories."
+            self.add_realizations( number_of_trajectories - self.number_of_trajectories, chunk_size=chunk_size, verbose=verbose, storage_mode=storage_mode, progress_bar=progress_bar)
+
+            self.step1_complete=True
         else:
-            # If we don't store the realizations (or use the stored ones)
-            if chunk_size is None:
-                chunk_size = self._determine_chunk_size(number_of_trajectories)
-            if not verbose:
-                progress_bar=False
-            else:
-                print "Generating {0} realizations of the model, running mapper & aggregator (chunk size={1})".format(number_of_trajectories,chunk_size)
+            if verbose:
+                print "Step 1: Done. Using previously computed simulation trajectories."
+
+        ######
+        # 2. Run Map function for the MapReduce post-processing
+        if (not self.step2_complete) or (store_realizations == True and self.number_of_results < self.number_of_trajectories):
+            if verbose:
+                print "Step 2: Running mapper & aggregator on the result objects (number of results={0}, chunk size={1})".format(self.number_of_trajectories*len(self.parameters), chunk_size)
+            self.run_mappers(mapper, aggregator, chunk_size=chunk_size, progress_bar=progress_bar, cache_results=cache_results)
+            self.step2_complete=True
+            self.step3_complete=False # To ensure that step 3 is always run after step 2 runs (as in the case of a re-execution).
+        else:
+            if verbose:
+                print "Step 2: Done. Using previously computed results."
 
+        ######
+        # 3. Run Reduce function for the MapReduce post-processing
+        if not self.step3_complete:
+            if verbose:
+                print "Step 3: Running reducer on mapped and aggregated results (size={0})".format(len(self.mapped_results[0]))
+            if reducer is None:
+                reducer = builtin_reducer_default
+            self.reduced_results = self.run_reducer(reducer)
+            self.step3_complete = True
+            self.save_state()
+        else:
+            if verbose:
+                print "Step 3: Done. Using previously computed results."
+
+        ######
+        # Clean up
+        if not store_realizations:
+            self.delete_realizations()
+            self.delete_results()
+            self.save_state()
+
+        ######
+        # Return results
+        return self.reduced_results
+
+    #-----------------------------------------------------------------------------------
+    def run_reducer(self, reducer):
+        """ Inside the run() function, apply the reducer to all of the map'ped-aggregated result values. """
+        return reducer(self.mapped_results[0], parameters=self.parameters[0])
+
+    #-----------------------------------------------------------------------------------
+    def run_mappers(self, mapper, aggregator, chunk_size=None, progress_bar=True, cache_results=False):
+        """ Run the mapper and aggrregator function on each of the simulation trajectories. """
+        num_results_to_compute = self.number_of_trajectories - self.number_of_results
+        if self.running_MapReduceTask is None:
+            #sys.stderr.write('run_mappers(): Starting MapReduce Task\n')
+            # If number_of_trajectories > 0 and number_of_results < number_of_trajectories, only run mappers on the 'new' trajectories.
+            offset = self.number_of_results
             # chunks per parameter
-            num_chunks = int(math.ceil(number_of_trajectories/float(chunk_size)))
+            if chunk_size is None:
+                chunk_size = self._determine_chunk_size(num_results_to_compute)
+            #sys.stderr.write('chunk_size={0}, num_results_to_compute={1}\n'.format(chunk_size, num_results_to_compute))
+            if chunk_size < 1: chunk_size = 1
+            num_chunks = int(math.ceil(num_results_to_compute/float(chunk_size)))
             chunks = [chunk_size]*(num_chunks-1)
-            chunks.append(number_of_trajectories-chunk_size*(num_chunks-1))
+            chunks.append(num_results_to_compute-chunk_size*(num_chunks-1))
             # total chunks
             pchunks = chunks*len(self.parameters)
             num_pchunks = num_chunks*len(self.parameters)
-            pparams = []
+            #pparams = []
             param_set_ids = []
-            for id, param in enumerate(self.parameters):
-                param_set_ids.extend( [id]*num_chunks )
-                pparams.extend( [param]*num_chunks )
+            presult_list = []
+            try:  #'try' is for Debugging
+                for id, param in enumerate(self.parameters):
+                    param_set_ids.extend( [id]*num_chunks )
+                    #pparams.extend( [param]*num_chunks )
+                    for i in range(num_chunks):
+                        presult_list.append( self.result_list[id][(i*chunk_size+offset):((i+1)*chunk_size+offset)] )
+            except Exception as e:
+                #sys.stderr.write('run_mappers(): caught exception while trying start MapReduce Task: {0}\n'.format(e))
+                #sys.stderr.write('run_mappers(): self.parameters={0}\n'.format(self.parameters))
+                #sys.stderr.write('run_mappers(): self.result_list={0}\n'.format(self.result_list))
+                raise
 
-            seed_list = []
-            for _ in range(len(self.parameters)):
-                #need to do it this way cause the number of run per chunk might not be even
-                seed_list.extend(range(self.seed_base, self.seed_base+number_of_trajectories, chunk_size))
-                self.seed_base += number_of_trajectories
-            #def run_ensemble_map_and_aggregate(model_class, parameters, seed_base, number_of_trajectories, mapper, aggregator=None):
-            results  = self.lv.map_async(run_ensemble_map_and_aggregate, [self.model_class]*num_pchunks, pparams, param_set_ids, seed_list, pchunks, [mapper]*num_pchunks, [aggregator]*num_pchunks)
+            self.running_MapReduceTask = self.lv.map_async(map_and_aggregate, presult_list, param_set_ids, [mapper]*num_pchunks,[aggregator]*num_pchunks,[cache_results]*num_pchunks)
+            self.save_state()
+        else:
+            #sys.stderr.write('run_mappers(): MapReduce Task already running\n')
+            pass
+        
 
 
         if progress_bar:
-            # This should be factored out somehow.
             divid = str(uuid.uuid4())
             pb = HTML("""
                           <div style="border: 1px solid black; width:500px">
@@ -697,38 +887,44 @@ def run(self, mapper, aggregator=None, reducer=None, number_of_trajectories=None
                           </div>
                           """.format(divid))
             display(pb)
+            
+            #sys.stderr.write('run_mappers(): running_MapReduceTask.ready()={0}\n'.format(self.running_MapReduceTask.ready()))
+            while not self.running_MapReduceTask.ready():
+                self.running_MapReduceTask.wait(timeout=1)
+                progress = 100.0 * self.running_MapReduceTask.progress / len(self.running_MapReduceTask)
+                display(Javascript("$('div#%s').width('%f%%')" % (divid, 100.0*(+1)/len(self.running_MapReduceTask))))
+            display(Javascript("$('div#%s').width('%f%%')" % (divid, 100.0)))
+            #sys.stderr.write('run_mappers(): running_MapReduceTask.ready()={0}\n'.format(self.running_MapReduceTask.ready()))
+        else:
+            #sys.stderr.write('run_mappers(): waiting for MapReduce Task to complete (no progress bar)\n')
+            self.running_MapReduceTask.wait()
 
-        # We process the results as they arrive.
-        mapped_results = {}
-        for i,rset in enumerate(results):
+    
+        #sys.stderr.write("\n\nself.running_MapReduceTask.result={0}\n\n".format(self.running_MapReduceTask.result))
+        # Process the results.
+        cnt=0
+        self.mapped_results = {}
+        for i,rset in enumerate(self.running_MapReduceTask.result):
+            
             param_set_id = rset['param_set_id']
             r = rset['result']
-            if param_set_id not in mapped_results:
-                mapped_results[param_set_id] = []
+            if param_set_id not in self.mapped_results:
+                self.mapped_results[param_set_id] = []
             if type(r) is type([]):
-                mapped_results[param_set_id].extend(r) #if a list is returned, extend that list
+                self.mapped_results[param_set_id].extend(r) #if a list is returned, extend that list
             else:
-                mapped_results[param_set_id].append(r)
-            if progress_bar:
-                display(Javascript("$('div#%s').width('%f%%')" % (divid, 100.0*(i+1)/len(results))))
-
-        if verbose:
-            print "Running reducer on mapped and aggregated results (size={0})".format(len(mapped_results[0]))
-        if reducer is None:
-            reducer = builtin_reducer_default
-        # Run reducer
-        return self.run_reducer(reducer, mapped_results)
-
-
-
-    def run_reducer(self, reducer, mapped_results):
-        """ Inside the run() function, apply the reducer to all of the map'ped-aggregated result values. """
-        return reducer(mapped_results[0], parameters=self.parameters[0])
-
-
-
+                self.mapped_results[param_set_id].append(r)
+            cnt+=len(r)
+            
+        if cnt != num_results_to_compute*len(self.parameters):
+            raise MolnsUtilException('run_mappers() num_results_to_compute={0} len(parameters)={1}, got {2} results'.format(num_results_to_compute, len(self.parameters), cnt))
+        
+        self.number_of_results = self.number_of_trajectories
+        # Set the state to not running
+        self.running_MapReduceTask = None
+        self.save_state()
 
-    #--------------------------
+    #-----------------------------------------------------------------------------------
     def add_realizations(self, number_of_trajectories=None, chunk_size=None, verbose=True, progress_bar=True, storage_mode="Shared"):
         """ Add a number of realizations to the ensemble. """
         if number_of_trajectories is None:
@@ -739,37 +935,38 @@ def add_realizations(self, number_of_trajectories=None, chunk_size=None, verbose
         if chunk_size is None:
             chunk_size = self._determine_chunk_size(number_of_trajectories)
 
-        if not verbose:
-            progress_bar=False
-        else:
+        if verbose:
             if len(self.parameters) > 1:
                 print "Generating {0} realizations of the model at {1} parameter points (chunk size={2})".format(number_of_trajectories, len(self.parameters), chunk_size)
             else:
                 print "Generating {0} realizations of the model (chunk size={1})".format(number_of_trajectories,chunk_size)
+            
+        if self.running_SimulationTask is None:
+            #sys.stderr.write('add_realizations(): Starting Simulation Task\n')
+            num_chunks = int(math.ceil(number_of_trajectories/float(chunk_size)))
+            chunks = [chunk_size]*(num_chunks-1)
+            chunks.append(number_of_trajectories-chunk_size*(num_chunks-1))
+            # total chunks
+            pchunks = chunks*len(self.parameters)
+            num_pchunks = num_chunks*len(self.parameters)
+            pparams = []
+            param_set_ids = []
+            for id, param in enumerate(self.parameters):
+                param_set_ids.extend( [id]*num_chunks )
+                pparams.extend( [param]*num_chunks )
 
-        self.number_of_trajectories += number_of_trajectories
-
-        num_chunks = int(math.ceil(number_of_trajectories/float(chunk_size)))
-        chunks = [chunk_size]*(num_chunks-1)
-        chunks.append(number_of_trajectories-chunk_size*(num_chunks-1))
-        # total chunks
-        pchunks = chunks*len(self.parameters)
-        num_pchunks = num_chunks*len(self.parameters)
-        pparams = []
-        param_set_ids = []
-        for id, param in enumerate(self.parameters):
-            param_set_ids.extend( [id]*num_chunks )
-            pparams.extend( [param]*num_chunks )
-
-        seed_list = []
-        for _ in range(len(self.parameters)):
-            #need to do it this way cause the number of run per chunk might not be even
-            seed_list.extend(range(self.seed_base, self.seed_base+number_of_trajectories, chunk_size))
-            self.seed_base += number_of_trajectories
-        results  = self.lv.map_async(run_ensemble, [self.model_class]*num_pchunks, pparams, param_set_ids, seed_list, pchunks, [storage_mode]*num_pchunks)
+            seed_list = []
+            for _ in range(len(self.parameters)):
+                #need to do it this way cause the number of run per chunk might not be even
+                seed_list.extend(range(self.seed_base, self.seed_base+number_of_trajectories, chunk_size))
+                self.seed_base += number_of_trajectories
+            self.running_SimulationTask  = self.lv.map_async(run_ensemble, [self.model_class]*num_pchunks, pparams, param_set_ids, seed_list, pchunks, [storage_mode]*num_pchunks)
+            self.save_state()
+        else:
+            #sys.stderr.write('add_realizations(): Simulation Task already running\n')
+            pass
 
         if progress_bar:
-            # This should be factored out somehow.
             divid = str(uuid.uuid4())
             pb = HTML("""
                           <div style="border: 1px solid black; width:500px">
@@ -777,23 +974,49 @@ def add_realizations(self, number_of_trajectories=None, chunk_size=None, verbose
                           </div>
                           """.format(divid))
             display(pb)
+            
+            #sys.stderr.write("add_realizations(): running_SimulationTask.ready()={0}\n".format(self.running_SimulationTask.ready()))
+            while not self.running_SimulationTask.ready():
+                self.running_SimulationTask.wait(timeout=1)
+                progress = 100.0 * self.running_SimulationTask.progress / len(self.running_SimulationTask)
+                display(Javascript("$('div#%s').width('%f%%')" % (divid, 100.0*(+1)/len(self.running_SimulationTask))))
+            display(Javascript("$('div#%s').width('%f%%')" % (divid, 100.0)))
+            #sys.stderr.write("add_realizations(): running_SimulationTask.ready()={0}\n".format(self.running_SimulationTask.ready()))
+        else:
+            #sys.stderr.write('add_realizations(): waiting for Simulation Task to complete (no progress bar)\n')
+            self.running_SimulationTask.wait()
 
+
+        #sys.stderr.write("\n\nself.running_SimulationTask.result={0}\n\n".format(self.running_SimulationTask.result))
         # We process the results as they arrive.
-        for i,ret in enumerate(results):
+        cnt=0
+        for i,ret in enumerate(self.running_SimulationTask.result):
             r = ret['filenames']
             param_set_id = ret['param_set_id']
             if param_set_id not in self.result_list:
                 self.result_list[param_set_id] = []
             self.result_list[param_set_id].extend(r)
-            if progress_bar:
-                display(Javascript("$('div#%s').width('%f%%')" % (divid, 100.0*(i+1)/len(results))))
+            cnt+=len(r)
+        #sys.stderr.write('add_realizations(): stored {0} simulation results\n'.format(cnt))
+        if cnt != number_of_trajectories*len(self.parameters):
+            raise MolnsUtilException('add_realizations() number_of_trajectories={0} len(parameters)={1}, got {2} results'.format(number_of_trajectories, len(self.parameters), cnt))
+            
+        wall_time = self.running_SimulationTask.wall_time
+        serial_time = self.running_SimulationTask.serial_time
+        self.number_of_trajectories += number_of_trajectories
+        self.running_SimulationTask = None
+        self.save_state()
 
+        #sys.stderr.write('add_realizations(): self.result_list has {0} keys\n'.format(len(self.result_list)))
 
-        return {'wall_time':results.wall_time,'serial_time':results.serial_time}
 
+        return {'wall_time':wall_time,'serial_time':serial_time}
 
 
+
+    #-----------------------------------------------------------------------------------
     #-------- Convenience functions with builtin mappers/reducers  ------------------
+    #-----------------------------------------------------------------------------------
 
     def mean_variance(self, mapper=None, number_of_trajectories=None, chunk_size=None, verbose=True, store_realizations=True, storage_mode="Shared", cache_results=False):
         """ Compute the mean and variance (second order central moment) of the function g(X) based on number_of_trajectories realizations
@@ -816,13 +1039,17 @@ def histogram_density(self, g=None, number_of_trajectories=None):
             in the ensemble. """
         raise Exception('TODO')
 
-    #--------------------------
+    #-----------------------------------------------------------------------------------
+    #-----------------------------------------------------------------------------------
 
     def _update_client(self, client=None):
+        """ Setup the IPython.parallel.Client() object. """
         if client is None:
-            self.c = IPython.parallel.Client()
+            self.c = IPython.parallel.Client(profile='default')
         else:
             self.c = client
+        if len(self.c.ids) == 0:
+            raise MolnsUtilException("No parallel engines detected.  Molnsutil will not work in this situation.  Do you need to start more MOLNs workers?")
         self.c[:].use_dill()
         if self.num_engines == None:
             self.lv = self.c.load_balanced_view()
@@ -836,19 +1063,37 @@ def _update_client(self, client=None):
                 engines = self.c.ids[:self.num_engines]
                 self.lv = self.c.load_balanced_view(engines)
 
+        # Start a spin thread
+        self.c.stop_spin_thread()
+        self.c.spin_thread(interval=10)
         # Set the number of times a failed task is retried. This makes it possible to recover
         # from engine failure.
         self.lv.retries=3
 
     def _determine_chunk_size(self, number_of_trajectories):
         """ Determine a optimal chunk size. """
-        return int(max(1, round(number_of_trajectories/float(self.num_engines))))
+        return min(int(max(1, round(number_of_trajectories/float(self.num_engines)))),1)
 
     def _clear_cache(self):
         """ Remove all cached result objects on the engines. """
         pass
         # TODO
 
+    def clear_results(self):
+        """ Remove the output from the Mappers and Reducers from the storge."""
+        self.delete_results()
+        self.mapper_fn = None
+        self.aggregator_fn = None
+        self.reducer_rn = None
+        self.step3_complete = False
+        self.save_state()
+
+    def delete_results(self):
+        """ Delete final results of the computation. """
+        self.step2_complete = False
+        self.mapped_results = {}
+        self.number_of_results = 0
+
     def delete_realizations(self):
         """ Delete realizations from the storage. """
         if self.storage_mode is None:
@@ -864,60 +1109,83 @@ def delete_realizations(self):
                     ss.delete(filename)
                 except OSError as e:
                     pass
+        self.step1_complete = False
+        self.result_list = {}
+        self.number_of_trajectories = 0
 
     def __del__(self):
         """ Deconstructor. """
-        try:
-            self.delete_realizations()
-        except Exception as e:
-            pass
-
+        pass
 
+############################################################################
 class ParameterSweep(DistributedEnsemble):
     """ Making parameter sweeps on distributed compute systems easier. """
+    
+    my_class_name = 'ParameterSweep'
 
-    def __init__(self, model_class, parameters, client=None, num_engines=None):
+    def __init__(self, name=None, model_class=None, parameters=None, client=None, num_engines=None, ignore_model_mismatch=False):
         """ Constructor.
         Args:
-          model_class: a class object of the model for simulation, must be a sub-class of URDMEModel
+          model_class: a class object of the model for simulation
           parameters:  either a dict or a list.
             If it is a dict, the keys are the arguments to the class constructions and the
               values are a list of values that argument should take.
               e.g.: {'arg1':[1,2,3],'arg2':[1,2,3]}  will produce 9 parameter points.
             If it is a list, where each element of the list is a dict
             """
+        self.parameters = None
 
-        DistributedEnsemble.__init__(self, model_class, parameters, client, num_engines)
-
-        self.my_class_name = 'ParameterSweep'
-        self.parameters = []
-        
-        # process the parameters
-        if type(parameters) is type({}):
-            vals = []
-            keys = []
-            for key, value in parameters.items():
-                 keys.append(key)
-                 vals.append(value)
-            pspace=itertools.product(*vals)
-
-            paramsets = []
-
-            for p in pspace:
-                pset = {}
-                for i,val in enumerate(p):
-                    pset[keys[i]] = val
-                paramsets.append(pset)
-
-            self.parameters = paramsets
-        elif type(parameters) is type([]):
-            self.parameters = parameters
-        else:
-            raise MolnsUtilException("parameters must be a dict.")
-
+        if not isinstance(name, str):
+            raise MolnsUtilException("name not specified")
+        self.name = name
+        if not inspect.isclass(model_class):
+            raise MolnsUtilException("model_class not a class")
+        self.model_class = cloudpickle.dumps(model_class)
         # Set the Ipython.parallel client
         self.num_engines = num_engines
-        self._update_client()
+        self._update_client(client)
+        if not self.load_state(ignore_model_mismatch=ignore_model_mismatch):
+            # State not found, set defaults
+            self.number_of_trajectories = 0
+            self.seed_base = self.generate_seed_base()
+            self.storage_mode = None
+            self.result_list = {}
+            self.running_MapReduceTask = None
+            self.running_SimulationTask = None
+            self.reduced_results = None
+            self.mapped_results = None
+            self.number_of_results = 0
+            self.step1_complete = False
+            self.step2_complete = False
+            self.step3_complete = False
+            self.mapper_fn = None
+            self.aggregator_fn  = None
+            self.reducer_fn = None
+
+        if self.parameters is None:
+            # process the parameters
+            if type(parameters) is type({}):
+                vals = []
+                keys = []
+                for key, value in parameters.items():
+                     keys.append(key)
+                     vals.append(value)
+                pspace=itertools.product(*vals)
+
+                paramsets = []
+
+                for p in pspace:
+                    pset = {}
+                    for i,val in enumerate(p):
+                        pset[keys[i]] = val
+                    paramsets.append(pset)
+
+                self.parameters = paramsets
+            elif type(parameters) is type([]):
+                self.parameters = parameters
+            else:
+                raise MolnsUtilException("parameters must be a dict.")
+
 
     def _determine_chunk_size(self, number_of_trajectories):
         """ Determine a optimal chunk size. """
@@ -926,18 +1194,21 @@ def _determine_chunk_size(self, number_of_trajectories):
             return number_of_trajectories
         return int(max(1, math.ceil(number_of_trajectories*num_params/float(self.num_engines))))
 
-    def run_reducer(self, reducer, mapped_results):
+    def run_reducer(self, reducer):
         """ Inside the run() function, apply the reducer to all of the map'ped-aggregated result values. """
         ret = ParameterSweepResultList()
         for param_set_id, param in enumerate(self.parameters):
-            ret.append(ParameterSweepResult(reducer(mapped_results[param_set_id], parameters=param), parameters=param))
+            ret.append(ParameterSweepResult(reducer(self.mapped_results[param_set_id], parameters=param), parameters=param))
         return ret
     #--------------------------
 
 
 
 class ParameterSweepResult():
-    """TODO"""
+    """ Result object to encapsulate the results of a parameter sweep at a particular parameter point along with the parameters used to compute it.  
+    
+    The object has two member variables, 'result' which contains the MapReduce output and 'parameters' which is a dict containing the parameters used to compute the simulation trajectories.
+    """
     def __init__(self, result, parameters):
         self.result = result
         self.parameters = parameters
@@ -946,6 +1217,7 @@ def __str__(self):
         return "{0} => {1}".format(self.parameters, self.result)
 
 class ParameterSweepResultList(list):
+    """ List encapsulating the ParameterSweepResult objects. """
     def __str__(self):
         l = []
         for i in self:
diff --git a/setup.py b/setup.py
index d3e55af..caa1a9a 100644
--- a/setup.py
+++ b/setup.py
@@ -1,5 +1,5 @@
 from setuptools import setup
 
 setup(name='molnsutil', 
-      version="0.0.2",
+      version="0.8",
       packages=['molnsutil'])