Subspace Freq domain statistic

appveyor.yml

@@ -36,7 +36,7 @@ install:
   - "conda create -q --yes -n test python=%PYTHON_VERSION%"
   - "activate test"
   # Install default dependencies
-  - "conda install -q --yes pip numpy scipy opencv3 matplotlib obspy mock flake8 pyflakes=0.9.0 cython h5py bottleneck"
+  - "conda install -q --yes pip numpy scipy opencv3 matplotlib==1.5.3 obspy==1.0.3 mock flake8 pyflakes=0.9.0 cython h5py bottleneck"
   # additional dependecies
   #  - "choco install opencv"
   #  - "powershell copy-item C:\\OpenCV\\opencv\\build\\python\\2.7\\x64\\cv2.pyd C:\\conda\\envs\\test\\lib\\site-packages\\."

eqcorrscan/core/match_filter.py

@@ -1349,8 +1349,7 @@ def __init__(self, name=None, st=None, lowcut=None, highcut=None,
         self.process_length = process_length
         self.prepick = prepick
         if event is not None:
-            if len(event.comments) > 0 and \
-               "eqcorrscan_template_" + temp_name not in \
+            if "eqcorrscan_template_" + temp_name not in \
                [c.text for c in event.comments]:
                 event.comments.append(Comment(
                         text="eqcorrscan_template_" + temp_name,
@@ -2888,7 +2887,10 @@ def _group_detect(templates, stream, threshold, threshold_type, trig_int,
     :return:
         :class:`eqcorrscan.core.match_filter.Party` of families of detections.
     """
-    ncores = cpu_count()
+    if parallel_process:
+        ncores = cpu_count()
+    else:
+        ncores = 1
     st = [Stream()]
     master = templates[0]
     # Check that they are all processed the same.
@@ -3306,8 +3308,9 @@ def extract_from_stream(stream, detections, pad=5.0, length=30.0):
                 print('No data in stream for pick:')
                 print(pick)
                 continue
-            cut_stream += tr.copy().trim(starttime=pick.time - pad,
-                                         endtime=pick.time - pad + length)
+            cut_stream += tr.slice(
+                starttime=pick.time - pad,
+                endtime=pick.time - pad + length).copy()
         streams.append(cut_stream)
     return streams
 

eqcorrscan/core/subspace.py

@@ -17,7 +17,6 @@
 from __future__ import unicode_literals
 
 import numpy as np
-import scipy
 import warnings
 import time
 import h5py
@@ -473,7 +472,6 @@ def _detect(detector, st, threshold, trig_int, moveout=0, min_trig=0,
     :return: list of detections
     :rtype: list of eqcorrscan.core.match_filter.Detection
     """
-    from eqcorrscan.core import subspace_statistic
     detections = []
     # First process the stream
     if process:
@@ -493,27 +491,25 @@ def _detect(detector, st, threshold, trig_int, moveout=0, min_trig=0,
         stream = [st]
         stachans = detector.stachans
     outtic = time.clock()
-    if debug > 0:
-        print('Computing detection statistics')
     # If multiplexed, how many samples do we increment by?
     if detector.multiplex:
-        inc = np.uint32(len(detector.stachans))
+        Nc = len(detector.stachans)
     else:
-        inc = np.uint32(1)
-    # Stats must be same size for multiplexed or non multiplexed!
+        Nc = 1
+    # Here do all ffts
+    fft_vars = do_ffts(detector, stream, Nc)
+    if debug > 0:
+        print('Computing detection statistics')
     if debug > 0:
         print('Preallocating stats matrix')
     stats = np.zeros((len(stream[0]),
-                      (len(stream[0][0]) // inc) -
-                      (len(detector.data[0].T[0]) // inc) + 1),
-                     dtype=np.float32)
-    # Hard typing in Cython loop requires float32 type.
-    for det_channel, in_channel, i in zip(detector.data, stream[0],
-                                          np.arange(len(stream[0]))):
-        stats[i] = subspace_statistic.\
-            det_statistic(detector=det_channel.astype(np.float32),
-                          data=in_channel.data.astype(np.float32),
-                          inc=inc)
+                      (len(stream[0][0]) // Nc) - (fft_vars[4] // Nc) + 1))
+    for det_freq, data_freq_sq, data_freq, i in zip(fft_vars[0], fft_vars[1],
+                                                    fft_vars[2],
+                                                    np.arange(len(stream[0]))):
+        # Calculate det_statistic in frequency domain
+        stats[i] = det_stat_freq(det_freq, data_freq_sq, data_freq,
+                                 fft_vars[3], Nc, fft_vars[4], fft_vars[5])
         if debug >= 1:
             print('Stats matrix is shape %s' % str(stats[i].shape))
         if debug >= 3:
@@ -525,8 +521,7 @@ def _detect(detector, st, threshold, trig_int, moveout=0, min_trig=0,
             ax.plot([min(t), max(t)], [threshold, threshold], color='r', lw=1,
                     label='Threshold')
             ax.legend()
-            plt.title('%s.%s' % (in_channel.stats.station,
-                                 in_channel.stats.channel))
+            plt.title('%s' % str(stream[0][i].stats.station))
             plt.show()
     trig_int_samples = detector.sampling_rate * trig_int
     if debug > 0:
@@ -586,6 +581,39 @@ def _detect(detector, st, threshold, trig_int, moveout=0, min_trig=0,
     return detections
 
 
+def do_ffts(detector, stream, Nc):
+    min_fftlen = int(stream[0][0].data.shape[0] +
+                     detector.data[0].shape[0] - Nc)
+    fftlen = 1 << min_fftlen.bit_length()
+    mplen = stream[0][0].data.shape[0]
+    ulen = detector.data[0].shape[0]
+    num_st_fd = [np.fft.fft(tr.data, n=fftlen)
+                 for tr in stream[0]]
+    denom_st_fd = [np.fft.fft(np.square(tr.data), n=fftlen)
+                   for tr in stream[0]]
+    # Frequency domain of boxcar
+    w = np.fft.fft(np.ones(detector.data[0].shape[0]),
+                   n=fftlen)
+    # This should go into the detector object as in Detex
+    detector_fd = []
+    for dat_mat in detector.data:
+        detector_fd.append(np.array([np.fft.fft(col[::-1], n=fftlen)
+                                     for col in dat_mat.T]))
+    return detector_fd, denom_st_fd, num_st_fd, w, ulen, mplen
+
+
+def det_stat_freq(det_freq, data_freq_sq, data_freq, w, Nc, ulen, mplen):
+    num_cor = np.multiply(det_freq, data_freq)  # Numerator convolution
+    den_cor = np.multiply(w, data_freq_sq)  # Denominator convolution
+    # Do inverse fft
+    # First and last Nt - 1 samples are invalid; clip them off
+    num_ifft = np.real(np.fft.ifft(num_cor))[:, ulen-1:mplen:Nc]
+    denominator = np.real(np.fft.ifft(den_cor))[ulen-1:mplen:Nc]
+    # Ratio of projected to envelope energy = det_stat across all channels
+    result = np.sum(np.square(num_ifft), axis=0) / denominator
+    return result
+
+
 def _subspace_process(streams, lowcut, highcut, filt_order, sampling_rate,
                       multiplex, align, shift_len, reject, no_missed=True,
                       stachans=None, parallel=False, plot=False):

eqcorrscan/tests/clustering_test.py

@@ -175,10 +175,10 @@ def test_svd(self):
                 tr.detrend('simple')
                 tr.filter('bandpass', freqmin=5.0, freqmax=15.0)
                 tr.trim(tr.stats.starttime + 40, tr.stats.endtime - 45)
-        SVectors, SValues, Uvectors, stachans = svd(stream_list=stream_list)
+        UVectors, SValues, SVectors, stachans = svd(stream_list=stream_list)
         self.assertEqual(len(SVectors), len(stachans))
         self.assertEqual(len(SValues), len(stachans))
-        self.assertEqual(len(Uvectors), len(stachans))
+        self.assertEqual(len(UVectors), len(stachans))
         for SVec in SVectors:
             self.assertEqual(len(SVec), len(stream_list))
         with warnings.catch_warnings(record=True) as w:
@@ -227,11 +227,11 @@ def test_svd_to_stream(self):
                 tr.resample(sampling_rate=samp_rate)
                 tr.trim(tr.stats.starttime + 40, tr.stats.endtime - 45)
         SVectors, SValues, Uvectors, stachans = svd(stream_list=stream_list)
-        svstreams = svd_to_stream(svectors=SVectors, stachans=stachans, k=4,
+        svstreams = svd_to_stream(uvectors=SVectors, stachans=stachans, k=4,
                                   sampling_rate=samp_rate)
         self.assertEqual(len(svstreams), 4)
         with warnings.catch_warnings(record=True) as w:
-            SVD_2_stream(SVectors=SVectors, stachans=stachans, k=4,
+            SVD_2_stream(uvectors=SVectors, stachans=stachans, k=4,
                          sampling_rate=samp_rate)
             self.assertEqual(len(w), 1)
             self.assertTrue('Depreciated' in str(w[0].message))

eqcorrscan/tests/lag_calc_test.py

@@ -17,6 +17,8 @@
 from eqcorrscan.core.match_filter import normxcorr2, Detection
 from eqcorrscan.utils.sfile_util import read_event
 
+warnings.simplefilter("always")
+
 
 class TestMethods(unittest.TestCase):
     @classmethod
@@ -233,6 +235,8 @@ def test_bad_interp(self):
         with warnings.catch_warnings(record=True) as w:
             warnings.simplefilter("always")
             _xcorr_interp(ccc, 0.1)
+        for _w in w:
+            print(_w.message)
         self.assertEqual(len(w), 2)
         self.assertTrue('Less than 5 samples' in str(w[0].message))
         self.assertTrue('Residual in quadratic fit' in str(w[1].message))

eqcorrscan/tests/subspace_test.py

@@ -14,7 +14,7 @@
 
 from obspy import Stream, read
 
-from eqcorrscan.core import subspace, subspace_statistic
+from eqcorrscan.core import subspace
 
 
 class SimpleSubspaceMethods(unittest.TestCase):
@@ -87,10 +87,13 @@ def test_stat(self):
         detector.partition(2)
         stream = read(os.path.join(os.path.abspath(os.path.dirname(__file__)),
                                    'test_data', 'subspace', 'test_trace.ms'))
-        tr_data = stream[0].data[:,None].astype(np.float32)
-        stat = subspace_statistic.det_statistic(
-            detector.data[0].astype(np.float32), tr_data, np.uint32(1))
-        self.assertEqual((stat.max().round(6) - 0.306929).round(6), 0)
+        st = [stream]
+        fft_vars = subspace.do_ffts(detector, st, len(detector.stachans))
+        stat = subspace.det_stat_freq(fft_vars[0][0], fft_vars[1][0],
+                                      fft_vars[2][0], fft_vars[3],
+                                      len(detector.stachans), fft_vars[4],
+                                      fft_vars[5])
+        self.assertEqual((stat.max().round(6) - 0.229755).round(6), 0)
 
 
 class SubspaceTestingMethods(unittest.TestCase):
@@ -327,7 +330,7 @@ def test_not_multiplexed(self):
         st = self.st
         detections = detector.detect(st=st, threshold=0.5, trig_int=4,
                                      debug=1, moveout=2, min_trig=5)
-        self.assertEqual(len(detections), 17)
+        self.assertEqual(len(detections), 16)
 
     def test_multi_detectors(self):
         """Test the efficient looping in subspace."""
@@ -349,14 +352,14 @@ def test_multi_detectors(self):
                                               trig_int=10, moveout=5,
                                               min_trig=5,
                                               parallel=False, num_cores=2)
-        self.assertEqual(len(detections), 14)
+        self.assertEqual(len(detections), 6)
         detections = subspace.subspace_detect(detectors=[detector1, detector2],
                                               stream=self.st.copy(),
                                               threshold=0.7,
                                               trig_int=10, moveout=5,
                                               min_trig=5,
                                               parallel=True, num_cores=2)
-        self.assertEqual(len(detections), 14)
+        self.assertEqual(len(detections), 6)
 
     def partition_fail(self):
         templates = copy.deepcopy(self.templates)

eqcorrscan/tests/synth_seis_test.py

@@ -52,7 +52,7 @@ def test_phaseout(self):
 
 class TestSVDSim(unittest.TestCase):
     def test_svd_sim(self):
-        V, s, U, stachans = SVD_sim(sp=15, lowcut=2, highcut=8, samp_rate=20)
+        U, s, V, stachans = SVD_sim(sp=15, lowcut=2, highcut=8, samp_rate=20)
         self.assertEqual(V[0].shape[0], s[0].shape[0])
         self.assertEqual(V[0].shape[0], U[0].shape[1])
         self.assertEqual(len(stachans), 1)

eqcorrscan/tests/tutorials_test.py

@@ -95,7 +95,7 @@ def test_lag_calc(self):
     def test_subspace(self):
         """Test the subspace tutorial."""
         detections = subspace.run_tutorial(plot=False)
-        self.assertEqual(len(detections), 2)
+        self.assertEqual(len(detections), 11)
 
 if __name__ == '__main__':
     """

eqcorrscan/tutorials/subspace.py

@@ -85,7 +85,7 @@ def run_tutorial(plot=False, multiplex=True, return_streams=False):
     # We set a very low threshold because the detector is not that great, we
     # haven't aligned it particularly well - however, at this threshold we make
     # two real detections.
-    detections, det_streams = detector.detect(st=st, threshold=0.005,
+    detections, det_streams = detector.detect(st=st, threshold=0.3,
                                               trig_int=2,
                                               extract_detections=True)
     if return_streams: