Skip to content

_align_two_time_arrays() adds an extra sample for some identical inputs, and always starts at dt instead of the records' start time #96

Open

Description

@jsh9

Found while reviewing #74. Pre-existing on main (b764057); #74 keeps the same behavior.

Problem

_align_two_time_arrays() builds the common time array as np.linspace(dt, tmax, num=n_time) with n_time = int(np.ceil(max(tmax1, tmax2) / dt)):

dt1 = t1[1] - t1[0]
dt2 = t2[1] - t2[0]
tmax1 = t1[-1]
tmax2 = t2[-1]
dt = min(dt1, dt2) # use the smaller time interval of the two time array
n_time = int(np.ceil(max(tmax1, tmax2) / dt))
tmax = dt * n_time # use the larger end time of the two as the end time
t_output = np.linspace(dt, tmax, num=n_time)
return t_output # noqa: RET504

This has two problems:

  1. One extra sample. max(tmax1, tmax2) / dt is often slightly above a whole number because of floating-point rounding (for example 0.30000000000000004 / 0.1). np.ceil() then adds one more sample, past the end of both inputs. For t = np.arange(1, n + 1) * 0.01 with n from 2 to 5000, aligning t with itself gives the wrong length in 342 of 4999 cases.
  2. The start time is ignored. The output always starts at dt. A record that starts at t = 0 (like tests/files/sample_accel.txt) loses its first sample. A record that starts later gets extra samples before its start.

compare_two_accel(), which Ground_Motion.compare() calls, interpolates both records onto this array with np.interp(). Outside a record, np.interp() repeats the record's first or last value, so the extra samples are filled with copies of those values:

t_ = _align_two_time_arrays(t_in, t_out)
a_in_ = np.interp(t_, t_in, a_in)
a_out_ = np.interp(t_, t_out, a_out)

fig, ax = sr.compare_two_accel(
accel_in,
accel_out,
smooth=smooth,
input_accel_label=input_accel_label,
output_accel_label=output_accel_label,
amplification_ylabel=amp_ylabel,
phase_shift_ylabel=phs_ylabel,
)

import numpy as np
import PySeismoSoil.helper_site_response as sr

# 1. One extra sample for two identical arrays
t = np.arange(1, 4) * 0.1  # [0.1, 0.2, 0.30000000000000004]
print(sr._align_two_time_arrays(t, t))  # [0.1 0.2 0.3 0.4]  (expected [0.1 0.2 0.3])

n_wrong = sum(
    len(sr._align_two_time_arrays(tt, tt)) != len(tt)
    for tt in (np.arange(1, n + 1) * 0.01 for n in range(2, 5001))
)
print(n_wrong)  # 342  (expected 0)

# 2. The output always starts at dt
t0 = np.array([0.0, 0.01, 0.02, 0.03])
print(sr._align_two_time_arrays(t0, t0))  # [0.01 0.02 0.03]  (expected [0.   0.01 0.02 0.03])

t5 = np.array([5.0, 5.5, 6.0, 6.5])
t_ = sr._align_two_time_arrays(t5, t5)
print(len(t_), t_[:3])  # 13 [0.5 1.  1.5]  (expected 4 [5.  5.5 6. ])

# What compare_two_accel() does next with these times
a5 = np.array([0.3, -1.0, 2.0, 0.7])
print(np.interp(t_, t5, a5))
# [ 0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3  0.3 -1.   2.   0.7]  (expected [ 0.3 -1.   2.   0.7])

The existing tests test_align_two_time_arrays__normal_case__not_identical_arrays and test_align_two_time_arrays__normal_case__identical_arrays on main (in #74: test_align_two_time_arrays__normal_case[not_identical_arrays] and [identical_arrays]) pass, because their arrays start at dt and their end times divide evenly.

Fix

  1. Start the output at the earlier start time, and count the samples with a small tolerance:

    dt = min(dt1, dt2)  # use the smaller time interval of the two time array
    t_start = min(t1[0], t2[0])
    t_end = max(tmax1, tmax2)
    n_time = int(np.ceil((t_end - t_start) / dt - 1e-6)) + 1
    return t_start + dt * np.arange(n_time)

    With only this change on main, the repro above gives the expected values (n_wrong is 0), the two existing cases still match their benchmarks, and all 240 tests pass.

  2. Add cases to the test (test_align_two_time_arrays__normal_case in Combine near-duplicate test functions into parametrized tests #74): aligning np.arange(1, 4) * 0.1, [0.0, 0.01, 0.02, 0.03] and [5.0, 5.5, 6.0, 6.5] with themselves returns the same array.

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions