Reject misaligned grids in waveform arithmetic

cxx/include/mspass/seismic/CoreSeismogram.h

@@ -223,8 +223,8 @@ class CoreSeismogram : public mspass::seismic::BasicTimeSeries,
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to add to this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreSeismogram &operator+=(const CoreSeismogram &d);
   /*! Addition operator.
@@ -250,8 +250,8 @@ class CoreSeismogram : public mspass::seismic::BasicTimeSeries,
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to subract from this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreSeismogram &operator-=(const CoreSeismogram &d);
   /*! Subtraction operator.

cxx/include/mspass/seismic/CoreTimeSeries.h

@@ -144,8 +144,8 @@ This method is complicated by the need to sync the changed value with
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to add to this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreTimeSeries &operator+=(const CoreTimeSeries &d);
   /*! Addition operator.
@@ -172,8 +172,8 @@ This method is complicated by the need to sync the changed value with
   the ends of the rhs are inside the range of the lhs.
 
   \param d is other signal to subract from this.
-  \exception MsPASSError can be thrown if lhs and rhs do not have matching
-  time standards.
+  \exception MsPASSError if lhs and rhs have incompatible time standards,
+  sample intervals, or start-time grids.
   **/
   CoreTimeSeries &operator-=(const CoreTimeSeries &d);
   /*! Subtraction operator.

cxx/src/lib/seismic/CoreSeismogram.cc

@@ -3,6 +3,7 @@
 #include "mspass/seismic/keywords.h"
 #include "mspass/utility/MsPASSError.h"
 #include "mspass/utility/SphericalCoordinate.h"
+#include "waveform_arithmetic.h"
 #include <boost/any.hpp>
 #include <float.h>
 #include <math.h>
@@ -913,82 +914,36 @@ CoreSeismogram &CoreSeismogram::operator*=(const double scale) {
   return (*this);
 }
 CoreSeismogram &CoreSeismogram::operator+=(const CoreSeismogram &d) {
-  int i, iend, jend;
-  size_t j, i0, j0;
   // Silently do nothing if d or lhs is marked dead
   if (d.dead() || (this->dead()))
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreSeismogram += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j) {
-    this->u(0, j) += d.u(0, i);
-    this->u(1, j) += d.u(1, i);
-    this->u(2, j) += d.u(2, i);
+  const auto overlap =
+      detail::arithmetic_overlap(*this, d, "CoreSeismogram::operator+=");
+  for (std::size_t i = 0; i < overlap.count; ++i) {
+    const std::size_t lhs_index = overlap.lhs_begin + i;
+    const std::size_t rhs_index = overlap.rhs_begin + i;
+    this->u(0, lhs_index) += d.u(0, rhs_index);
+    this->u(1, lhs_index) += d.u(1, rhs_index);
+    this->u(2, lhs_index) += d.u(2, rhs_index);
   }
   return (*this);
 }
 /* IMPORTANT:  this code is absolutely identical to that for operator+=
 except the += in the last loop becomes -=.  Any changes in operator+=
 must have exactly the same change here (other than a message with a
 tag to the function)*/
-CoreSeismogram &CoreSeismogram::operator-=(const CoreSeismogram &data) {
-  int i, iend, jend;
-  size_t j, i0, j0;
-  // Sun's compiler complains about const objects without this.
-  CoreSeismogram &d = const_cast<CoreSeismogram &>(data);
+CoreSeismogram &CoreSeismogram::operator-=(const CoreSeismogram &d) {
   // Silently do nothing if d is marked dead
-  if (!d.mlive)
+  if (d.dead())
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreSeismogram += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j) {
-    this->u(0, j) -= d.u(0, i);
-    this->u(1, j) -= d.u(1, i);
-    this->u(2, j) -= d.u(2, i);
+  const auto overlap =
+      detail::arithmetic_overlap(*this, d, "CoreSeismogram::operator-=");
+  for (std::size_t i = 0; i < overlap.count; ++i) {
+    const std::size_t lhs_index = overlap.lhs_begin + i;
+    const std::size_t rhs_index = overlap.rhs_begin + i;
+    this->u(0, lhs_index) -= d.u(0, rhs_index);
+    this->u(1, lhs_index) -= d.u(1, rhs_index);
+    this->u(2, lhs_index) -= d.u(2, rhs_index);
   }
   return (*this);
 }

cxx/src/lib/seismic/CoreTimeSeries.cc

@@ -3,6 +3,7 @@
 #include "mspass/seismic/keywords.h"
 #include "mspass/utility/Metadata.h"
 #include "mspass/utility/MsPASSError.h"
+#include "waveform_arithmetic.h"
 #include <vector>
 namespace mspass::seismic {
 using namespace std;
@@ -65,85 +66,27 @@ CoreTimeSeries &CoreTimeSeries::operator=(const CoreTimeSeries &tsi) {
 /*  Sum operator for CoreTimeSeries object */
 
 CoreTimeSeries &CoreTimeSeries::operator+=(const CoreTimeSeries &data) {
-  int i, iend, jend;
-  size_t i0;
-  size_t j, j0;
-  // Sun's compiler complains about const objects without this.
-  CoreTimeSeries &d = const_cast<CoreTimeSeries &>(data);
-  // Silently do nothing if d is marked dead
-  if (!d.mlive)
+  // Silently do nothing if data is marked dead
+  if (data.dead())
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreTimeSeries += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  // cout << "i0="<<i0<<" j0="<<j0<<" iend="<<iend<<" jend="<<jend<<endl;
-  /*  Now do the actual sum using the computed ranges */
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j)
-    this->s[j] += d.s[i];
+  const auto overlap =
+      detail::arithmetic_overlap(*this, data, "CoreTimeSeries::operator+=");
+  for (std::size_t i = 0; i < overlap.count; ++i)
+    this->s[overlap.lhs_begin + i] += data.s[overlap.rhs_begin + i];
   return (*this);
 }
 /* IMPORTANT:  this code is absolutely identical to that for operator+=
 except the += in the last loop becomes -=.  Any changes in operator+=
 must have exactly the same change here (other than a message with a
 tag to the function)*/
 CoreTimeSeries &CoreTimeSeries::operator-=(const CoreTimeSeries &data) {
-  int i, iend, jend;
-  size_t i0;
-  size_t j, j0;
-  // Sun's compiler complains about const objects without this.
-  CoreTimeSeries &d = const_cast<CoreTimeSeries &>(data);
-  // Silently do nothing if d is marked dead
-  if (!d.mlive)
+  // Silently do nothing if data is marked dead
+  if (data.dead())
     return (*this);
-  // Silently do nothing if d does not overlap with data to contain sum
-  if ((d.endtime() < mt0) || (d.mt0 > (this->endtime())))
-    return (*this);
-  if (d.tref != (this->tref))
-    throw MsPASSError("CoreTimeSeries += operator cannot handle data with "
-                      "inconsistent time base\n",
-                      ErrorSeverity::Invalid);
-  /* this defines the range of left and right hand sides to be summed */
-  i = d.sample_number(this->mt0);
-  if (i < 0) {
-    j0 = this->sample_number(d.t0());
-    i0 = 0;
-  } else {
-    j0 = 0;
-    i0 = i;
-  }
-  iend = d.sample_number(this->endtime());
-  jend = this->sample_number(d.endtime());
-  if (iend >= (d.npts())) {
-    iend = d.npts() - 1;
-  }
-  if (jend >= this->npts()) {
-    jend = this->npts() - 1;
-  }
-  // cout << "i0="<<i0<<" j0="<<j0<<" iend="<<iend<<" jend="<<jend<<endl;
-  /*  Now do the actual sum using the computed ranges */
-  for (i = i0, j = j0; i <= iend && j <= jend; ++i, ++j)
-    this->s[j] -= d.s[i];
+  const auto overlap =
+      detail::arithmetic_overlap(*this, data, "CoreTimeSeries::operator-=");
+  for (std::size_t i = 0; i < overlap.count; ++i)
+    this->s[overlap.lhs_begin + i] -= data.s[overlap.rhs_begin + i];
   return (*this);
 }
 const CoreTimeSeries

cxx/src/lib/seismic/waveform_arithmetic.h

@@ -0,0 +1,67 @@
+#ifndef MSPASS_SEISMIC_WAVEFORM_ARITHMETIC_H
+#define MSPASS_SEISMIC_WAVEFORM_ARITHMETIC_H
+
+#include "mspass/seismic/BasicTimeSeries.h"
+#include "mspass/utility/MsPASSError.h"
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <sstream>
+#include <string>
+
+namespace mspass::seismic::detail {
+struct ArithmeticOverlap {
+  std::size_t lhs_begin;
+  std::size_t rhs_begin;
+  std::size_t count;
+};
+
+inline ArithmeticOverlap arithmetic_overlap(const BasicTimeSeries &lhs,
+                                            const BasicTimeSeries &rhs,
+                                            const char *caller) {
+  using mspass::utility::ErrorSeverity;
+  using mspass::utility::MsPASSError;
+
+  if (lhs.timetype() != rhs.timetype())
+    throw MsPASSError(std::string(caller) +
+                          ": operands use inconsistent time references",
+                      ErrorSeverity::Invalid);
+
+  const double lhs_dt = lhs.dt();
+  const double rhs_dt = rhs.dt();
+  const double dt_tolerance =
+      1.0e-6 * std::max(std::abs(lhs_dt), std::abs(rhs_dt));
+  if (!std::isfinite(lhs_dt) || !std::isfinite(rhs_dt) ||
+      !(std::abs(lhs_dt - rhs_dt) <= dt_tolerance)) {
+    std::ostringstream message;
+    message << caller << ": sample intervals do not match: lhs dt=" << lhs_dt
+            << ", rhs dt=" << rhs_dt;
+    throw MsPASSError(message.str(), ErrorSeverity::Invalid);
+  }
+
+  const double offset_samples = (rhs.t0() - lhs.t0()) / lhs_dt;
+  const double rounded_offset = std::round(offset_samples);
+  if (!(std::abs(offset_samples - rounded_offset) <= 1.0e-6)) {
+    std::ostringstream message;
+    message << caller << ": start times are not aligned to the sample grid: "
+            << "offset=" << offset_samples << " samples";
+    throw MsPASSError(message.str(), ErrorSeverity::Invalid);
+  }
+
+  if (rounded_offset >= static_cast<double>(lhs.npts()) ||
+      rounded_offset <= -static_cast<double>(rhs.npts()))
+    return ArithmeticOverlap{0, 0, 0};
+
+  if (rounded_offset >= 0.0) {
+    const std::size_t lhs_begin = static_cast<std::size_t>(rounded_offset);
+    return ArithmeticOverlap{lhs_begin, 0,
+                             std::min(lhs.npts() - lhs_begin, rhs.npts())};
+  }
+
+  const std::size_t rhs_begin = static_cast<std::size_t>(-rounded_offset);
+  return ArithmeticOverlap{0, rhs_begin,
+                           std::min(lhs.npts(), rhs.npts() - rhs_begin)};
+}
+} // namespace mspass::seismic::detail
+
+#endif

cxx/test/CMakeLists.txt

@@ -1,6 +1,7 @@
 include(CTest)
 if(BUILD_TESTING)
   add_subdirectory(amap)
+  add_subdirectory(arithmetic)
   add_subdirectory(decon_serialization)
   add_subdirectory(data_directory)
   add_subdirectory(dmatrix)
@@ -19,6 +20,8 @@ if(BUILD_TESTING)
 
 # note this test needs to access a data file in the run directory. 
 # we set that with WORKING_DIRECTORY and the special symbol after the keyword
+  add_test(NAME test_waveform_arithmetic
+    COMMAND ${PROJECT_BINARY_DIR}/test/arithmetic/test_waveform_arithmetic)
   add_test(NAME test_decon_serialization COMMAND ${PROJECT_BINARY_DIR}/test/decon_serialization/test_decon_serialization WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/decon_serialization)
   add_test(NAME test_data_directory
     COMMAND ${PROJECT_BINARY_DIR}/test/data_directory/test_data_directory)

cxx/test/arithmetic/CMakeLists.txt

@@ -0,0 +1,8 @@
+add_executable(test_waveform_arithmetic test_waveform_arithmetic.cc)
+include_directories(
+  ${Boost_INCLUDE_DIRS}
+  ${pybind11_INCLUDE_DIR}
+  ${Python_INCLUDE_DIRS}
+  ${PROJECT_SOURCE_DIR}/include/)
+
+target_link_libraries(test_waveform_arithmetic PRIVATE mspass ${Boost_LIBRARIES} ${Python_LIBRARIES})