Accumulators that != reduce

[pavlis]

I am becoming convinced we have a general class of problems for which we need a generic solution in MsPASS. That problem is one that acts like a reduce operator, in some sense, but for which I am pretty sure neither dask nor spark can handle. I have run across two examples in recent weeks:

1. This one is actually a generic algorithm we need to consider adding to MsPASS. Given: bag of TimeSeries or Seismogram objects produced by a MongoDB query followed by read_distributed data. Want: an output TimeSeriesEnsemble or SeismogramEnsemble. The example where I hit this was a workflow where I wanted to read common source gathers, window the raw data around a phase, demean the data, and then convert the result to a TimeSeriesEnsemble that I can use as input to bundle. I'm not certain how to do this other than something like item 2 I'll relate next.
2. In working on pwmig this came up as part of the migration algorithm. That particular algorithm produces a bag of Seismogram like objects that are depth-migrated data with some computed weights used for amplitude normalization and weights used as part of the "Generalized Radon Transform" that is a geeky concept inside of pwmig. The step inside pwmig needs to take a bag of all the migrated data and stuff them into slots in a large 4D grid used to assemble one "plane wave" that is eventually summed in a more standard reduce operation (subject to a memory and efficiency issue noted elsewhere). The point is this also has the model of a bag (RDD in spark) of one type of object that need to be reduced to a different type. The operations satisfy the reduce requirement of being communitive and associatve, but I can'f see how to resolve the type conflicts in the api? i.e. the binary reduce operator in spark and dask both seem to dogmatically demand that both args used as input and the output must all be the same type.

First, is my understanding of this correct? Maybe there is some other function the dask or spark api that can handle this, but I've not found it.

The solution I used in pwmig is not ideal, but is a workaround. It is also a dask specific solution because I don't think spark has the equivalent of delayed. Rather than show the pwmig code for mspass it might be more informative to show a rough implementation of part of problem 1 above (take a bag of TimeSeries objects and produce a TimeSeriesEnsemble.)

# code above here defines mydata as a dask bag of TimeSeries objects
ensemble = TimeSeriesEnsemble()
delayed_data=mydata.to_delayed()   # this produces a list of length defined by partitioning
for dgroup in delayed_data:
   dlist=dgroup.compute()
   for d in dlist:
      ensemble.member.append(d)   # this is the actual accumulator - some parallelism with partitions but likely bottlenecked here

We could use that approach to build a generic bag to ensemble function for dask. I don't know how to do the equivalent in spark. That particular functionality probably should be included in MsPASS.

Is there a more generic solution to this problem?

[pavlis]

Correction to above. Could have edited that box, I suppose, but preserving this to emphasize this point is probably important.

Correction is that the generic issue above breaks spark and dask because the binary operator is not associative. In normal arithmetic, addition is associative, for example, because c=b+a=a+b. When a and b are not simple types that cannot be converted, as in the examples above, that rule is broken.

We can cast this as a different abstract rule. We need a generalization of c+=a where c and a can be different types but for which the c+=a operation makes sense (could literally mean it is defined). For example 1 above c is an ensemble object and a is one of our atomic data types (TimeSeries or Seismogram). For example 2, c is the grid object into which data are accumulated and "a" is the seismogram like data migrated from time to depth.

[pavlis]

This discussion has been inactive for months and I want to bring it to the front of our discussion. In working on updates to the User Manual I realized the functionality discussed in this page defines some essential tools we need to have in MsPASS.

1. We need to solve the problem I raised initially here. How can you take a bag/RDD of TimeSeries or Seismogram objects and "reduce" them to a (always smaller) bag/RDD of TimeSeriesEnsemble or SeismogramEnsemble objects respectively? I have not figured out a way to do that. I think the answer may be in "GroupBy" constructions but I'm not sure.
2. A simpler, but related problem is how to most efficiently create a bag of ensembles as a read operation? I have a solution I'm putting in the user's manual, but we may want to make a generic function to implement the approach. Here is what I do in pseudocode and know this works:

1) run a function to create a list of queries (e.g. source_id's to create common source gathers)
2) convert the query list to a bag/RDD
3) map using function to run queries and return one ensemble per query (empty, dead ensembles if the query fails)
4)  rest of the workflow

We need to come to closure on both of these issues. Propose we start in our next weekly meeting and preserve our decision here for the record.

[pavlis]

For the record in our call yesterday I agreed to add to this discussion to design the api for problem 2 in the box above: standard methods for parallel loading of a data set of ensembles driven by a list of something small enough to always fit in memory. The type example would be common source gathers defined at startup by a list of source_ids.

I think first this problem can be made more generic in the same way we have implemented editing and normalization: using a base class and a standard api for a class that defines the query generator. Here is an initial suggestion of the basic content of the base class:

def query_generator:
   """
   abstract base class.
   """
  def __init(self, db,collection="wf_TimeSeries"): 
    """
    With this design the base class would probably at least need a definition of what collection 
    the class should access.   Probably best to pass the top level database handle, db, and 
    have the constructor cache the collection handle.  Done below for this design discussion.
    """
    self.dbcol = db[collection]
  def __call__(self, x):
    """
    would be useful here for simplified syntax as with editor and normalize.   
    """
    self.read_ensemble(x)
  @abstractmethod
  def generate(self,x):
    """
    Takes content of input x, which would normally be an element in a map call or a loop for 
    the comparable serial job, and returns a query string.  
    """
 @abstractmethod
  def read_ensemble(self,x):
    """
    This method takes the value of x, calls the generate method, and then goes one step further 
    and applies the query and read and returns the ensemble created by Database.read_ensemble 
    with the resulting cursor.  Note we should not have the possible intermediate method that 
    return a cursor.  I think such an algorithm could create cursors that would timeout. 
    """
  @abstractmethod
  def ensemble_size(self,x):
    """
    This method returns the ensmble size that would result from x.   Before starting any long running   
     workflow driven by this approach users should be advised to run this method on the dataset 
     to check for huge ensembles and unexpected null ensembles.  This method would run the 
    generate method and then run count_documents instead of calling the ensemble reader.  
   """

This design allows a generic algorithm to load based on this pseudocode using a dask construct:

1) generator = example_generator_instance(some_constructor_input)
2)  generate list of data items expected as input to query_geneator instance (e.g. list of source_ids)=dlist
3) mybag= dask.bag.from_sequence(dlist)
4) mybag = mybad.map(generator.read_ensemble)

The concrete implementations of this base class I think we should produce immediately are:

1. Common source gather generator. It would expect x to be a source_id. Have done this one in
   form that isn't this generic already and will be easy to implement.
2. Common receiver gather by id. Very similar to common source algorithm but would likely need some
   way to switch between channel_id and site_id.
3. Common receiver gather by seed code. This would use net, sta, chan, and (when defined) loc to
   generate the query. May want this to only function with wf_miniseed but it would absolutely
   need to default to wf_miniseed.
4. Time interval query. This would be used for finding all or a subset of data within a specified
   time interval. This is a key algorithm for people doing things like noise correlation or phased
   array processing of continuous data. We currently have nothing like this but it is a fundamental
   part of the framework that is needed.

Looking forward to feedback on this design. Also think about what other implementations might be helpful. In writing this I realzed item 4 above (the time interval query) is a huge hole for a large class of algorithms.