GATHER = subset of Ensemble design

[pavlis]

Purpose of this discussion page

We have a number of issues floating around about what we have been calling “ensembles” in MsPASS, but which are often also called a “gather” in seismic reflection processing. The purpose of this discussion page is to help the group (and community who might read this) understand some key ideas and conceptual models of data related to how an ensemble/gather should be defined. i.e. this is a higher level design document that we should use to flesh out the abstractions we should use of those “conceptual models” of an ensemble/gather.

Definitions

Let me do a couple internal definitions that might help make this discussion more tractable. Humans name things to distinguish different concepts so let me make this pair of definitions for this discussion:

• I will define a gather as the data model used in seismic reflection processing wherein the data can be represented by a matrix.
• I will define an ensemble as the more general definition we have used till now in MsPASS

Each section below expands what I mean in these definitions.

gather concepts

The way a gather is treated in seismic reflection processing comes from a combination of the history of instrumentation and the data model that emerged in the field for efficiency. All seismic reflection systems have relics of the “multichannel” data recorders that were universal for all production systems until about 30 years ago. Multichannel systems used a single digitizers that cycled through “channels” in a single sample cycle. The natural data output was equivalent to writing a matrix in row order (multichannel data). A big issue in the earliest days of seismic reflection processing (i.e. 1960 and early 1970s) that is almost laughable in today’s world of GBy computers was “demultiplexing”. Demultiplexing is mathematically identical to a matrix transpose of a Fortran storage array. Demux algorithms were more or less that but were a challenge when a single shot gather could not fit in a computer’s memory. The point is that a field “shot gather” was universally conceptualized as a matrix from the 1960s until today. An individual “trace” was one column in the matrix. For a “shot gather” each columns was the output of one sensor for one shot. It was thus natural to organize all other “gathers” into a matrix. e.g. a “CMP gather” is a set of column vectors from a set of shot gather matrices created by just sorting the vectors into “CMP bins”. Furthermore, because the only language of the day was FORTRAN the data order was naturally frozen as FORTRAN order (row index varies fastest).

There are three fundamental assumptions in gather processing that make that approach possible:

1. All data have exactly the same sample rate.
2. All data have the same time standard. In reflection seismology that is a relative time standard where component 0 of each vector is time 0 and that time is always assumed to be “shot time” that in earthquake seismology we would call the “origin time”. A geeky detail, as an important aside I will return to below, is that multichannel data always had a one sample ambiguity in the timing and with a granularity of dt/Number_channels ambiguity in the timing created by cyclic sampling.
3. All data have the same time span = number of rows in the data acquisition matrix. Note the segy file format has this assumption for an entire segy data set.

ensemble concepts

For the record, the roots of MsPASS are my seispp library that students and I developed from around 2005 till around 2018. That library provided a prototype for the mspasspy.ccore.seismic modules. We implemented an ensemble in a more generic way that does not require the the assumptions above. i.e. THE fundamental difference between a MsPASS ensemble and a seismic reflection "gather" is that it does not require ANY of the three assumptions above that are implicit in seismic reflection gather processing.

The way we accomplished that is to define two atomic data objects: (1) TimeSeries objects and (2) Seismogram objects. The former is used for scalar data while the later are bundled three-component data. Each atomic object can have a different sampling interval or number of samples than any other ensemble member. Further, there is no assumption about the time standard. Data can be defined by unix epoch times or as some relative time standard and can be mixed up (although that would usually be a bad idea). Readers should note that an ensemble is then stored in a C++ std::vector container defined by the symbol member in the class definition. The same name is bound to python using pybind11 in a way that the member symbol can be manipulated with subscripts and iterators.

To be clear, this is the class template prototype for an ensemble object:

template <typename Tdata> class Ensemble : public mspass::utility::Metadata
{
public:
  /*! \brief Container holding data objects.

  An ensemble is defined as a collection of data objects linked by
  some property that is normally defined by data sharing a comomon
  data-value pair.  The data are stored here in a vector container.
  We make this vector public because the assumption is the member
  objects can be large and copying with a getter would be a serious
  inefficiently.
  */
  std::vector<Tdata> member;

  ... rest of api definition follows
}

member is the std::vector container holding atomic data objects. The template argument, Tdata, is TimeSeries for a TimeSeriesEnsemble and Seismogram for a SeismogramEnsemble. Note also our Ensemble has a Metadata container used to hold metadata common to all atomic data "members".

It is also noteworthy that pybind11 has allowed us to bind the data arrays in both TimeSeries and Seismogram objects to be accessible from numpy functions. TimeSeries sample data is made to look like a vector (one-dimensional ndarray) while the sample data of a Seismogram is bound to a 3xnsample numpy array in Fortran order (each column of the matrix is a 3-component vector data sample at a time that can be computed by the vector index).

Finally, for this design discussion it is important to recognize that a patch we made after a year+ development in MsPASS is to redefine an ensemble to implement two useful concepts implemented in atomic data objects: (1) a boolean that allows the entire ensemble to be "killed" and (2) an ErrorLogger attribute (symbol elog) that allows ensemble related error messages to be posted without aborting a parallel job. For ease of reading this is portion of the class template that defines our LoggingEnsemble that pybind11 binds to the two functional names TimeSeriesEnsemble and SeismogramEnsemble:

template <typename T> class LoggingEnsemble : public Ensemble<T>
{
public:
  /*! Standard mspass container for holding error logs. */
  mspass::utility::ErrorLogger elog;

  ... class methods are defined here ...

private:
  bool ensemble_is_live;
};

ensembles versus gather issues

The concepts defining a gather versus a MsPASS ensemble can be summarized in one phrase: a gather is a special case of an ensemble. That is, a standard reflection seismology gather could be represented by aTimeSeriesEnsemble if we impose three restrictions matching the assumptions
noted above:

1. All members must have a common sample rate.
2. The first sample of each member must well-defined time standard. The simplest is the reflection model where "time 0" (time of the first sample) is the shot/origin time. We have previously used two other standards: (a) a fixed UTC time (a constant time window) and (b) what I seispp called an "arrival time reference" meaning the start time of each member is a constant time offset from some measured or theoretical phase arrival time. There may be others we haven't thought of but the common concept for a gather is it must be a reproducible time standard.
3. All ensemble members must have the same number of samples. With restriction 1 (common dt) that also means the time duration of each member is fixed.

Only if a TimeSeriesEnsemble satisfies the above restrictions can its sample data array be abstracted as a matrix.

A driving question is this: why is there a need for a gather abstraction given the now solid code base for TimeSeriesEnsemble data objects? There are at least three reasons I can think of:

1. Efficiency for some algorithms. Because of the importance of linear algebra inside the data model of so many computing problems there are many algorithms that are highly optimized to work on matrices and vectors. A prime example is numpy that is more or less a set of python bindings on well-tuned linear algebra operators.
2. Dask has a parallel array type that can be used to parallelize large memory operations efficiently an automatically. Many imaging operations in seismology would likely benefit by an abstraction that would allow the user to use dask arrays when appropriate.
3. There are a class of open-source algorithms using the gather model that would be easier to adapt to the framework with an appropriate API. A type example is the large seismic unix package. It has a number of C functions that might be relatively easy to adapt to MsPASS with the right API. There are probably others in the noise-correlation community.
4. Export and import of data from any package that uses the gather model could be standardized if we had a Gather object implemented.

A final issue to address is how a ThreeComponentEnsemble could be mapped into a gather concept? I would suggest that the answer is actually fairly simple: a ThreeComponentEnsemble that satisfies the gather
assumption could readily be represented as a three-dimensional array. The only question would be which index varies fastest? Addressing that question is best left to later, but it is important for efficiency and how hard it would be to adapt any legacy code. Also numpy may simply the whole thing as they have stock methods to reorganize array data. I would also note that implementing a gather model for three-component data should probably be viewed as a low priority as I know of no standard package that it needs to be adapted to the framework that would profit from the distinction between an ensemble and gather as defined here.

Functional Specifications for Gather Implementation

First, I propose we agree on my definition of a gather versus and ensemble. If so, we should aim to use a namespace with the top level working class called Gather. If, as I propose below, we choose to utilize inheritance as a key design element we can expect some intermediate class names like BasicGather or CoreGather to be useful. There might also be subclasses of Gather that are specialized. e.g. as I'll suggest below we might want something like an SUGather that has frozen features mated with Seismic Unix
code.

I would propose the following functional requirements for the new Gather API:

• A Gather should match as much as possible the API for ensemble objects. i.e. all ensemble methods that match in concept with a gather should be implemented in the Gather object. That can assure maximum interoperability. It is to be determined where the mismatch in concept occurs.
• A Gather should be less dogmatic about sample data type than ensemble. A bit of a flaw in our design for ensemble objects is the dogmatic use of 64 bit floats for sample data. To simplify adapting legacy code we should allow support for other numeric data types.
• It must be easy for the user to convert a TimeSeriesEnsmble into a Gather and vice-versa. That is probably best done as a functions rather than as class constructors to make it easy to run the function in a map operation.
• We need a way to easily switch a gather from an numpy array to a dask array. Some design tests will be needed to determine the best default partitioning for dask arrays.
• The class definition will probably need to have a way to repartition a dask array for different algorithms. Imaging algorithms, for example, would likely benefit from block partitions while most operations in the current MsPASS that are signal-by-signal oriented would be best blocked just that way (i.e. a block is defined as on column or n blocks per column for truly huge data block like a continuous data set).
• (optional) We may want to allow a way to efficiently define gaps. Front and tail gaps of a few samples could be a common feature in a Gather.

Implementation Proposal

First a couple of axioms:

• We need to first prototype most if not all of the ideas here in python. We should implement parts in C/C++ only if necessary for performance.
• This is an object design problem. Hence, a key question is the inheritance structure we should use to maximize reuse.

I'm not exactly sure if the following would work, but let me suggest the following as a base class for Gather objects:

from mspasspy.ccore.seismic import BasicTimeSeries
from mspasspy.ccore.utility import Metadata,ProcessingHistory

class BasicGather(BasicTimeSeries,Metadata,ProcessingHistory):

I really have no idea if pybind11 objects can be used in an inheritance tree like this, but it at least defines the base functionality I think a base class for a Gather should contain. Metadata and ProcessingHistory create a functionality that is automatically compatible with our exixting LoggingEnsemble aka TimeSeriesEnsemble and SeismogramEnsemble. Using BasicTimeSeries in the above way is a concise, although perhaps confusing, way to add the functionality in BasicTimeSeries of concepts of "kill" and "time standard". The methods of BasicTimeSeries are all appropriate for a gather object - it is just the name that is potentially confusing. The t0, npts, and all the related methods are appropriate for a gather because of the assumptions that define a gather as a subset of an ensemble. The only concept collision I see is that there are cases where the time standard could be ambiguous. If time is UTC then t0 is unambiguous, but if time is relative two gathers may not be equivalent with no way to tell the difference. e.g. a gather in relative time could have t0 defined by event origin time (shot time in the reflection world) or something like the "arrival time reference". Those two would need to be treated very differently. I think we could work around that by either defining some metadata key-value pair with the value defining some time standard or add to this enum class defined in the C++ code base:

enum class TimeReferenceType {
        UTC, /*!< Use an absolute (usually UTC) time base - previously absolute in SEISPP. */
        Relative /*! Time is relative to some other standard like shot time */
};

to make it something like this

enum class TimeReferenceType {
        UTC, /*!< Use an absolute (usually UTC) time base - previously absolute in SEISPP. */
        Relative, /*! Time is relative to some other standard like shot time */
        MeasuredArrivals,  /*! Time is relative to a measured arrival time. */
        ModelArrival  /*! Time is relative to a model-based arrival time.  Use Metadata to define model:method*/
};

I haven't looked to see if that would break anything in the C++ code base but it might. It shouldn't, however, be hard to handle if it does - some if-else constructs might need to be switch-case construct. The bigger unknown is if one can build a python class from multiple inheritance of 3 classes defined by pybind11? I would think it would work, but that is an unknown. If we agree on this idea one of the first things to do is see if that can be done. The key point is the three classes I suggested above get us a long way to building a functional object with a long list of useful methods we wouldn't need to implement.

I next suggest the working Gather object would be defined as a simple subclass of BasicGather but adding a data array and associated stuff. Here is a rough prototype:

class Gather(BasicGather):
   def __init__(self,use_dask=False, dtype=numpy.double, otherargs):
      """
      Simple constructor to build a skeleton gather to which data is added
      later. dtype arg defines the numpy data type and would allow data
      stored as other numberic types.
      """
      if  use_dask:
         self.data = call to dask array constructor
      else:
          self.data = call to numpy array constructor
  def __init__(self,ensemble):
     """
     Overloaded constructor to create from a TimeSeriesEnsemble - details to flesh out.  As noted 
      above this might not even be desirable and this might be better set up as a function.   TBD
     """
  def change_dtype(dtype=numpy.single):
     """
     In place change of data type to type defined by dtype arg.
     """
  def convert_to_ensemble()->TimeSeriesEnsemble:

Other constructors are nearly guaranteed to be needed. Others please add comments to suggest other methods that would be necessary core methods.

As suggested above a goal of our design should be to maximize reusability and allow specialization of the more generic Gather object. Here is a rough version of the seismic unix gather idea noted earlier:

class SUGather(Gather):
   def __init__(self,filename):
     """
     Construct from an SU file
     """
     ...
  def __init__(self,Gather):
     """
     Construct from a generic Gather.
     """
  def __init__(self,ensemble):
     """
     Constuct from an TimeSeriesSensemble - again maybe better as a function?
     """

A couple special features of a SUgather that would be locked in are:

• the metadata namespace (Schema) would match that in Seismic Unix. That would require some complexity for namespace mapping in the constructors but I think we have tools for doing that.
• Last I knew SU's core data type was C floats (32 bit floats) so any constructor would need to guarantee that coversion.
• Clean interaction with SU C functions might require some way to get a C pointer to the self.data array.
• How this could be used to adapt SU algorithms is a totally open question. There may be some dark alleys here we don't want to step into.

[pavlis]

Quick addition to this discussion. I just poked around the source code for Seismic Unix and don't think there is any merit to creating an SUgather object. If someone wanted to adapt an SU algorithm they would likely be best advised to produce pybind11 wrapper for the SU C code and build a python wrapper around the C code to interact with a gather. SU has way too many interwoven libraries and a long list of main programs using the unix filter model that would make adapting it to a python world clumsy to impossible in my opinio.

What might be useful, however, is thinking about a generic interface to run SU unix filters on a Gather. Such a filter would only need a metadata schema translator to set required SU "header" values. Then the model I'd use to run an SU program would be to make an SU file image, create a fifo, launch an SU string of programs reading from the fifo, push the image of the file into the fifo, and pull the output out and put it into a Gather container.

@wangyinz what other packages were you thinking of that drove you to suggest we needed to implement what I've here called the gather concept?

[pavlis]

For anyone reading this outside the mspass-team you should be aware this conversation was held for a while around a google doc an discussions during weekly meeting. I (@pavlis ) sat down and turned current ideas into a sketch of a prototype python API. I think this forum will make it easier to refine suggestions than the google doc.

First, this design makes use of inheritance. A problem I faced in trying to define how to do that is a huge impedance mismatch between how inheritance works in python and C++. Why that matters is that I was trying to make this interface as much like the C++/pybind11 API for atomic data in MsPASS (TimeSeries and Seismogram) and commonality with MsPASS ensemble objects that are actually implemented as C++ templates. Things that would make sense in C++ would be crazy to impossible in python. For that reason I am here proposing the following inheritance structure for this new data object:

1. I choose to create a base class with the name BasicGather. A BasicGather is abstract and cannot be instantiated but actually contains most methods in my proposal. More fundamental is that the array implementation should be hidden in this base class. As I suggest the job of creating the work space be delegated to this base constructor. Then concrete class constructors should fill in the data array.
2. I propose a Gather be a subclass of BasicGather. Its constructor should take input of either another Gather or a TimeSeriesEnsemble to be created.
3. Similarly is propose a SeismogramGather be a subclass of BasicGather. Its constructor should take input as either another SeismogramGather or a SeismogramEnsemble.

A few clarifications that are side issues from the above that this design includes:

• I've included what are essentially copy constructor functionality because I am not sure if large matrix implementtions like dask work work with something like "deepcopy". I doubt it, but I haven't researched the topic.
• You will see I propose another feature in this design for efficiency in creating one of these things in pieces. That is, I think the append method should behave like C++ stdlib containers. That is things like std::vector containers will accept an open ended number of calls to the push_back method and the container will be reallocated and grow as needed. One at time operators like that on a large array could be truly horrible so I suggest a blocksize argument wherein appends will not force a resize until that size is passed. Not sure if that is possible with the various implementations for large matrix arrays we have been discussing but this design has that. I'm putting it here as a topic for discussion in this forum as we finalize this design.

That said, here is my sketched design for the base class:

from abc import ABC, abstractmethod

class BasicGather(ABC):
    """
    Base class for new MsPASS ensemble class that utilizes array 
    implementations to store the data.  The implementtion may or may not 
    provide for arrays that are too large to fit in memory.  
    That behavior should be controlled by parameters on constructin.  
    
    The data object defined by this class and it's children needs to be
    understood as appropriate only for ensemble data that satisfy several 
    restrictions that the current MsPASS ensemble objects 
    (TimeSeriesEnsemble and SeismogramEnsemble) do not rqequire:
        1.  All the data have the same sample rate.
        2.  Each signal has the same number of samples.   
        3.  The data cannot have gaps.   Data with gaps neeed to be discarded 
            or patched (e.g. with a gap fill) before being loaded into 
            one of these contaienrs. 
    
    Methods for this base class are largely getters and setters for 
    univeral parameters that match the above.  It must be understood that 
    with this desig the base class is an incomplete skeleton that 
    is fleshed out by subclasses.   This clsss has some virtual methods 
    (defined with the ABC @abstractmethod decorator) and a set of methods 
    that would never work if only the base class were contructed
    (not allowed because of the virtual meethods).  In particular, the 
    way the base class handle metadata is potnetially confusing because
    the base class constructor only defines stubs for the data used to 
    hold member and ensemble metadata  
    
    A new design feature suggested here not in our earlier documents is the 
    idea of dynamic, periodic resizes for the append method.  The append
    method and the description of resize_blocksize in the constructor.  
     
    
    A design issue for discussionis is whether not a gather should implement 
    the ator and rota functionality of BasicTimeSeries.   The cost is 
    tiny up front and retrofitting after that fact might cause a lot of 
    headaches.  REcommend we implement them as defined below.
    
    """
    
    def __init__(self,
                 number_members,
                 npts,
                 number_components=1,
                 dt=None,
                 array_type="numpy",
                 resize_blocksize=10
               ):
        """
        This base class constructor should take the view that the base 
        class is the fraework with no data.  That is in line with 
        BasicTimeSeries and an appropriate launching for subclasses that
        define scalar and 3C data as conrete implementations.  
        A difference here, however, is that constuctor will create the 
        work array but just not populate it.  
    .
              
        :param number_members: Expected number of signals in the ensemble
          to be loaded.  That is, the number of columns in the array 
          container.   The idea is that we standardize creation of the 
          work array in this base class.   Concrete implementations are
          expected to estimte the size of array needed and then call 
          this constuctor before filling the array with data based on 
          the input to the subclass constructor.  e.g. for MsPASS 
          ensembles that would be set by len(ensemble.member).  
          The array implementtion is set by the related parameter 
          array_size.  This parameter is required and there is no
          default due to the expected use only by subclasses
        :type number_members:  integer
        
        :param npts:  length in samples of all members.  This parameter 
          is required and there is no default due to the expected use only 
          by subclasses
        :type npts:  integer
        
        :param dt:  sample interval of all data in ensemble.  Default is 
          None which taken to mean it should be derived from the data.  
        :type dt:  float
        
        :param array_type:  sets the way large arrays are implemented.  
          Default should by "numpy" which means the sample array will 
          be defined as a nptsXnumber_members FORTRAN order numpy array.
          Other obvious choices are "dask" or "zarr"  List of acceptable 
          options to be determined
         :type array_type:  string that must match keywords TBD
         
        :param resize_blocksize:  Resizing a container like ths is always 
          going to be an expensive operation.  The container needs to behave 
          like C++ stdlib containers where there is a hidden blocksize 
          that triggers periodic reallocs when the data buffer fills.   
          This parameter is a proposed way to implement the same concept.
          It would mainly be used for the append method.   Whenever append 
          hits the array wall it would resize the array adding resize_blocksize
          columns (signal slots), copy the current data, append the new, 
          and destroy the old.  Obviously a potentially very expenseive 
          operation but a very useful feature. This name sucks but the 
          concept is what is important
        :type resize_blocksize:  integer
          
        """
        # these need additional code to handle autosetting but for now 
        # define them like this to ake below clearer
        
        self.npts = npts
        self.dt = dt
        self.array_type == array_type
        # we need an error log
        self.elog = ErrorLogger()   
        # this a global boolean used to mark the entire ensemble invalid 
        # when false.
        self.ensemble_is_marked_live=False
        self.live=[]  # boolean array for members
        self.t0=[] # starttime array (maybe should be an ndarray)
        # This illustrates the idea of how number_members should be used
        # Not complete but illustrtes the idea
        if input_data is None:
            if number_members is None:
                raise RuntimeError("usage message")
            self.number_members=0
            self.array_member_size=number_members
            self.data = _private_create_function(self.npts,self.array_member_size)
        elif isinstance(input_data,TimeSeriesEnsmeble):
            self.number_members = len(input_data.member)
            self.array_member_size = self.number_members
            self.data = _private_create_function(elf.npts,self.array_member_size)
        elif isinstance(input_data,SeismogramEnsemble):
            elf.number_members = len(input_data.member)
            self.array_member_size = self.number_members
            self.data = _private_create_function_for_3C(elf.npts,self.array_member_size)
        else:
            raise TypeError("Input is not a vaid type")
            
        
    def append(self,mspass_object,otherargs):
        """
        Appends data in mspass_object to internal array object.  
        Needs to handle sizing similar to std containers in C++ 
        where allocs (in this case a resize) is triggered by filling.
        Constructors need to also handle this concept - see above.  
        The resize increment is resize_blocksize.  When append filsl the 
        current space (self.number_members == self.array_member_size)
        a resize will be triggered.  
        """
        pass
    def starttime(self)->list:
        """
        Return a list (pythona array) of starttimes.  We should assume 
        times are always stored as doubles.
        """
        pass
    def column_values(self)->list:
        """
        Returns the values assigned to the column axis.  These should 
        default to integers matching the column index values.   They 
        should be allowed to be anything, however, to handle things 
        like variable offset record sections. 
        """
        pass
    def set_column_values(self,x):
        """
        Set the column values array to the values in x. 
        """
        pass

    
    # The following names match BasicTimeSeries because the match in 
    # concept.  starttime above doesn't really because in a single object 
    # t0==starttime is constant.  
    def time(self,i,j):
        """
        Return time of array position i,j (Note works the same for 3c and 
        scalar data)
        """
        pass
    def sample_number(self,time,x)->tuple:
        """
        Return a tuple of ints of the index position of an ordered pair defined by a time 
        (row) value and a column value that can be mapped to the 
        column index.   
        """
        pass
    def time_is_UTC(self):
        """
        Returns True if the time standard for the data is set as 
        UTC
        """
        pass
    def time_is_relative(self):
        """
        Returns True if the time standard for the data is set as 
        relative
        """
        pass
    def live(self,j)->bool:
        """
        Return true if jth member is set live.
        """
    def dead(self,j)->bool:
        """
        Return true if jth member is marked dead (opposite of live).
        """
    # This next set are s in BasicTimeSeries but all
    def samprate(self):
        return 1.0/self.dt
    def rtoa(self):
        pass
    def ator(self,shift):
        pass
    def shift(self,timeshift):
        pass
    # The following methods in BasicTimeSeries are optional.  It isn't 
    # clear to me if their use is required - at least initially.  With 
    # pythona ddingg methods is not as probematic as a Java or C++
    # They are;  time_reference, timetype, set_to, set_tref, and 
    # set_npts.  
    
    # An alternative design would be to put these in an intermediate 
    # class but don't think that would be useful for a python implementation
    def metadata(self,j,return_type="dict"):
        """
        Return the Metadata components of member j.  return_type should 
        allow optional return as dict or Metadata.  Not clear which should 
        be default.
        """
        pass
    def set_metadata(self,j,md):
        """
        Setter for metadata of member j.   md is the new continer that would 
        replace current conent.  Most useful for constructors.
        """
        pass
    def edit_metadata(self,j,md):
        """
        Differs from set_metadata in that the contents of md are added to 
        the current and do not fully repalce them.   Needed for updating 
        metadata after costruction
        """
        pass
        
    @abstractmethod
    def member(self,j,data_only=False):
        """
        Returns the data associated with member j.   Virtual method 
        returns different type for scalar ersus 3c data.  
        data_only in this context is only a hint to a concrete 
        implementation to have an option to return just the data vector/matrix
        """
        pass
    @abstractmethod
    def subset(self):
        """
        Return a subset of the data defined by either a collection of 
        members or a restricted sample range (a form of windowing).
        Concrete implementtions can define what this means.   Minimum 
        usage would be something functionally equivalent to the f90/matlab
        synatx of ensemble.member[i:j].   We may also want to be less 
        generic and require addition methods with other names like 
        "extract_members" or "window" (time range).  
        """
        pass
    def ensemble_metadata(self,return_type="dict"):
        """
        Like metadata method but returns the ensemble metadata.  
        """
        pass
    def sync_metadata(self):
        """
        This would act like the one in Ensemble in C++ copying all the 
        ensemble key-value pairs to the members.  It is debatable that 
        this would be needed but would be trivial to implement.
        """
        pass
    # not sure if we can make this member abstract.  It needs to 
    # be in this design to allow different signatures for scalar and 3c data
    @abstractmethod
    def __getitem__(self):
        """
        Python version of operator[].   Not sure the declaration as 
        abstractmethod will work or is even deirable.  I put it hat 
        way for now as we do want to require thiat method.  For scalar
        data it should get two indices while for 3C data it needs 3.
        """
        pass
    @abstractmethod
    def __setitem__(self):
        """
        Setter by index.   Same issue that is probably isn't necesessy
        in he implementtion to use the abstractmethod eecorator. 
        """
        pass

I emphasize this is mostly just method signatures and docstrings, but it is more concrete than what we've been discussing on github.

I next propose this for single component data that I suggest we call a Gather:

class Gather(BasciEnsembleArray):
    """
    Concrete implementtion for a scalar gather.  This is the array 
    equivalent of a TimeSeriesEnsemble appropriate when the input matches
    the concept of BasicGather.    It follows the OOP stndard 
    pardigm that creation is initialization.  This puts a lot of features 
    in the constructor for the class.
    """
    def __init__(self,mspass_object,
                 number_members=None,
                 npts=None,
                 dt=None,
                 array_type="numpy",
                 resize_blocksize=10):
        """
        Because in this design the base class only sets up the workspace
        most of the nitty gritty work of building this thig will go here.  
        
        :param mspss_object:  Should contain the data to be loaded.  
        :type mspass_object:  should accept TimeSeriesEnsemble or another 
          Gather object.  Would advise a different class if we wanted to 
          allow soemthig like a raw matrix as input.  We perhaps should 
          accept a None and in that situation just create the workspace
          of a speciied size and assume the gather will be build by 
          many calls to append.
        
        :param number_members:  expected number of signals the gather 
          will eventually contain.  If None (default) the size is 
          determined by the input ensemble.  If nonzero it should be 
          a positive integer larger than the size of the ensemble passed 
          as mspass_object.  
        :type:  integer
        
        :param npts:  expected number of samples for all data in the gather 
          If None, which is the default, it should be estimated from the 
          input data.  Note mspass_object beng set None and this parameter 
          or number_members set None should cause an exception to be 
          raised.  Either that or it should be ignored and a warning 
          message posted - probably a better idea.
          
        dt, array_type, and resize_blocksize are as defined in the 
        base clsss.  They are passed directly to the base class 
        construtor 
        """
        # this is a rough prototype showing how the base clsss constructor 
        # should be used.
        if mspass_object is None:
            if number_members is None or npts is None:
                raise RuntimeError("illegal input ")
            super().__init__(number_members,
              npts,
              number_components=1,
              dt=dt,
              array_type=array_type,
              resize_blocksize=resize_blocksize
            )
        elif isinstance(mspass_object,TimeSeriesEnsemble):
            enssize=len(mspass_object.member)
            if number_members is None:
                nmem_to_use=enssize
            else:
                if number_members<enssize:
                    nmem_to_use=enssize
                else:
                    nmem_to_use=enssize
            for d in mspass_object.member:
                if d.live:
                    npts_to_use = d.npts
                    break
            super().__init__(nmem_to_use,
              npts_to_use,
              number_components=1,
              dt=dt,
              array_type=array_type,
              resize_blocksize=resize_blocksize
            )
            # elo is assumed created in base clsss consructor
            self.elog.log_error("npts and input disageement",ErrorSeverity.Warning)
            
            # Additional code to initialie here to load metadata and any
            # other stuff
        elif isinstance(mspass_object, Gather):
            # this maybe could do with deepcopy but suspect the disk 
            # arrays might not work
            super().__init__(
                  mspass_object.number_members,
                  mspass_object.npts,
                  number_components=1,
                  dt=mspass_object.dt,
                  array_type=mspass_object.array_type,
                  resize_blocksize=mspass_object.resize_blocksize
                )
            #additional code here to load other attributes to copy
        else:
            raise RuntimeError("illegal input error")
            
    def member(self,j)->TimeSeries:
        """
        Returns the data associated with member j.  Unlike above I suggest
        we not do the idea of "data_only".  Aways return a TimeSeries.
        Suggest data method for that purpose
        """
        pass  # stub
    def data(self,j)->ndarray:
        """
        Return the raw data vector associated with column j.   Defined here 
        as an ndarray return but probably should be any iterable
        container that acts like a vector.  Not clear what will be needed to 
        support large matrix formats.
        """
        pass

    def subset(self,start,end)->Gather:
        """
        Return a subset of the Gather with signals from start to end 
        (like start:end in F90 or matlab).  
        
        TBD is if subset should be overloaded to allow other selections
        """
        pass
    
    def __getitem__(self,i,j):
        """
        Index operator to return data value at sample index i,j.
        """
        pass
    def __setitem__(self,i,j,newvalue):
        """
        Index setter - I think this is the right syntax for two indices.
        """
        pass

And this painfully similar for this rough outline for Seismogram input:

class SeismogramGather(BasicGather):
    """
    Gather for three-component data.  Can only be created from inputs that 
    satisfy the restrictions of BasicGather.  It is much like Gather 
    but with a 3D instead of a 2D array holding sample data. 
    """
    def __init__(self,mspass_object,
                 number_members=None,
                 npts=None,
                 dt=None,
                 array_type="numpy",
                 resize_blocksize=10):
        """
        Because in this design the base class only sets up the workspace
        most of the nitty gritty work of building this thig will go here.  
        
        :param mspss_object:  Should contain the data to be loaded.  
        :type mspass_object:  should accept TimeSeriesEnsemble or another 
          Gather object.  Would advise a different class if we wanted to 
          allow soemthig like a raw matrix as input.  We perhaps should 
          accept a None and in that situation just create the workspace
          of a speciied size and assume the gather will be build by 
          many calls to append.
        
        :param number_members:  expected number of signals the gather 
          will eventually contain.  If None (default) the size is 
          determined by the input ensemble.  If nonzero it should be 
          a positive integer larger than the size of the ensemble passed 
          as mspass_object.  
        :type:  integer
        
        :param npts:  expected number of samples for all data in the gather 
          If None, which is the default, it should be estimated from the 
          input data.  Note mspass_object beng set None and this parameter 
          or number_members set None should cause an exception to be 
          raised.  Either that or it should be ignored and a warning 
          message posted - probably a better idea.
          
        dt, array_type, and resize_blocksize are as defined in the 
        base clsss.  They are passed directly to the base class 
        construtor 
        """
        # this is a rough prototype showing how the base clsss constructor 
        # should be used.
        if mspass_object is None:
            if number_members is None or npts is None:
                raise RuntimeError("illegal input ")
            super().__init__(number_members,
              npts,
              number_components=3,
              dt=dt,
              array_type=array_type,
              resize_blocksize=resize_blocksize
            )
        elif isinstance(mspass_object,SeismogramEnsemble):
            enssize=len(mspass_object.member)
            if number_members is None:
                nmem_to_use=enssize
            else:
                if number_members<enssize:
                    nmem_to_use=enssize
                else:
                    nmem_to_use=enssize
            for d in mspass_object.member:
                if d.live:
                    npts_to_use = d.npts
                    break
            super().__init__(nmem_to_use,
              npts_to_use,
              number_components=3,
              dt=dt,
              array_type=array_type,
              resize_blocksize=resize_blocksize
            )
            # elo is assumed created in base clsss consructor
            self.elog.log_error("npts and input disageement",ErrorSeverity.Warning)
            
            # Additional code to initialie here to load metadata and any
            # other stuff
        elif isinstance(mspass_object, SeismogramGather):
            # this maybe could do with deepcopy but suspect the disk 
            # arrays might not work
            super().__init__(
                  mspass_object.number_members,
                  mspass_object.npts,
                  number_components=3,
                  dt=mspass_object.dt,
                  array_type=mspass_object.array_type,
                  resize_blocksize=mspass_object.resize_blocksize
                )
            #additional code here to load other attributes to copy
        else:
            raise RuntimeError("illegal input error")
            
    def member(self,j)->Seismogram:
        """
        Returns the data associated with member j.  Unlike above I suggest
        we not do the idea of "data_only".  Aways return a TimeSeries.
        Suggest data method for that purpose
        """
        pass  # stub
    def data(self,j)->ndarray:
        """
        Return the raw data matrix associated with column j.   Defined here 
        as an ndarray return but probably should be any iterable
        container that acts like a matrix.  Not clear what will be needed to 
        support large matrix formats.
        """
        pass

    def subset(self,start,end)->SeismogramGather:
        """
        Return a subset of the Gather with signals from start to end 
        (like start:end in F90 or matlab).  
        
        TBD is if subset should be overloaded to allow other selections
        """
        pass
    
    def __getitem__(self,i,j,k):
        """
        Index operator to fetch sample from time axis at i, member axis at j, 
        and component number k.   
        """
        pass
    def __setitem__(self,i,j,k,newvalue):
        """
        Indexing setter.   Not sure this is the right syntax but should show 
        the idea.
        """

I think the strongest thing about this design is it cleanly separates scalar and 3C data in a manner as consistently as I can come up with. I've tried to model the methods to be as consistent with BasicTimeSeries and our existing Ensemble classes as possible.

I'm going to close this comment. Please put comments on this basic design below this box as replies I'm going to add two more boxes that are independent topics you should comment on separately.

[Yangzhengtang]

Some add-on on the implementation:

1. About resize_blocksize: there are two concepts invovled: 1) in a distributed/parallel setting, the array data are distributed into smaller chunks, we need a mechanism to control the chunksize 2) we need a mechanism to control the growth of the capacity to avoid expensive copy. For 1), I think we can have a field "num_partition" so that chunksize=total_size/num_partition. For 2), the explicit control over the growth might be redundant, since the underlying array implementations (dask, zarr, numpy...) should have some optimization for those situations. Further research is needed for the performance in different implementations, and for now we can have a field "growth_rate" to let the capacity grow expontentially just like the std::vector.
2. The copy constructor that accepts a Gather object in the Gather class might be unnecessary. In Python, the shallow/deep copying can be controlled by the "copy" module (here).

[pavlis]

In composing the above and reviewing various notes I had an idea this encompasses. It build on an earlier comment about seismic unix. Let me just put down the class with a rough docstring that I think introduces the idea reasonably well:

class SUgather(Gather):
     """
     This is an example of a subclass that is based of a reflection 
     format that matches the restrctions of a BasicGather. 
       
     
     This needs some work, but this could be a slick way to provide a
     mspass interface with seismic unix.   It could contain methods to 
     exec a line of an SU program.   Worth considering
     """
     def __init__(self,
                  file_like_handle,
                  array_type="dask"):
         """
         This type of gather could be construted directly from an SU 
         data file defined with file_like_handle.  
         """
         
     def exec(arg):
         """
         An idea here is we could use some variat of exec to execute a 
         single SU comamdn or a pipeline of commands through this method. 
         I am not sure if thatis possible but I put it out there as an 
         idea.
         """

This is very rough but it boils down to two ideas:

1. An SUgather is a natural child of the more Gatherobject. It's metadta namespace would be frozen as SU segy header fields. More importantly it would be constucted from a seismic unix file. Writing that constructor would not be trivial as SU now uses an XDR float format. It would probably require some C++ functions and pybind11 to implement. The key point, howeer, is it shows the merit of this inheritance structure. An SUgather is a Gather. It just has a special, fixed namespace and is created from a file.
2. I stub in an idea of aexec method that could be used to pass the data array in the gather through an SU processing pipeline. It is one idea for how to have the framework work with legacy package like SU

[pavlis]

There are two low-level algorithms we have defined for current MsPASS data objects for which we need the comparable functionality for these new data objects: (1) time windowing and (2) ExtractComponent for 3C data. I first thought both might be more appropriate as methods but should remain in algorithms. Neither are trivial. If you accept my proposal to allow the new family of gather objects to handle relative and UTC time like BasicTimeSeries, time windowing has complexities we need to discuss. I'm pretty certain ExtractComponent would have a dependence on the large matrix package implementation. Point is these two are core algorithms we would need to implement as part of this development. A discussion point is what other things currently in algorithms should accept one of these objects directly and not require a translation to something like a loop over atomic objects?

[pavlis]

Have one last point for this. We perhaps need mechanisms for things like the member method to return obspy Trace or Stream objects instead of MsPASS data objects. Could streamline obspy algorithms at the cost of losing the history and error logging capability of the MsPASS data objects. I can see both sides of the argument. Suggest we not consider that until someone out there specifically asks for it. I don't think it would be a difficult retrofit. Constructing a TimeSeries and a Trace from the contents of a Gather are not really very different.