[revert] changed slope calculation in AICPicker back to older state (probably causing problems changing results in test_autopickstation.py)

pylot/core/pick/charfuns.py

@@ -59,7 +59,7 @@ class CharacteristicFunction(object):
         self.setOrder(order)
         self.setFnoise(fnoise)
         self.setARdetStep(t2)
-        self.calcCF(self.getDataArray())
+        self.calcCF()
         self.arpara = np.array([])
         self.xpred = np.array([])
 
@@ -211,17 +211,15 @@ class CharacteristicFunction(object):
             data = self.orig_data.copy()
             return data
 
-    def calcCF(self, data=None):
-        self.cf = data
+    def calcCF(self):
+        pass
 
 
 class AICcf(CharacteristicFunction):
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate the Akaike Information Criterion (AIC) after Maeda (1985).
-        :param data: data, time series (whether seismogram or CF)
-        :type data: tuple
         :return: AIC function
         :rtype:
         """
@@ -259,13 +257,11 @@ class HOScf(CharacteristicFunction):
         """
         super(HOScf, self).__init__(data, cut, pickparams["tlta"], pickparams["hosorder"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate skewness (statistics of order 3) or kurtosis
         (statistics of order 4), using one long moving window, as published
         in Kueperkoch et al. (2010), or order 2, i.e. STA/LTA.
-        :param data: data, time series (whether seismogram or CF)
-        :type data: tuple
         :return: HOS cf
         :rtype:
         """
@@ -280,47 +276,28 @@ class HOScf(CharacteristicFunction):
         elif self.getOrder() == 4:  # this is kurtosis
             y = np.power(xnp, 4)
             y1 = np.power(xnp, 2)
-        elif self.getOrder() == 2:  # this is variance, used for STA/LTA processing
-            y = np.power(xnp, 2)
-            y1 = np.power(xnp, 2)
 
         # Initialisation
         # t2: long term moving window
         ilta = int(round(self.getTime2() / self.getIncrement()))
-        ista = int(round((self.getTime2() / 10) / self.getIncrement()))  # TODO: still hard coded!!
         lta = y[0]
         lta1 = y1[0]
-        sta = y[0]
         # moving windows
         LTA = np.zeros(len(xnp))
-        STA = np.zeros(len(xnp))
         for j in range(0, len(xnp)):
             if j < 4:
                 LTA[j] = 0
-                STA[j] = 0
-            elif j <= ista and self.getOrder() == 2:
-                lta = (y[j] + lta * (j - 1)) / j
-                if self.getOrder() == 2:
-                    sta = (y[j] + sta * (j - 1)) / j
-            # elif j < 4:
             elif j <= ilta:
                 lta = (y[j] + lta * (j - 1)) / j
                 lta1 = (y1[j] + lta1 * (j - 1)) / j
-                if self.getOrder() == 2:
-                    sta = (y[j] - y[j - ista]) / ista + sta
             else:
                 lta = (y[j] - y[j - ilta]) / ilta + lta
                 lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1
-                if self.getOrder() == 2:
-                    sta = (y[j] - y[j - ista]) / ista + sta
             # define LTA
             if self.getOrder() == 3:
                 LTA[j] = lta / np.power(lta1, 1.5)
             elif self.getOrder() == 4:
                 LTA[j] = lta / np.power(lta1, 2)
-            else:
-                LTA[j] = lta
-                STA[j] = sta
 
         # remove NaN's with first not-NaN-value,
         # so autopicker doesnt pick discontinuity at start of the trace
@@ -329,10 +306,7 @@ class HOScf(CharacteristicFunction):
             first = ind[0]
             LTA[:first] = LTA[first]
 
-        if self.getOrder() > 2:
         self.cf = LTA
-        else:  # order 2 means STA/LTA!
-            self.cf = STA / LTA
         self.xcf = x
 
 
@@ -342,12 +316,10 @@ class ARZcf(CharacteristicFunction):
         super(ARZcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Parorder"],
                                     fnoise=pickparams["addnoise"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         function used to calculate the AR prediction error from a single vertical trace. Can be used to pick
         P onsets.
-        :param data:
-        :type data: ~obspy.core.stream.Stream
         :return: ARZ cf
         :rtype:
         """
@@ -478,14 +450,12 @@ class ARHcf(CharacteristicFunction):
         super(ARHcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
                                     fnoise=pickparams["addnoise"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate a characteristic function using autoregressive modelling of the waveform of
         both horizontal traces.
         The waveform is predicted in a moving time window using the calculated AR parameters. The difference
         between the predicted and the actual waveform servers as a characteristic function.
-        :param data: wavefor stream
-        :type data: ~obspy.core.stream.Stream
         :return: ARH cf
         :rtype:
         """
@@ -634,14 +604,12 @@ class AR3Ccf(CharacteristicFunction):
         super(AR3Ccf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
                                      fnoise=pickparams["addnoise"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate a characteristic function using autoregressive modelling of the waveform of
         all three traces.
         The waveform is predicted in a moving time window using the calculated AR parameters. The difference
         between the predicted and the actual waveform servers as a characteristic function
-        :param data: stream holding all three traces
-        :type data: ~obspy.core.stream.Stream
         :return: AR3C cf
         :rtype:
         """

pylot/core/pick/picker.py

@@ -173,7 +173,7 @@ class AICPicker(AutoPicker):
         nn = np.isnan(self.cf)
         if len(nn) > 1:
             self.cf[nn] = 0
-        # taper AIC-CF to get rid off side maxima
+        # taper AIC-CF to get rid of side maxima
         tap = np.hanning(len(self.cf))
         aic = tap * self.cf + max(abs(self.cf))
         # smooth AIC-CF
@@ -316,16 +316,7 @@ class AICPicker(AutoPicker):
                             plt.close(fig)
                     return
                 iislope = islope[0][0:imax + 1]
-            # MP MP change slope calculation
-            # get all maxima of aicsmooth
-            iaicmaxima = argrelmax(aicsmooth)[0]
-            # get first index of maximum after pickindex (indices saved in iaicmaxima)
-            aicmax = iaicmaxima[np.where(iaicmaxima > pickindex)[0]]
-            if len(aicmax) > 0:
-                iaicmax = aicmax[0]
-            else:
-                iaicmax = -1
-            dataslope = aicsmooth[pickindex: iaicmax]
+            dataslope = self.Data[0].data[iislope]
             # calculate slope as polynomal fit of order 1
             xslope = np.arange(0, len(dataslope), 1)
             try:
@@ -336,7 +327,7 @@ class AICPicker(AutoPicker):
                 else:
                     self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
                     # normalize slope to maximum of cf to make it unit independent
-                    self.slope /= aicsmooth[iaicmax]
+                    self.slope /= self.Data[0].data[icfmax]
             except Exception as e:
                 print("AICPicker: Problems with data fitting! {}".format(e))
 
@@ -376,7 +367,7 @@ class AICPicker(AutoPicker):
                             label='Signal Window')
                 ax2.axvspan(self.Tcf[iislope[0]], self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0,
                             label='Slope Window')
-                ax2.plot(self.Tcf[pickindex: iaicmax], datafit, 'g', linewidth=2,
+                ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2,
                          label='Slope')  # MP MP changed temporarily!
 
                 if self.slope is not None:

tests/test_autopickstation/test_autopickstation.py

@@ -272,7 +272,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=wfstream, pickparam=self.pickparam_taupy_disabled,
                                               metadata=(None, None))
-        compare_dicts(result, expected)
+        compare_dicts(result=result['P'], expected=expected['P'])
+        compare_dicts(result=result['S'], expected=expected['S'])
         self.assertEqual('GRA1', station)
 
     def test_autopickstation_a106_taupy_enabled(self):
@@ -290,7 +291,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.a106, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        compare_dicts(result=result, expected=expected)
+        compare_dicts(result=result['P'], expected=expected['P'])
+        compare_dicts(result=result['S'], expected=expected['S'])
 
 
     def test_autopickstation_station_missing_in_metadata(self):
@@ -312,7 +314,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.a005a, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        compare_dicts(result, expected)
+        compare_dicts(result=result['P'], expected=expected['P'])
+        compare_dicts(result=result['S'], expected=expected['S'])
 
 
 def run_dict_comparison(result, expected):
@@ -332,8 +335,8 @@ def compare_dicts(result, expected):
         run_dict_comparison(result, expected)
     except AssertionError:
         raise AssertionError(f'Dictionaries not equal.'
-                             f'\n<<Expected>>: \n\n{pretty_print_dict(expected)}'
-                             f'\n<<Result>>: \n{pretty_print_dict(result)}')
+                             f'\n\n<<Expected>>\n{pretty_print_dict(expected)}'
+                             f'\n\n<<Result>>\n{pretty_print_dict(result)}')
 
 
 def pretty_print_dict(dct):