Using NumPy and SciPy with MsPASS#

Overview#

The NumPy and SciPy packages are heavily used for scientific computing in a wide range of fields. One of the reasons for that is performance. Built-in Python containers are not designed for high-performance numerical linear algebra. NumPy provides typed multidimensional arrays and compiled numerical kernels; SciPy builds additional scientific algorithms on that foundation. ObsPy also uses NumPy arrays for waveform samples. As a result, an early design constraint in MsPASS was to provide as clean an interface to NumPy and SciPy as possible. However, because we chose to implement our own core data objects in C++, that presents a disconnect that causes a few anomalies in the MsPASS API. This section of the User’s Manual documents these issues and provides guidance for mixing NumPy and SciPy algorithms with TimeSeries, Seismogram, and ensemble objects. The experimental regular-grid Gather and SeismogramGather types use array-oriented storage directly and are outside this page’s scope.

TimeSeries Data#

Recall that the MsPASS TimeSeries class represents a single channel of fixed sample rate seismic data spanning a given time period. As such we used the standard signal processing abstraction of storing the data as a vector of numbers.

Every computing language that has any role in numerical computing defines a syntax for arrays. A vector in this context is a one-dimensional array. A Python list is not the same thing as an array in Fortran or C. Both Fortran and C define an array as a fixed block of memory of a specified length containing a sequence of numbers of the same type. (e.g. a C double array of length 10 is a block of \(8 \times 10\) bytes.) NumPy performance is achieved by redefining an array to assume the data are in a contiguous memory block and using wrappers that allow low-level operations to actually be performed by compiled C or FORTRAN functions. We do the same thing in MsPASS but our data objects do not contain NumPy array classes directly. Instead, we use the same approach to allow low-level operations to be done in C++. Because we use the same concept (arrays are contiguous blocks of memory) the arrays can mostly be used interchangeably, but not always.

So what is the difference? If you execute this code fragment in MsPASS:

import numpy as np
from mspasspy.ccore.seismic import DoubleVector, TimeSeries
x_numpy = np.zeros(100)
x_mspass = TimeSeries(100)
print("Type of numpy object=",type(x_numpy))
print("Type of mspass data vector=",type(x_mspass.data))

You should get this output:

Type of numpy object= <class 'numpy.ndarray'>
Type of mspass data vector= <class 'mspasspy.ccore.seismic.DoubleVector'>

The point is that both define two conceptually similar things that are a vector of sample data. In this example these are stored with the symbols x_numpy and x_mspass.data. The print statements emphasize they are different data types. They are not, however, an apples-to-oranges kind of comparison but more of a Jonathan versus Gala apple comparison. Python APIs often use “duck typing”: an operation accepts objects that provide the required behavior rather than demanding one concrete class. Because of “duck typing” we can, in most cases, interchange the use of the two classes ndarray and DoubleVector. For example, you can usually send DoubleVector data to a NumPy or SciPy function that accepts array-like input. Here are a few examples:

import numpy as np
from scipy.signal import hilbert
# Setting data to nonzero values - not essential but otherwise
# trivial zero vector will cause problems for the log example
for i in range(x_mspass.npts):
  x_mspass.data[i]=float(i+1)
x = np.log10(x_mspass.data)
y = np.sqrt(x_mspass.data)
z=np.cos(x_mspass.data)
yh = hilbert(x_mspass.data)

That works because DoubleVector is what NumPy’s documentation calls array-like. In fact, any NumPy or SciPy function that accepts an array-like argument can generally accept the data member of a TimeSeries. In reality, DoubleVector is a C++ “vector container” that in C++ would be declared as:

std::vector<double> x;

In MsPASS we use stock “bindings” from pybind11 that allow us to do this behavior. The buffer protocol also lets np.asarray(x_mspass.data) create a writable NumPy view without copying. Mutating that view therefore mutates the TimeSeries samples; call np.array(x_mspass.data, copy=True) when isolation is required.

This works for most, but not all, vector operators that mix DoubleVector and ndarray types. For example, with the symbols above defined all of the following will work and do what is expected:

z=x_numpy-x_mspass.data
z=x_mspass.data+x_numpy
x_mspass.data -= z
z=x_numpy*x_mspass.data  # a peculiar numpy operation - not a dot product
x_mspass.data *= 10.0

The main thing to avoid is an assignment where the left hand side and right hand side resolve to different types. For instance, the following will fail with a TypeError exception:

z=np.ones(100)
x_mspass.data = x_numpy
x_mspass.data = x_numpy + z

The solution is to construct a DoubleVector explicitly:

x_mspass.data = DoubleVector(x_numpy)
x_mspass.data = DoubleVector(x_numpy + z)

That is necessary because the right-hand side is an ndarray. The DoubleVector constructor copies and converts the values.

Be aware that if you mix types in a vector operation you can get some surprising results. For example, the following code generates a TypeError exception:

ts=TimeSeries(100)
z=np.ones(100)
ts.data = z + ts.data

while reversing the order of z and ts.data

ts=TimeSeries(100)
z=np.ones(100)
ts.data = ts.data + z

works and sets ts.data to the vector sum of the two symbols on the right-hand side. The result type follows the implementation selected by Python’s binary-operator dispatch, so operand order matters. Avoid relying on that detail; make the conversion explicit:

ts=TimeSeries(100)
z = np.ones(100)
x = z + np.asarray(ts.data)
ts.data = DoubleVector(x)

Current DoubleVector bindings support ordinary Python slices. A slice returns a new DoubleVector:

ts = TimeSeries(100)
x = ts.data[5:10]
assert isinstance(x, DoubleVector)
assert len(x) == 5

NumPy is still more convenient for advanced indexing. The important invariant is that ts.npts must equal len(ts.data). Assigning a shorter vector does not update npts automatically, so update both:

ts = TimeSeries(100)
x = np.ones(100)[5:10]
ts.set_npts(len(x))
ts.data = DoubleVector(x)
assert ts.npts == len(ts.data) == 5

Seismogram Data#

MsPASS uses a different data type to hold data from three-component sensors. The motivation is described in the section titled Data Object Design Concepts. The result is that the array holding the sample data is a matrix instead of a vector. Some implementation details of note are:

  1. In MsPASS the matrix has 3 rows and npts columns. That means the rows of the matrix can be extracted to create a TimeSeries subset of the data while the columns are single “vector” samples of ground motion.

  2. Currently MsPASS uses a simple, lightweight matrix class to store the sample data called dmatrix. A key reason we made that choice was control over methods defined for the class. A case in point is the shape attribute, which is important when working with NumPy.

  3. For performance dmatrix is implemented in C++. The class has methods for standard matrix operations with Python bindings for those operators. When working with Seismogram data, use the native operators when possible.

  4. A dmatrix stores the array data in a contiguous block of memory in what NumPy calls Fortran order. Fortran stores a matrix with the row index “varying fastest”. Exchanging the sample data stored this way with NumPy is like that with the std::vector used in TimeSeries. np.asarray(seismogram.data) provides a writable view of that buffer.

  5. We use a variation of the same “magic” in the pybind11 binding code that allows the Seismogram.data matrix to interact cleanly with NumPy and SciPy functions that require a matrix. Like the scalar case there are impedance mismatches, however, that can complicate that exchange. We discuss these below.

Almost all of what was discussed above about using TimeSeries arrays with NumPy are similar. For example, although the result is meaningless, you can run the following code snippet with MsPASS and it will run and produce the expected result:

import numpy as np
from mspasspy.ccore.seismic import Seismogram
from mspasspy.ccore.utility import dmatrix
x_mspass = Seismogram(100)
# initialize the matrix to some meaningless, nonzero values
for i in range(3):
  for j in range(x_mspass.npts):
     x_mspass.data[i,j]=i+j+1
# these apply the math operation for each number in the data matrix
z=np.cos(x_mspass.data)
z=np.exp(x_mspass.data)
z=np.sqrt(x_mspass.data)
# matrix operation that multiplies all samples by 2.0
x_mspass.data *= 2.0
# matrix operator adding content of numpy matrix z to data matrix
y=x_mspass.data + z
# matrix operation adding y and x_mspass.data matrices
x_mspass.data += y
# Same thing done with a binary + operator
# use a copy to show more likely use
z=Seismogram(x_mspass)
x_mspass.data = z.data + y

As with TimeSeries, a mismatch occurs if the operation yields a type mismatch. For example, the following minor variant of above will run:

z=Seismogram(x_mspass)
y=np.ones([3,100])
x_mspass.data = z.data + y

In contrast, the following that might seem completely equivalent raises a TypeError exception:

z=Seismogram(x_mspass)
y=np.ones([3,100])
x_mspass.data = y + z.data

The reason is that for the right hand side resolves to a NumPy ndarray and the left-hand side requires a dmatrix. Exactly like above, this simpler assignment will fail exactly the same way:

x_mspass.data = y

The solution is similar to that we used above for TimeSeries:

x_mspass.data = dmatrix(y)

That is, DoubleVector is replaced by dmatrix for matrix data.

dmatrix also supports ordinary row and column slices through its NumPy interface; those slices are returned as ndarray objects rather than as another dmatrix. When a result is assigned back, construct a dmatrix explicitly, require exactly three rows, and keep npts equal to the number of columns:

window = np.array(x_mspass.data, copy=True)[:, 10:20]
x_mspass.set_npts(window.shape[1])
x_mspass.data = dmatrix(window)
assert x_mspass.data.shape == (3, x_mspass.npts)