version 1.1.2

* More accuracy in `calc_surface_energy` - now interpolates between time steps. More tests added.

HISTORY.rst

@@ -2,6 +2,12 @@
 History
 =======
 
+1.1.2 (2019-10-31)
+-------------------
+
+* More accuracy in `calc_surface_energy` - now interpolates between time steps. More tests added.
+
+
 1.1.1 (2019-10-29)
 -------------------
 

eqsig/__about__.py

@@ -1,4 +1,4 @@
 __project__ = "eqsig"
 __author__ = "Maxim Millen"
-__version__ = "1.1.1"
+__version__ = "1.1.2"
 __license__ = "MIT"

eqsig/surface.py

@@ -30,15 +30,14 @@ def trim_to_length(values, npts, surf2depth_travel_times, dt, trim=False, start=
     if start:  # Trim front
         sis = start_shift - surf_to_depth_shifts
         if not trim:
-            extras = np.max([sis, 0]) - np.min([np.min(2 * surf_to_depth_shifts), 0])
+            extras = np.max([np.max(sis), 0]) - np.min([np.min(2 * surf_to_depth_shifts), 0])
             npts = npts + extras  # plus the zero padding in front and back
     else:
         if trim:
             sis = np.zeros_like(surf_to_depth_shifts)
         else:  # no changes required
             return values
     outs = np.zeros((len(surf_to_depth_shifts), npts))
-    print('len_outs: ', outs.shape)
     for i in range(len(surf_to_depth_shifts)):
         if sis[i] < 0:
             outs[i] = values[i, -sis[i]: npts - sis[i]]
@@ -47,13 +46,13 @@ def trim_to_length(values, npts, surf2depth_travel_times, dt, trim=False, start=
     return outs
 
 
-def calc_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, stt=0.0, trim=False, start=False):
+def calc_surface_energy(asig, travel_times, nodal=True, up_red=1., down_red=1., stt=0.0, trim=False, start=False):
     """
     Calculates the energy at different travel times from a surface
 
     Parameters
     ----------
-    asig
+    asig: eqsig.AccSignal object
     travel_times: array_like
         Travel times from surface to depth of interest
     nodal: bool
@@ -73,9 +72,15 @@ def calc_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, st
     -------
 
     """
-    shifts = np.array(2 * travel_times / asig.dt, dtype=int)
-    up_wave = np.pad(asig.values, (0, np.max(shifts)), mode='constant', constant_values=0) * up_red  # 1d
-    down_waves = fn.put_array_in_2d_array(asig.values, shifts) * down_red
+    if not hasattr(up_red, '__len__'):
+        up_red = up_red * np.ones(len(travel_times))
+    if not hasattr(down_red, '__len__'):
+        down_red = down_red * np.ones(len(travel_times))
+    shifts = np.array(2 * travel_times / asig.dt)
+    max_shift = int(max(shifts))
+    up_wave = np.pad(asig.values, (0, max_shift), mode='constant', constant_values=0)[np.newaxis, :] * up_red[:, np.newaxis]  # 1d
+    dshifted = np.arange(asig.npts + max_shift)[np.newaxis, :] - shifts[:, np.newaxis]
+    down_waves = np.interp(dshifted, np.arange(asig.npts), asig.values, left=0) * down_red[:, np.newaxis]
     if nodal:
         acc_series = - down_waves + up_wave
     else:
@@ -85,8 +90,47 @@ def calc_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, st
     return trim_to_length(e, asig.npts, travel_times, asig.dt, trim=trim, start=start, s2s_travel_time=stt)
 
 
-def calc_cum_abs_surface_energy(asig, time_shifts, nodal=True, up_red=1, down_red=1, stt=0.0, trim=False, start=False):
-    energy = calc_surface_energy(asig, time_shifts, nodal=nodal, up_red=up_red, down_red=down_red, stt=stt,
+# def dep_calc_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, stt=0.0, trim=False, start=False):
+#     """
+#     Calculates the energy at different travel times from a surface
+#
+#     Parameters
+#     ----------
+#     asig
+#     travel_times: array_like
+#         Travel times from surface to depth of interest
+#     nodal: bool
+#         If true then surface is nodal (minima)
+#     up_red: float or array_like
+#         upward wave reduction factors
+#     down_red: float or array_like
+#         downward wave reduction factors
+#     trim: bool
+#         if true then forces array to be same length as values array
+#     start: bool
+#         if true then forces array to have the same start time as values array
+#     stt: float
+#         Travel time from input motion location to the surface
+#
+#     Returns
+#     -------
+#
+#     """
+#
+#     shifts = np.array(2 * travel_times / asig.dt, dtype=int)
+#     up_wave = np.pad(asig.values, (0, np.max(shifts)), mode='constant', constant_values=0) * up_red  # 1d
+#     down_waves = fn.put_array_in_2d_array(asig.values, shifts) * down_red
+#     if nodal:
+#         acc_series = - down_waves + up_wave
+#     else:
+#         acc_series = down_waves + up_wave
+#     velocity = scipy.integrate.cumtrapz(acc_series, dx=asig.dt, initial=0, axis=1)
+#     e = 0.5 * velocity * np.abs(velocity)
+#     return trim_to_length(e, asig.npts, travel_times, asig.dt, trim=trim, start=start, s2s_travel_time=stt)
+
+
+def calc_cum_abs_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, stt=0.0, trim=False, start=False):
+    energy = calc_surface_energy(asig, travel_times, nodal=nodal, up_red=up_red, down_red=down_red, stt=stt,
                                  trim=trim, start=start)
     diff = np.zeros_like(energy)
     diff[:, 1:] = np.diff(energy, axis=1)
@@ -98,6 +142,8 @@ def calc_cum_abs_surface_energy(asig, time_shifts, nodal=True, up_red=1, down_re
     b = np.sin(a)
     import eqsig
     accsig = eqsig.AccSignal(b, 0.1)
-    tshifts = np.array([0.1, 0.2])
-    ke = calc_cum_abs_surface_energy(accsig, tshifts, stt=0.3, nodal=True, up_red=1, down_red=1, trim=True)
+    t_times = np.array([0.1, 0.2])
+    ke = calc_cum_abs_surface_energy(accsig, t_times, stt=0.3, nodal=True, up_red=1, down_red=1, trim=True)
+
+    ke = calc_cum_abs_surface_energy(accsig, t_times, stt=0.3, nodal=True, up_red=tshifts, down_red=tshifts, trim=True)
 

tests/test_surface.py

@@ -76,6 +76,115 @@ def test_calc_cum_abs_surface_energy():
     assert np.isclose(cases[2][-1], expected_cases[2][-1])
 
 
+def test_calc_cum_abs_surface_energy_start_check_same_energy():
+    accsig = conftest.t_asig()
+    tshifts = np.array([0.01, 0.2, 0.5])
+    cases = eqsig.surface.calc_cum_abs_surface_energy(accsig, tshifts, nodal=True, up_red=1, down_red=1, stt=0.0,
+                                                      trim=False, start=False)
+    cases_w_s = eqsig.surface.calc_cum_abs_surface_energy(accsig, tshifts, nodal=True, up_red=1, down_red=1, stt=0.0,
+                                                      trim=False, start=True)
+    diff = np.sum(abs(cases[:, -1] - cases_w_s[:, -1]))
+    assert diff < 2.0e-5, diff
+
+
+def test_calc_cum_abs_surface_energy_start_from_base():
+    accsig = conftest.t_asig()
+    travel_times = np.array([0.01, 0.2, 0.5])
+    stt = 1.0
+    cases_w_s = eqsig.surface.calc_cum_abs_surface_energy(accsig, travel_times, nodal=True, up_red=1, down_red=1,
+                                                          stt=stt, trim=True, start=True)
+
+    cases = eqsig.surface.calc_cum_abs_surface_energy(accsig, travel_times, nodal=True, up_red=1, down_red=1, stt=stt,
+                                                      trim=True, start=False)
+    case_interp_0 = np.interp(accsig.time + (stt - travel_times[0]), accsig.time, cases_w_s[0])
+    diff0 = np.sum(abs(case_interp_0 - cases[0])) / cases[0][-1]
+    assert np.isclose(diff0, 0., atol=1.0e-7), diff0
+    case_interp_1 = np.interp(accsig.time + (stt - travel_times[1]), accsig.time, cases_w_s[1])
+    diff1 = np.sum(abs(case_interp_1 - cases[1])) / cases[1][-1]
+    assert np.isclose(diff1, 0., atol=1.0e-7), diff1
+    case_interp_2 = np.interp(accsig.time + (stt - travel_times[2]), accsig.time, cases_w_s[2])
+    diff2 = np.sum(abs(case_interp_2 - cases[2])) / cases[2][-1]
+    assert np.isclose(diff2, 0., atol=1.0e-7), diff2
+
+
+def test_calc_cum_abs_surface_energy_start_from_top():
+    accsig = conftest.t_asig()
+    travel_times = np.array([0.01, 0.2, 0.5])
+    stt = 0.0
+    cases_w_s = eqsig.surface.calc_cum_abs_surface_energy(accsig, travel_times, nodal=True, up_red=1, down_red=1,
+                                                          stt=stt, trim=True, start=True)
+
+    cases = eqsig.surface.calc_cum_abs_surface_energy(accsig, travel_times, nodal=True, up_red=1, down_red=1, stt=stt,
+                                                      trim=True, start=False)
+    case_interp_0 = np.interp(accsig.time + (stt - travel_times[0]), accsig.time, cases_w_s[0])
+    diff0 = np.sum(abs(case_interp_0 - cases[0])) / cases[0][-1]
+    assert np.isclose(diff0, 0., atol=5.0e-3), diff0
+    case_interp_1 = np.interp(accsig.time + (stt - travel_times[1]), accsig.time, cases_w_s[1])
+    diff1 = np.sum(abs(case_interp_1 - cases[1])) / cases[1][-1]
+    assert np.isclose(diff1, 0., atol=5.0e-3), diff1
+    case_interp_2 = np.interp(accsig.time + (stt - travel_times[2]), accsig.time, cases_w_s[2])
+    diff2 = np.sum(abs(case_interp_2 - cases[2])) / cases[2][-1]
+    assert np.isclose(diff2, 0., atol=8.0e-2), diff2
+
+
+
+def skip_plot():
+    import matplotlib.pyplot as plt
+    accsig = eqsig.load_asig(conftest.TEST_DATA_DIR + 'test_motion_dt0p01.txt')
+    accsig = eqsig.AccSignal(accsig.values[100:], accsig.dt)
+    travel_times = np.linspace(0, 2, 44)
+    travel_times = np.array([0.2, 0.5])
+    stt = 0.0
+    cases_w_s = eqsig.surface.calc_cum_abs_surface_energy(accsig, travel_times, nodal=True, up_red=1, down_red=1, stt=stt,
+                                                      trim=True, start=True)
+    # plt.plot(tshifts, cases[:, -1])
+    # plt.plot(tshifts, expected_cases[:, -1])
+    from bwplot import cbox
+    plt.plot(accsig.time, cases_w_s[0], c=cbox(0))
+    plt.plot(accsig.time, cases_w_s[1], c='r')
+    cases = eqsig.surface.calc_cum_abs_surface_energy(accsig, travel_times, nodal=True, up_red=1, down_red=1, stt=stt,
+                                                      trim=True, start=False)
+    plt.plot(accsig.time, cases[0], c=cbox(1), ls='--')
+    plt.plot(accsig.time, cases[1], c=cbox(1), ls='--')
+    plt.plot(accsig.time + (stt - travel_times[0]), cases[0], c='k', ls=':')
+    plt.plot(accsig.time + (stt - travel_times[1]), cases[1], c='k', ls=':')
+    case_interp_0 = np.interp(accsig.time + (stt - travel_times[0]), accsig.time, cases_w_s[0])
+    diff0 = np.sum(abs(case_interp_0 - cases[0])) / cases[0][-1]
+    assert np.isclose(diff0, 0., atol=5.0e-3), diff0
+    case_interp_1 = np.interp(accsig.time + (stt - travel_times[1]), accsig.time, cases_w_s[1])
+    diff1 = np.sum(abs(case_interp_1 - cases[1])) / cases[1][-1]
+    assert np.isclose(diff1, 0., atol=5.0e-3), diff1
+    case_interp_2 = np.interp(accsig.time, accsig.time + (stt - travel_times[2]), cases_w_s[2])
+    diff2 = np.sum(abs(case_interp_2 - cases[2])) / cases[2][-1]
+    assert np.isclose(diff2, 0., atol=5.0e-3), diff2
+
+    plt.show()
+
+def skip_plot2():
+    import matplotlib.pyplot as plt
+    asig = eqsig.load_asig(conftest.TEST_DATA_DIR + 'test_motion_dt0p01.txt')
+    travel_times = np.array([0.2, 0.5])
+    cases = eqsig.surface.calc_cum_abs_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, stt=1.0,
+                                                      trim=False, start=True)
+    # plt.plot(tshifts, cases[:, -1])
+    # plt.plot(tshifts, expected_cases[:, -1])
+    plt.plot(cases[0])
+    plt.plot(cases[1])
+
+    cases = eqsig.surface.calc_cum_abs_surface_energy(asig, travel_times, nodal=True, up_red=1, down_red=1, stt=1.0,
+                                                      trim=False, start=False)
+    plt.plot(cases[0])
+    plt.plot(cases[1])
+
+    # plt.plot(asig.time - 2 * travel_times[0], cases[0])
+    # plt.plot(asig.time - 2 * travel_times[1], cases[1])
+    # plt.plot(asig.time - 2 * travel_times[2], cases[2])
+    plt.show()
+
+
 if __name__ == '__main__':
-    test_calc_cum_abs_surface_energy()
+    # test_calc_cum_abs_surface_energy_start_check_same_energy()
+    test_calc_cum_abs_surface_energy_start_from_top()
+    # skip_plot()
+    # test_calc_cum_abs_surface_energy()