From 6a5caedc55336477ee3a5d60a049920f4a396801 Mon Sep 17 00:00:00 2001 From: wangyinz Date: Wed, 12 Aug 2026 11:18:51 -0500 Subject: [PATCH 1/4] Validate waveform arithmetic grids --- cxx/include/mspass/seismic/CoreSeismogram.h | 8 +- cxx/include/mspass/seismic/CoreTimeSeries.h | 8 +- cxx/src/lib/seismic/CoreSeismogram.cc | 83 ++----- cxx/src/lib/seismic/CoreTimeSeries.cc | 83 +------ cxx/src/lib/seismic/waveform_arithmetic.h | 67 +++++ cxx/test/CMakeLists.txt | 3 + cxx/test/arithmetic/CMakeLists.txt | 8 + .../arithmetic/test_waveform_arithmetic.cc | 235 ++++++++++++++++++ python/tests/test_waveform_arithmetic.py | 139 +++++++++++ 9 files changed, 492 insertions(+), 142 deletions(-) create mode 100644 cxx/src/lib/seismic/waveform_arithmetic.h create mode 100644 cxx/test/arithmetic/CMakeLists.txt create mode 100644 cxx/test/arithmetic/test_waveform_arithmetic.cc create mode 100644 python/tests/test_waveform_arithmetic.py diff --git a/cxx/include/mspass/seismic/CoreSeismogram.h b/cxx/include/mspass/seismic/CoreSeismogram.h index 247414d67a..685df5359e 100644 --- a/cxx/include/mspass/seismic/CoreSeismogram.h +++ b/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. diff --git a/cxx/include/mspass/seismic/CoreTimeSeries.h b/cxx/include/mspass/seismic/CoreTimeSeries.h index b96b2dd25e..5a3a6d792a 100644 --- a/cxx/include/mspass/seismic/CoreTimeSeries.h +++ b/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. diff --git a/cxx/src/lib/seismic/CoreSeismogram.cc b/cxx/src/lib/seismic/CoreSeismogram.cc index e4cedee251..93a0b4253e 100644 --- a/cxx/src/lib/seismic/CoreSeismogram.cc +++ b/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 #include #include @@ -913,39 +914,17 @@ 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); } @@ -953,42 +932,18 @@ CoreSeismogram &CoreSeismogram::operator+=(const CoreSeismogram &d) { 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(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); } diff --git a/cxx/src/lib/seismic/CoreTimeSeries.cc b/cxx/src/lib/seismic/CoreTimeSeries.cc index b4924e4fa6..d557d6c8dd 100644 --- a/cxx/src/lib/seismic/CoreTimeSeries.cc +++ b/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 namespace mspass::seismic { using namespace std; @@ -65,42 +66,13 @@ 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(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="<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="< +#include +#include +#include +#include + +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(lhs.npts()) || + rounded_offset <= -static_cast(rhs.npts())) + return ArithmeticOverlap{0, 0, 0}; + + if (rounded_offset >= 0.0) { + const std::size_t lhs_begin = static_cast(rounded_offset); + return ArithmeticOverlap{lhs_begin, 0, + std::min(lhs.npts() - lhs_begin, rhs.npts())}; + } + + const std::size_t rhs_begin = static_cast(-rounded_offset); + return ArithmeticOverlap{0, rhs_begin, + std::min(lhs.npts(), rhs.npts() - rhs_begin)}; +} +} // namespace mspass::seismic::detail + +#endif diff --git a/cxx/test/CMakeLists.txt b/cxx/test/CMakeLists.txt index 4a3b795346..e972974af7 100644 --- a/cxx/test/CMakeLists.txt +++ b/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) diff --git a/cxx/test/arithmetic/CMakeLists.txt b/cxx/test/arithmetic/CMakeLists.txt new file mode 100644 index 0000000000..1cd4b5659a --- /dev/null +++ b/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}) diff --git a/cxx/test/arithmetic/test_waveform_arithmetic.cc b/cxx/test/arithmetic/test_waveform_arithmetic.cc new file mode 100644 index 0000000000..3c565dcb5c --- /dev/null +++ b/cxx/test/arithmetic/test_waveform_arithmetic.cc @@ -0,0 +1,235 @@ +#include "mspass/seismic/Seismogram.h" +#include "mspass/seismic/TimeSeries.h" +#include "mspass/utility/MsPASSError.h" +#include +#include +#include +#include +#include +#include +#include + +using mspass::seismic::Seismogram; +using mspass::seismic::TimeReferenceType; +using mspass::seismic::TimeSeries; +using mspass::utility::ErrorSeverity; +using mspass::utility::MsPASSError; + +namespace { +int failure_count{0}; + +void check(const bool condition, const std::string &message) { + if (!condition) { + std::cerr << message << '\n'; + ++failure_count; + } +} + +template struct Samples; + +template <> struct Samples { + static constexpr std::size_t component_count{1}; + static double get(const TimeSeries &data, const std::size_t, + const std::size_t sample) { + return data.s[sample]; + } + static void set(TimeSeries &data, const std::size_t, const std::size_t sample, + const double value) { + data.s[sample] = value; + } +}; + +template <> struct Samples { + static constexpr std::size_t component_count{3}; + static double get(const Seismogram &data, const std::size_t component, + const std::size_t sample) { + return data.u(component, sample); + } + static void set(Seismogram &data, const std::size_t component, + const std::size_t sample, const double value) { + data.u(component, sample) = value; + } +}; + +template +Waveform make_waveform(const std::size_t sample_count, const double t0, + const double dt, const double base) { + Waveform data(sample_count); + data.set_t0(t0); + data.set_dt(dt); + data.set_tref(TimeReferenceType::Relative); + data.set_live(); + data.put("sentinel", std::string("unchanged")); + for (std::size_t component = 0; + component < Samples::component_count; ++component) + for (std::size_t sample = 0; sample < sample_count; ++sample) + Samples::set(data, component, sample, + base + 10.0 * static_cast(component) + + static_cast(sample)); + return data; +} + +template struct Snapshot { + std::size_t sample_count; + double t0; + double dt; + TimeReferenceType time_reference; + bool live; + int error_count; + std::string sentinel; + std::vector samples; +}; + +template Snapshot snapshot(const Waveform &data) { + Snapshot result{data.npts(), + data.t0(), + data.dt(), + data.timetype(), + data.live(), + data.elog.size(), + data.get_string("sentinel"), + {}}; + result.samples.reserve(Samples::component_count * data.npts()); + for (std::size_t component = 0; + component < Samples::component_count; ++component) + for (std::size_t sample = 0; sample < data.npts(); ++sample) + result.samples.push_back(Samples::get(data, component, sample)); + return result; +} + +template +void check_state(const Waveform &actual, const Snapshot &expected, + const std::string &context) { + check(actual.npts() == expected.sample_count, context + ": npts changed"); + check(actual.t0() == expected.t0, context + ": t0 changed"); + check(actual.dt() == expected.dt, context + ": dt changed"); + check(actual.timetype() == expected.time_reference, + context + ": time reference changed"); + check(actual.live() == expected.live, context + ": live state changed"); + check(actual.elog.size() == expected.error_count, + context + ": error log changed"); + check(actual.get_string("sentinel") == expected.sentinel, + context + ": metadata changed"); + std::size_t index{0}; + for (std::size_t component = 0; + component < Samples::component_count; ++component) + for (std::size_t sample = 0; sample < actual.npts(); ++sample, ++index) + check(Samples::get(actual, component, sample) == + expected.samples[index], + context + ": sample changed at component " + + std::to_string(component) + ", index " + + std::to_string(sample)); +} + +template +void combine(Waveform &lhs, const Waveform &rhs, const bool add) { + if (add) + lhs += rhs; + else + lhs -= rhs; +} + +template +void verify_valid(const bool add, const double rhs_t0, const double rhs_dt, + const int expected_offset, const std::string &context, + const double lhs_dt = 1.0) { + Waveform lhs = make_waveform(5, 0.0, lhs_dt, 100.0); + const Waveform rhs = make_waveform(4, rhs_t0, rhs_dt, 10.0); + Waveform expected(lhs); + for (std::size_t rhs_index = 0; rhs_index < rhs.npts(); ++rhs_index) { + const int lhs_index = expected_offset + static_cast(rhs_index); + if (lhs_index < 0 || lhs_index >= static_cast(lhs.npts())) + continue; + for (std::size_t component = 0; + component < Samples::component_count; ++component) { + const double value = Samples::get(lhs, component, lhs_index) + + (add ? 1.0 : -1.0) * Samples::get( + rhs, component, rhs_index); + Samples::set(expected, component, lhs_index, value); + } + } + + combine(lhs, rhs, add); + check_state(lhs, snapshot(expected), context); +} + +template +void verify_rejected(const bool add, const double rhs_t0, const double rhs_dt, + const std::string &context, const double lhs_dt = 1.0) { + Waveform lhs = make_waveform(5, 0.0, lhs_dt, 100.0); + Waveform rhs = make_waveform(4, rhs_t0, rhs_dt, 10.0); + rhs.elog.log_error("rhs", "must not merge on rejection", + ErrorSeverity::Complaint); + const auto before = snapshot(lhs); + bool threw_invalid{false}; + try { + combine(lhs, rhs, add); + } catch (const MsPASSError &error) { + threw_invalid = error.severity() == ErrorSeverity::Invalid; + } + check(threw_invalid, context + ": did not throw MsPASSError Invalid"); + check_state(lhs, before, context + " rejection atomicity"); +} + +template +void run_suite(const bool add, const std::string &type_name) { + const std::string operation = add ? "+=" : "-="; + const std::string prefix = type_name + " " + operation; + + verify_valid(add, 0.0, 1.0, 0, prefix + " equal grid"); + verify_valid(add, 2.0, 1.0, 2, prefix + " positive offset"); + verify_valid(add, -2.0, 1.0, -2, prefix + " negative offset"); + verify_valid(add, 5.0, 1.0, 5, prefix + " positive no overlap"); + verify_valid(add, -4.0, 1.0, -4, prefix + " negative no overlap"); + + constexpr double tolerance = 1.0e-6; + const double positive_offset_at_tolerance = 1.0 + tolerance; + const double positive_offset_beyond = std::nextafter( + positive_offset_at_tolerance, std::numeric_limits::infinity()); + const double negative_offset_at_tolerance = -1.0 - tolerance; + const double negative_offset_beyond = std::nextafter( + negative_offset_at_tolerance, -std::numeric_limits::infinity()); + check(std::abs(positive_offset_at_tolerance - + std::round(positive_offset_at_tolerance)) <= tolerance, + "positive offset boundary fixture is invalid"); + check(std::abs(positive_offset_beyond - std::round(positive_offset_beyond)) > + tolerance, + "positive offset rejection fixture is invalid"); + check(std::abs(negative_offset_at_tolerance - + std::round(negative_offset_at_tolerance)) <= tolerance, + "negative offset boundary fixture is invalid"); + check(std::abs(negative_offset_beyond - std::round(negative_offset_beyond)) > + tolerance, + "negative offset rejection fixture is invalid"); + verify_valid(add, positive_offset_at_tolerance, 1.0, 1, + prefix + " positive offset at tolerance"); + verify_valid(add, negative_offset_at_tolerance, 1.0, -1, + prefix + " negative offset at tolerance"); + verify_rejected(add, positive_offset_beyond, 1.0, + prefix + " positive offset beyond tolerance"); + verify_rejected(add, negative_offset_beyond, 1.0, + prefix + " negative offset beyond tolerance"); + + constexpr double lhs_dt = 1.0e6; + constexpr double dt_at_tolerance = lhs_dt - 1.0; + const double dt_beyond = std::nextafter(dt_at_tolerance, 0.0); + check(std::abs(lhs_dt - dt_at_tolerance) == + tolerance * std::max(std::abs(lhs_dt), std::abs(dt_at_tolerance)), + "dt boundary fixture is invalid"); + check(std::abs(lhs_dt - dt_beyond) > + tolerance * std::max(std::abs(lhs_dt), std::abs(dt_beyond)), + "dt rejection fixture is invalid"); + verify_valid(add, 0.0, dt_at_tolerance, 0, + prefix + " dt at relative tolerance", lhs_dt); + verify_rejected(add, 0.0, dt_beyond, + prefix + " dt beyond relative tolerance", lhs_dt); +} +} // namespace + +int main() { + run_suite(true, "TimeSeries"); + run_suite(false, "TimeSeries"); + run_suite(true, "Seismogram"); + run_suite(false, "Seismogram"); + return failure_count == 0 ? 0 : 1; +} diff --git a/python/tests/test_waveform_arithmetic.py b/python/tests/test_waveform_arithmetic.py new file mode 100644 index 0000000000..611014a5b9 --- /dev/null +++ b/python/tests/test_waveform_arithmetic.py @@ -0,0 +1,139 @@ +import math + +import numpy as np +import pytest + +from mspasspy.ccore.seismic import Seismogram, TimeReferenceType, TimeSeries +from mspasspy.ccore.utility import ErrorSeverity, MsPASSError + + +def _make_waveform(waveform_type, sample_count, t0, dt, base): + waveform = waveform_type(sample_count) + waveform.t0 = t0 + waveform.dt = dt + waveform.tref = TimeReferenceType.Relative + waveform.set_live() + waveform["sentinel"] = "unchanged" + if waveform_type is TimeSeries: + for sample in range(sample_count): + waveform.data[sample] = base + sample + else: + for component in range(3): + waveform.data[component, :] = ( + base + 10.0 * component + np.arange(sample_count) + ) + return waveform + + +def _snapshot(waveform): + return { + "npts": waveform.npts, + "t0": waveform.t0, + "dt": waveform.dt, + "tref": waveform.tref, + "live": waveform.live, + "error_count": waveform.elog.size(), + "sentinel": waveform["sentinel"], + "data": np.array(waveform.data, copy=True), + } + + +def _assert_state(waveform, expected): + assert waveform.npts == expected["npts"] + assert waveform.t0 == expected["t0"] + assert waveform.dt == expected["dt"] + assert waveform.tref == expected["tref"] + assert waveform.live == expected["live"] + assert waveform.elog.size() == expected["error_count"] + assert waveform["sentinel"] == expected["sentinel"] + np.testing.assert_array_equal(waveform.data, expected["data"]) + + +def _combine(lhs, rhs, operation): + if operation == "add": + lhs += rhs + else: + lhs -= rhs + + +def _verify_valid(waveform_type, operation, rhs_t0, rhs_dt, offset, lhs_dt=1.0): + lhs = _make_waveform(waveform_type, 5, 0.0, lhs_dt, 100.0) + rhs = _make_waveform(waveform_type, 4, rhs_t0, rhs_dt, 10.0) + expected = _snapshot(lhs) + sign = 1.0 if operation == "add" else -1.0 + for rhs_index in range(rhs.npts): + lhs_index = offset + rhs_index + if 0 <= lhs_index < lhs.npts: + if waveform_type is TimeSeries: + expected["data"][lhs_index] += sign * rhs.data[rhs_index] + else: + expected["data"][:, lhs_index] += sign * rhs.data[:, rhs_index] + + _combine(lhs, rhs, operation) + _assert_state(lhs, expected) + + +def _verify_rejected(waveform_type, operation, rhs_t0, rhs_dt, lhs_dt=1.0): + lhs = _make_waveform(waveform_type, 5, 0.0, lhs_dt, 100.0) + rhs = _make_waveform(waveform_type, 4, rhs_t0, rhs_dt, 10.0) + rhs.elog.log_error("rhs", "must not merge on rejection", ErrorSeverity.Complaint) + before = _snapshot(lhs) + + with pytest.raises(MsPASSError) as exc_info: + _combine(lhs, rhs, operation) + + assert exc_info.value.severity == ErrorSeverity.Invalid + _assert_state(lhs, before) + + +@pytest.mark.parametrize("waveform_type", [TimeSeries, Seismogram]) +@pytest.mark.parametrize("operation", ["add", "subtract"]) +def test_waveform_arithmetic_grid_contract(waveform_type, operation): + valid_cases = [ + (0.0, 1.0, 0), + (2.0, 1.0, 2), + (-2.0, 1.0, -2), + (5.0, 1.0, 5), + (-4.0, 1.0, -4), + ] + + tolerance = 1.0e-6 + positive_offset_at_tolerance = 1.0 + tolerance + negative_offset_at_tolerance = -1.0 - tolerance + assert ( + abs(positive_offset_at_tolerance - round(positive_offset_at_tolerance)) + <= tolerance + ) + assert ( + abs(negative_offset_at_tolerance - round(negative_offset_at_tolerance)) + <= tolerance + ) + valid_cases.extend( + [ + (positive_offset_at_tolerance, 1.0, 1), + (negative_offset_at_tolerance, 1.0, -1), + ] + ) + + for rhs_t0, rhs_dt, offset in valid_cases: + _verify_valid(waveform_type, operation, rhs_t0, rhs_dt, offset) + + lhs_dt = 1.0e6 + dt_at_tolerance = lhs_dt - 1.0 + dt_beyond = math.nextafter(dt_at_tolerance, 0.0) + assert abs(lhs_dt - dt_at_tolerance) == tolerance * max( + abs(lhs_dt), abs(dt_at_tolerance) + ) + assert abs(lhs_dt - dt_beyond) > tolerance * max(abs(lhs_dt), abs(dt_beyond)) + _verify_valid(waveform_type, operation, 0.0, dt_at_tolerance, 0, lhs_dt=lhs_dt) + + positive_offset_beyond = math.nextafter(positive_offset_at_tolerance, math.inf) + negative_offset_beyond = math.nextafter(negative_offset_at_tolerance, -math.inf) + assert abs(positive_offset_beyond - round(positive_offset_beyond)) > tolerance + assert abs(negative_offset_beyond - round(negative_offset_beyond)) > tolerance + for rhs_t0, rhs_dt in [ + (positive_offset_beyond, 1.0), + (negative_offset_beyond, 1.0), + ]: + _verify_rejected(waveform_type, operation, rhs_t0, rhs_dt) + _verify_rejected(waveform_type, operation, 0.0, dt_beyond, lhs_dt=lhs_dt) From 2493c2ce50b00f9b63a128c1e2c8f468564bb4c6 Mon Sep 17 00:00:00 2001 From: wangyinz Date: Thu, 13 Aug 2026 02:59:48 -0500 Subject: [PATCH 2/4] Align waveform arithmetic callers to sample grids --- python/mspasspy/algorithms/MCXcorStacking.py | 65 +++++++++------ .../tests/algorithms/test_MCXcorStacking.py | 80 +++++++++++++++++++ .../test_manager_no_spark_dask.py | 4 +- .../global_history/test_manager_spark_dask.py | 5 +- python/tests/helper.py | 18 +++-- python/tests/test_reduce.py | 29 +++++-- 6 files changed, 161 insertions(+), 40 deletions(-) diff --git a/python/mspasspy/algorithms/MCXcorStacking.py b/python/mspasspy/algorithms/MCXcorStacking.py index d526215b02..5d14dd0519 100644 --- a/python/mspasspy/algorithms/MCXcorStacking.py +++ b/python/mspasspy/algorithms/MCXcorStacking.py @@ -1010,14 +1010,14 @@ def _dbxcor_stacker( for i in range(maxiterations): # this requires stack0 not be altered in this function wts = dbxcor_weights(ensemble, stack0, residual_norm_floor=residual_norm_floor) - # newstack = np.zeros(N) - # this is just a fast way to initalize to 0s + # this is just a fast way to initialize to zeros stack.set_npts(N) sumwts = 0.0 for j in range(M): if ensemble.member[j].live and wts[j] > 0.0: d = TimeSeries(ensemble.member[j]) d *= wts[j] + _snap_start_time_to_grid(d, stack) stack += d sumwts += wts[j] if sumwts > 0.0: @@ -1031,7 +1031,9 @@ def _dbxcor_stacker( # order may matter here. In this case delta becomes a numpy # array which is cleaner in this context # delta = newstack - stack.data - delta = stack - stack0 + aligned_stack0 = TimeSeries(stack0) + _snap_start_time_to_grid(aligned_stack0, stack) + delta = stack - aligned_stack0 relative_delta = np.linalg.norm(delta.data) / np.linalg.norm(stack0.data) # normalize by sample size or there is a window length dependency relative_delta /= N @@ -1117,15 +1119,20 @@ def beam_align(ensemble, beam, window=None, time_shift_limit=10.0): timelag = _xcor_shift(d, beam) # apply a ceiling/floor to allowed time shift via # the time_shift_limit arg + lag_was_limited = False if timelag > time_shift_limit: timelag = time_shift_limit + lag_was_limited = True elif timelag < (-time_shift_limit): timelag = -time_shift_limit + lag_was_limited = True # We MUST use this method instead of dithering t0 to keep # absolute time right. This will fail if the inputs were # not shifted from UTC times # also note a +lag requires a - shift ensemble.member[i].shift(timelag) + if not lag_was_limited: + _snap_start_time_to_grid(ensemble.member[i], beam) return ensemble @@ -1600,7 +1607,9 @@ def align_and_stack( timespan_method="ensemble_median", stack_md=Metadata(rbeam0), ) - delta_rbeam = rbeam - rbeam0 + aligned_rbeam0 = TimeSeries(rbeam0) + _snap_start_time_to_grid(aligned_rbeam0, rbeam) + delta_rbeam = rbeam - aligned_rbeam0 nrm_delta = np.linalg.norm(delta_rbeam.data) if nrm_rbeam == 0.0: if nrm_delta == 0.0: @@ -1645,6 +1654,10 @@ def align_and_stack( # this method alters the t0 values of the ensemble members # when used plots will tend to be aligned with 0 relative time ensemble.member[i].shift(tshift_mean) + output_reference = next(d for d in ensemble.member if d.live) + for datum in ensemble.member: + if datum.live: + _snap_start_time_to_grid(datum, output_reference, tolerance=1.0e-6) for i in range(len(ensemble.member)): if ensemble.member[i].live: # in this context it0_key should always be defined @@ -1654,24 +1667,14 @@ def align_and_stack( tshift = ensemble.member[i].t0 - initial_starttime ensemble.member[i].put_double(time_shift_key, tshift) if output_stack_window: - # this will clone the beam trace metadata automatically - # using pad option assures t0 will be output_stack_window.start - # and npts is consistent with window requested - output_stack = WindowData( - beam, - output_stack_window.start, - output_stack_window.end, - short_segment_handling="pad", - ) - # this is an obscure but fast way to initialize the data vector to all 0s - output_stack.set_npts(output_stack.npts) + requested_output_window = TimeWindow(output_stack_window) else: - # also clones beam metadata but in this case we get the size from the ensemble time span - output_stack = TimeSeries(beam) - output_stack_window = TimeWindow(ensemble_timespan) - output_stack.set_t0(output_stack_window.start) - npts = int((output_stack_window.end - output_stack_window.start) / beam.dt) + 1 - output_stack.set_npts(npts) + requested_output_window = ensemble_time_range(ensemble, metric="median") + output_start_index = output_reference.sample_number(requested_output_window.start) + output_end_index = output_reference.sample_number(requested_output_window.end) + output_stack = TimeSeries(beam) + output_stack.set_t0(output_reference.time(output_start_index)) + output_stack.set_npts(output_end_index - output_start_index + 1) if robust_stack_method == "dbxcor": # We need to post the final weights to all live members wts = dbxcor_weights(rens, rbeam, residual_norm_floor=residual_norm_floor) @@ -2452,6 +2455,20 @@ def _xcor_shift(ts, beam): return lagtime +def _snap_start_time_to_grid(datum, reference, tolerance=None): + """Snap ``datum.t0`` to ``reference`` without changing any samples. + + With ``tolerance`` omitted this applies the explicit nearest-sample + indexing used by the stacking operations. A tolerance limits the snap to + floating-point residue on an already common grid, so a deliberately + clipped cross-correlation shift remains visible to later iterations. + """ + sample_number = reference.sample_number(datum.t0) + grid_time = reference.time(sample_number) + if tolerance is None or abs((datum.t0 - grid_time) / reference.dt) <= tolerance: + datum.t0 = grid_time + + def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries: """ Internal method used to update the beam signal used for cross correlation @@ -2503,13 +2520,11 @@ def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries: xcorens.member[i], stime, etime, short_segment_handling="pad" ) d *= wt - # TimeSeries::operator+= handles time correctly so indexing is not needed here + _snap_start_time_to_grid(d, beam) beam += d sumwt += wt if sumwt > 0.0: - # TimeSeries does not have operator /= defined but it does have *= defined - scale = 1.0 / sumwt - beam *= scale + beam *= 1.0 / sumwt else: beam.elog.log_error( "_update_xcor_beam", diff --git a/python/tests/algorithms/test_MCXcorStacking.py b/python/tests/algorithms/test_MCXcorStacking.py index b8e412d230..9b9ad8e0e2 100755 --- a/python/tests/algorithms/test_MCXcorStacking.py +++ b/python/tests/algorithms/test_MCXcorStacking.py @@ -24,6 +24,7 @@ ) from mspasspy.ccore.algorithms.basic import TimeWindow from mspasspy.ccore.algorithms.amplitudes import MADAmplitude +from mspasspy.ccore.utility import ErrorSeverity, MsPASSError from mspasspy.algorithms.window import WindowData from mspasspy.algorithms.basic import ExtractComponent from mspasspy.algorithms.window import scale @@ -273,6 +274,10 @@ def test_align_and_stack(): robust_stack_window=rwin, correlation_window=xcorwin, ) + for datum in eo.member: + if datum.live: + grid_offset = (datum.t0 - beam.t0) / beam.dt + assert abs(grid_offset - round(grid_offset)) <= 1.0e-6 # test if this is right validate_lags(eo, lag_in_samples) # Repeat with one signal replaced by random noise @@ -391,6 +396,81 @@ def test_align_and_stack(): validate_lags(eo, lag_in_samples) +def test_align_and_stack_preserves_limited_offsets_and_strict_arithmetic(): + """A clipped lag remains visible while stacking uses explicit grid indices.""" + [_, ensemble, _] = load_test_data() + for datum in ensemble.member: + datum.ator(1000.0) + noisy_member = 5 + ensemble.member[noisy_member].data = DoubleVector( + np.random.default_rng(20260813).normal(size=ensemble.member[noisy_member].npts) + ) + samples_before = [np.array(datum.data, copy=True) for datum in ensemble.member] + + output, beam = align_and_stack( + ensemble, + ensemble.member[0], + window_beam=True, + robust_stack_window=TimeWindow(-1.0, 3.0), + correlation_window=TimeWindow(-2.0, 10.0), + time_shift_limit=0.03, + ) + + assert output.live + assert beam.live + for before, datum in zip(samples_before, output.member): + np.testing.assert_array_equal(datum.data, before) + + grid_residuals = [ + abs((datum.t0 - beam.t0) / beam.dt - round((datum.t0 - beam.t0) / beam.dt)) + for datum in output.member + if datum.live + ] + assert max(grid_residuals) > 1.0e-3 + + misaligned_member = next( + datum + for datum in output.member + if datum.live + and abs((datum.t0 - beam.t0) / beam.dt - round((datum.t0 - beam.t0) / beam.dt)) + > 1.0e-6 + ) + arithmetic_lhs = TimeSeries(beam) + lhs_before = { + "t0": arithmetic_lhs.t0, + "dt": arithmetic_lhs.dt, + "npts": arithmetic_lhs.npts, + "live": arithmetic_lhs.live, + "elog_size": arithmetic_lhs.elog.size(), + "data": np.array(arithmetic_lhs.data, copy=True), + } + with pytest.raises(MsPASSError) as exc_info: + arithmetic_lhs += misaligned_member + assert exc_info.value.severity == ErrorSeverity.Invalid + assert arithmetic_lhs.t0 == lhs_before["t0"] + assert arithmetic_lhs.dt == lhs_before["dt"] + assert arithmetic_lhs.npts == lhs_before["npts"] + assert arithmetic_lhs.live == lhs_before["live"] + assert arithmetic_lhs.elog.size() == lhs_before["elog_size"] + np.testing.assert_array_equal(arithmetic_lhs.data, lhs_before["data"]) + + expected_stack = np.zeros(beam.npts) + weight_sum = 0.0 + for datum in output.member: + if datum.live and datum["robust_stack_weight"] > 0.0: + weight = datum["robust_stack_weight"] + offset = beam.sample_number(datum.t0) + beam_begin = max(0, offset) + datum_begin = max(0, -offset) + count = min(beam.npts - beam_begin, datum.npts - datum_begin) + expected_stack[beam_begin : beam_begin + count] += ( + weight * np.asarray(datum.data)[datum_begin : datum_begin + count] + ) + weight_sum += weight + expected_stack /= weight_sum + np.testing.assert_allclose(beam.data, expected_stack, rtol=0.0, atol=1.0e-15) + + def test_align_and_stack_error_handlers(): """ Like above but only tests error handlers. Uses the diff --git a/python/tests/global_history/test_manager_no_spark_dask.py b/python/tests/global_history/test_manager_no_spark_dask.py index 40033f959e..d931377520 100644 --- a/python/tests/global_history/test_manager_no_spark_dask.py +++ b/python/tests/global_history/test_manager_no_spark_dask.py @@ -2,7 +2,7 @@ import sys sys.path.append("python/tests") -from helper import get_live_timeseries +from helper import get_live_timeseries, get_live_timeseries_list with mock.patch.dict(sys.modules, {"pyspark": None, "dask": None}): from mspasspy.reduce import stack @@ -39,7 +39,7 @@ def test_normal_map(): assert test_map_res[i].data == t[i].data + t[i].data def test_normal_reduce(): - t = [get_live_timeseries() for i in range(5)] + t = get_live_timeseries_list(5) s = t[0] for i in range(1, 5): s += t[i] diff --git a/python/tests/global_history/test_manager_spark_dask.py b/python/tests/global_history/test_manager_spark_dask.py index 195d56654b..8a501f470f 100644 --- a/python/tests/global_history/test_manager_spark_dask.py +++ b/python/tests/global_history/test_manager_spark_dask.py @@ -29,6 +29,7 @@ from helper import ( get_live_seismogram, get_live_timeseries, + get_live_timeseries_list, get_live_timeseries_ensemble, get_live_seismogram_ensemble, ) @@ -161,7 +162,7 @@ def test_normal_map(self): assert test_map_res[i].data == t[i].data + t[i].data def test_normal_reduce(self): - t = [get_live_timeseries() for i in range(5)] + t = get_live_timeseries_list(5) s = t[0] for i in range(1, 5): s += t[i] @@ -501,7 +502,7 @@ def test_mspass_reduce(self, spark_context): manager_db = Database(self.client, "test_manager") manager_db["history_global"].delete_many({}) - l = [get_live_timeseries() for i in range(5)] + l = get_live_timeseries_list(5) # test mspass_reduce for spark spark_res = spark_reduce(l, self.manager, spark_context) assert ( diff --git a/python/tests/helper.py b/python/tests/helper.py index 3ac65afc29..cf7eab8b9e 100644 --- a/python/tests/helper.py +++ b/python/tests/helper.py @@ -37,7 +37,7 @@ def get_live_seismogram(ts_size=255, sampling_rate=20.0): return seis -def get_live_timeseries(ts_size=255, sampling_rate=20.0): +def get_live_timeseries(ts_size=255, sampling_rate=20.0, start_time=None): ts = TimeSeries() ts.set_live() ts.dt = 1 / sampling_rate @@ -46,7 +46,9 @@ def get_live_timeseries(ts_size=255, sampling_rate=20.0): ts.put("npts", ts_size) ts.put("sampling_rate", sampling_rate) ts.tref = TimeReferenceType.UTC - ts.t0 = datetime.datetime.now().timestamp() + if start_time is None: + start_time = datetime.datetime.now().timestamp() + ts.t0 = start_time ts["delta"] = 0.1 ts["calib"] = 0.1 ts["site_id"] = bson.objectid.ObjectId() @@ -90,19 +92,23 @@ def get_live_seismogram_ensemble(n): return seis_e -def get_live_timeseries_ensemble(n): +def get_live_timeseries_ensemble(n, start_time=None): tse = TimeSeriesEnsemble() + if start_time is None: + start_time = datetime.datetime.now().timestamp() for i in range(n): - ts = get_live_timeseries() + ts = get_live_timeseries(start_time=start_time) tse.member.append(ts) tse.set_live() return tse -def get_live_timeseries_list(n, ts_size=255): +def get_live_timeseries_list(n, ts_size=255, start_time=None): ts_e = [] + if start_time is None: + start_time = datetime.datetime.now().timestamp() for i in range(n): - ts = get_live_timeseries(ts_size) + ts = get_live_timeseries(ts_size, start_time=start_time) ts_e.append(ts) return ts_e diff --git a/python/tests/test_reduce.py b/python/tests/test_reduce.py index c40a7f4c37..6f7f47eb3f 100644 --- a/python/tests/test_reduce.py +++ b/python/tests/test_reduce.py @@ -25,6 +25,7 @@ from helper import ( get_live_seismogram, get_live_timeseries, + get_live_timeseries_list, get_live_timeseries_ensemble, get_live_seismogram_ensemble, get_live_seismogram_list, @@ -84,14 +85,15 @@ def test_reduce_stack(): for i in range(3): assert np.isclose(seis1.data[i], res[i]).all() # fixme - ts1 = get_live_timeseries() - ts2 = get_live_timeseries() + start_time = 1234567890.0 + ts1 = get_live_timeseries(start_time=start_time) + ts2 = get_live_timeseries(start_time=start_time) ts1_cp = np.array(ts1.data) stack(ts1, ts2) assert np.isclose(ts1.data, (np.array(ts1_cp) + np.array(ts2.data))).all() - tse1 = get_live_timeseries_ensemble(2) - tse2 = get_live_timeseries_ensemble(2) + tse1 = get_live_timeseries_ensemble(2, start_time=start_time) + tse2 = get_live_timeseries_ensemble(2, start_time=start_time) tse1_cp = TimeSeriesEnsemble(tse1) stack(tse1, tse2) for i in range(2): @@ -112,6 +114,23 @@ def test_reduce_stack(): assert np.isclose(seis_e1.member[i].data[j], res[j]).all() # fixme +def test_timeseries_helpers_share_explicit_start_time(): + start_time = 0.0 + timeseries_list = get_live_timeseries_list(3, start_time=start_time) + timeseries_ensemble = get_live_timeseries_ensemble(3, start_time=start_time) + + assert [datum.t0 for datum in timeseries_list] == [start_time] * 3 + assert [datum.t0 for datum in timeseries_ensemble.member] == [start_time] * 3 + list_other_values = [datum.data[0] for datum in timeseries_list[1:]] + ensemble_other_values = [datum.data[0] for datum in timeseries_ensemble.member[1:]] + timeseries_list[0].data[0] = 1.0 + timeseries_ensemble.member[0].data[0] = 1.0 + assert [datum.data[0] for datum in timeseries_list[1:]] == list_other_values + assert [ + datum.data[0] for datum in timeseries_ensemble.member[1:] + ] == ensemble_other_values + + def test_reduce_stack_exception(): tse1 = get_live_timeseries_ensemble(2) tse2 = get_live_timeseries_ensemble(3) @@ -144,7 +163,7 @@ def spark_reduce(input, sc): def test_reduce_dask_spark(spark_context): - l = [get_live_timeseries() for i in range(5)] + l = get_live_timeseries_list(5) res = np.zeros(255) for i in range(5): for j in range(255): From 4973c2e8eb538aa168282ed3f7d1993dfeb7ea86 Mon Sep 17 00:00:00 2001 From: wangyinz Date: Fri, 14 Aug 2026 03:10:17 -0500 Subject: [PATCH 3/4] Keep waveform grid validation from rewriting MCXcor --- python/mspasspy/algorithms/MCXcorStacking.py | 65 ++++++--------- .../tests/algorithms/test_MCXcorStacking.py | 80 ------------------- 2 files changed, 25 insertions(+), 120 deletions(-) diff --git a/python/mspasspy/algorithms/MCXcorStacking.py b/python/mspasspy/algorithms/MCXcorStacking.py index 5d14dd0519..d526215b02 100644 --- a/python/mspasspy/algorithms/MCXcorStacking.py +++ b/python/mspasspy/algorithms/MCXcorStacking.py @@ -1010,14 +1010,14 @@ def _dbxcor_stacker( for i in range(maxiterations): # this requires stack0 not be altered in this function wts = dbxcor_weights(ensemble, stack0, residual_norm_floor=residual_norm_floor) - # this is just a fast way to initialize to zeros + # newstack = np.zeros(N) + # this is just a fast way to initalize to 0s stack.set_npts(N) sumwts = 0.0 for j in range(M): if ensemble.member[j].live and wts[j] > 0.0: d = TimeSeries(ensemble.member[j]) d *= wts[j] - _snap_start_time_to_grid(d, stack) stack += d sumwts += wts[j] if sumwts > 0.0: @@ -1031,9 +1031,7 @@ def _dbxcor_stacker( # order may matter here. In this case delta becomes a numpy # array which is cleaner in this context # delta = newstack - stack.data - aligned_stack0 = TimeSeries(stack0) - _snap_start_time_to_grid(aligned_stack0, stack) - delta = stack - aligned_stack0 + delta = stack - stack0 relative_delta = np.linalg.norm(delta.data) / np.linalg.norm(stack0.data) # normalize by sample size or there is a window length dependency relative_delta /= N @@ -1119,20 +1117,15 @@ def beam_align(ensemble, beam, window=None, time_shift_limit=10.0): timelag = _xcor_shift(d, beam) # apply a ceiling/floor to allowed time shift via # the time_shift_limit arg - lag_was_limited = False if timelag > time_shift_limit: timelag = time_shift_limit - lag_was_limited = True elif timelag < (-time_shift_limit): timelag = -time_shift_limit - lag_was_limited = True # We MUST use this method instead of dithering t0 to keep # absolute time right. This will fail if the inputs were # not shifted from UTC times # also note a +lag requires a - shift ensemble.member[i].shift(timelag) - if not lag_was_limited: - _snap_start_time_to_grid(ensemble.member[i], beam) return ensemble @@ -1607,9 +1600,7 @@ def align_and_stack( timespan_method="ensemble_median", stack_md=Metadata(rbeam0), ) - aligned_rbeam0 = TimeSeries(rbeam0) - _snap_start_time_to_grid(aligned_rbeam0, rbeam) - delta_rbeam = rbeam - aligned_rbeam0 + delta_rbeam = rbeam - rbeam0 nrm_delta = np.linalg.norm(delta_rbeam.data) if nrm_rbeam == 0.0: if nrm_delta == 0.0: @@ -1654,10 +1645,6 @@ def align_and_stack( # this method alters the t0 values of the ensemble members # when used plots will tend to be aligned with 0 relative time ensemble.member[i].shift(tshift_mean) - output_reference = next(d for d in ensemble.member if d.live) - for datum in ensemble.member: - if datum.live: - _snap_start_time_to_grid(datum, output_reference, tolerance=1.0e-6) for i in range(len(ensemble.member)): if ensemble.member[i].live: # in this context it0_key should always be defined @@ -1667,14 +1654,24 @@ def align_and_stack( tshift = ensemble.member[i].t0 - initial_starttime ensemble.member[i].put_double(time_shift_key, tshift) if output_stack_window: - requested_output_window = TimeWindow(output_stack_window) + # this will clone the beam trace metadata automatically + # using pad option assures t0 will be output_stack_window.start + # and npts is consistent with window requested + output_stack = WindowData( + beam, + output_stack_window.start, + output_stack_window.end, + short_segment_handling="pad", + ) + # this is an obscure but fast way to initialize the data vector to all 0s + output_stack.set_npts(output_stack.npts) else: - requested_output_window = ensemble_time_range(ensemble, metric="median") - output_start_index = output_reference.sample_number(requested_output_window.start) - output_end_index = output_reference.sample_number(requested_output_window.end) - output_stack = TimeSeries(beam) - output_stack.set_t0(output_reference.time(output_start_index)) - output_stack.set_npts(output_end_index - output_start_index + 1) + # also clones beam metadata but in this case we get the size from the ensemble time span + output_stack = TimeSeries(beam) + output_stack_window = TimeWindow(ensemble_timespan) + output_stack.set_t0(output_stack_window.start) + npts = int((output_stack_window.end - output_stack_window.start) / beam.dt) + 1 + output_stack.set_npts(npts) if robust_stack_method == "dbxcor": # We need to post the final weights to all live members wts = dbxcor_weights(rens, rbeam, residual_norm_floor=residual_norm_floor) @@ -2455,20 +2452,6 @@ def _xcor_shift(ts, beam): return lagtime -def _snap_start_time_to_grid(datum, reference, tolerance=None): - """Snap ``datum.t0`` to ``reference`` without changing any samples. - - With ``tolerance`` omitted this applies the explicit nearest-sample - indexing used by the stacking operations. A tolerance limits the snap to - floating-point residue on an already common grid, so a deliberately - clipped cross-correlation shift remains visible to later iterations. - """ - sample_number = reference.sample_number(datum.t0) - grid_time = reference.time(sample_number) - if tolerance is None or abs((datum.t0 - grid_time) / reference.dt) <= tolerance: - datum.t0 = grid_time - - def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries: """ Internal method used to update the beam signal used for cross correlation @@ -2520,11 +2503,13 @@ def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries: xcorens.member[i], stime, etime, short_segment_handling="pad" ) d *= wt - _snap_start_time_to_grid(d, beam) + # TimeSeries::operator+= handles time correctly so indexing is not needed here beam += d sumwt += wt if sumwt > 0.0: - beam *= 1.0 / sumwt + # TimeSeries does not have operator /= defined but it does have *= defined + scale = 1.0 / sumwt + beam *= scale else: beam.elog.log_error( "_update_xcor_beam", diff --git a/python/tests/algorithms/test_MCXcorStacking.py b/python/tests/algorithms/test_MCXcorStacking.py index 9b9ad8e0e2..b8e412d230 100755 --- a/python/tests/algorithms/test_MCXcorStacking.py +++ b/python/tests/algorithms/test_MCXcorStacking.py @@ -24,7 +24,6 @@ ) from mspasspy.ccore.algorithms.basic import TimeWindow from mspasspy.ccore.algorithms.amplitudes import MADAmplitude -from mspasspy.ccore.utility import ErrorSeverity, MsPASSError from mspasspy.algorithms.window import WindowData from mspasspy.algorithms.basic import ExtractComponent from mspasspy.algorithms.window import scale @@ -274,10 +273,6 @@ def test_align_and_stack(): robust_stack_window=rwin, correlation_window=xcorwin, ) - for datum in eo.member: - if datum.live: - grid_offset = (datum.t0 - beam.t0) / beam.dt - assert abs(grid_offset - round(grid_offset)) <= 1.0e-6 # test if this is right validate_lags(eo, lag_in_samples) # Repeat with one signal replaced by random noise @@ -396,81 +391,6 @@ def test_align_and_stack(): validate_lags(eo, lag_in_samples) -def test_align_and_stack_preserves_limited_offsets_and_strict_arithmetic(): - """A clipped lag remains visible while stacking uses explicit grid indices.""" - [_, ensemble, _] = load_test_data() - for datum in ensemble.member: - datum.ator(1000.0) - noisy_member = 5 - ensemble.member[noisy_member].data = DoubleVector( - np.random.default_rng(20260813).normal(size=ensemble.member[noisy_member].npts) - ) - samples_before = [np.array(datum.data, copy=True) for datum in ensemble.member] - - output, beam = align_and_stack( - ensemble, - ensemble.member[0], - window_beam=True, - robust_stack_window=TimeWindow(-1.0, 3.0), - correlation_window=TimeWindow(-2.0, 10.0), - time_shift_limit=0.03, - ) - - assert output.live - assert beam.live - for before, datum in zip(samples_before, output.member): - np.testing.assert_array_equal(datum.data, before) - - grid_residuals = [ - abs((datum.t0 - beam.t0) / beam.dt - round((datum.t0 - beam.t0) / beam.dt)) - for datum in output.member - if datum.live - ] - assert max(grid_residuals) > 1.0e-3 - - misaligned_member = next( - datum - for datum in output.member - if datum.live - and abs((datum.t0 - beam.t0) / beam.dt - round((datum.t0 - beam.t0) / beam.dt)) - > 1.0e-6 - ) - arithmetic_lhs = TimeSeries(beam) - lhs_before = { - "t0": arithmetic_lhs.t0, - "dt": arithmetic_lhs.dt, - "npts": arithmetic_lhs.npts, - "live": arithmetic_lhs.live, - "elog_size": arithmetic_lhs.elog.size(), - "data": np.array(arithmetic_lhs.data, copy=True), - } - with pytest.raises(MsPASSError) as exc_info: - arithmetic_lhs += misaligned_member - assert exc_info.value.severity == ErrorSeverity.Invalid - assert arithmetic_lhs.t0 == lhs_before["t0"] - assert arithmetic_lhs.dt == lhs_before["dt"] - assert arithmetic_lhs.npts == lhs_before["npts"] - assert arithmetic_lhs.live == lhs_before["live"] - assert arithmetic_lhs.elog.size() == lhs_before["elog_size"] - np.testing.assert_array_equal(arithmetic_lhs.data, lhs_before["data"]) - - expected_stack = np.zeros(beam.npts) - weight_sum = 0.0 - for datum in output.member: - if datum.live and datum["robust_stack_weight"] > 0.0: - weight = datum["robust_stack_weight"] - offset = beam.sample_number(datum.t0) - beam_begin = max(0, offset) - datum_begin = max(0, -offset) - count = min(beam.npts - beam_begin, datum.npts - datum_begin) - expected_stack[beam_begin : beam_begin + count] += ( - weight * np.asarray(datum.data)[datum_begin : datum_begin + count] - ) - weight_sum += weight - expected_stack /= weight_sum - np.testing.assert_allclose(beam.data, expected_stack, rtol=0.0, atol=1.0e-15) - - def test_align_and_stack_error_handlers(): """ Like above but only tests error handlers. Uses the From 9c267f37e48514404ef82ea8c8fd7c12e358218b Mon Sep 17 00:00:00 2001 From: wangyinz Date: Fri, 21 Aug 2026 15:32:44 -0500 Subject: [PATCH 4/4] Align MCXcor shifts to the nearest sample grid --- python/mspasspy/algorithms/MCXcorStacking.py | 21 +++++++++++++++--- .../tests/algorithms/test_MCXcorStacking.py | 22 +++++++++++++++++++ 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/python/mspasspy/algorithms/MCXcorStacking.py b/python/mspasspy/algorithms/MCXcorStacking.py index d526215b02..11c0e9d7df 100644 --- a/python/mspasspy/algorithms/MCXcorStacking.py +++ b/python/mspasspy/algorithms/MCXcorStacking.py @@ -1084,8 +1084,10 @@ def beam_align(ensemble, beam, window=None, time_shift_limit=10.0): because the value can to sent to a C++ method that it type sensitive) :return: the input ensemble with its members time shifted in place to align - with the time base of beam. If a window is defined it is used only for - correlation and is not applied to the ensemble members. + with the time base of beam. Each shifted start time is placed on the + nearest sample of the beam grid without changing any sample values. + If a window is defined it is used only for correlation and is not + applied to the ensemble members. """ # this may not be necessary for internal use but if used @@ -1126,6 +1128,7 @@ def beam_align(ensemble, beam, window=None, time_shift_limit=10.0): # not shifted from UTC times # also note a +lag requires a - shift ensemble.member[i].shift(timelag) + _snap_start_time_to_grid(ensemble.member[i], beam) return ensemble @@ -1645,6 +1648,7 @@ def align_and_stack( # this method alters the t0 values of the ensemble members # when used plots will tend to be aligned with 0 relative time ensemble.member[i].shift(tshift_mean) + _snap_start_time_to_grid(ensemble.member[i], beam) for i in range(len(ensemble.member)): if ensemble.member[i].live: # in this context it0_key should always be defined @@ -2452,6 +2456,17 @@ def _xcor_shift(ts, beam): return lagtime +def _snap_start_time_to_grid(datum, reference): + """Move ``datum.t0`` to the nearest sample on ``reference``'s grid. + + This is an explicit MCXcor alignment operation: it changes only the time + assigned to sample zero and never changes the sample values, ``dt``, or + ``npts``. General waveform arithmetic remains strict and rejects operands + that have not first been placed on a common grid. + """ + datum.t0 = reference.time(reference.sample_number(datum.t0)) + + def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries: """ Internal method used to update the beam signal used for cross correlation @@ -2503,7 +2518,7 @@ def _update_xcor_beam(xcorens, beam0, robust_stack_method, wts) -> TimeSeries: xcorens.member[i], stime, etime, short_segment_handling="pad" ) d *= wt - # TimeSeries::operator+= handles time correctly so indexing is not needed here + _snap_start_time_to_grid(d, beam) beam += d sumwt += wt if sumwt > 0.0: diff --git a/python/tests/algorithms/test_MCXcorStacking.py b/python/tests/algorithms/test_MCXcorStacking.py index b8e412d230..6864271a9b 100755 --- a/python/tests/algorithms/test_MCXcorStacking.py +++ b/python/tests/algorithms/test_MCXcorStacking.py @@ -40,6 +40,7 @@ _set_phases, regularize_sampling, ensemble_time_range, + _snap_start_time_to_grid, ) @@ -391,6 +392,27 @@ def test_align_and_stack(): validate_lags(eo, lag_in_samples) +def test_mcxcor_nearest_grid_alignment_preserves_samples(): + reference = TimeSeries(8) + reference.dt = 0.05 + reference.t0 = -1.0 + reference.set_live() + + datum = TimeSeries(6) + datum.dt = reference.dt + datum.t0 = reference.t0 + 2.49 * reference.dt + datum.set_live() + datum.data = DoubleVector(np.arange(datum.npts, dtype=float)) + samples_before = np.array(datum.data, copy=True) + + _snap_start_time_to_grid(datum, reference) + + assert datum.t0 == reference.time(2) + assert datum.dt == reference.dt + assert datum.npts == len(samples_before) + np.testing.assert_array_equal(datum.data, samples_before) + + def test_align_and_stack_error_handlers(): """ Like above but only tests error handlers. Uses the