Update test_unitary_event_analysis.py with validation test (#152)

* Update test_unitary_event_analysis.py

test added to reproduce the results of Riehle et al. 1997 science using unitary event method (implemented in Elephant) and check the consistency with the original publication (see Rostami et al. 1997 [Re] science for a detailed explanation)

* Update test_unitary_event_analysis.py

used mock library to exclude importing matplotlib

* Update test_unitary_event_analysis.py

Added neo loading routine to the test file.

* Update test_unitary_event_analysis.py

changed indentation and `numpy` to `np`

* Update test_unitary_event_analysis.py

removed mock

* Update test_unitary_event_analysis.py

Removed argument `encoding`

* Update test_unitary_event_analysis.py

added `try except` block around `np.load` to support loading  for python 2 and 3

* Update test_unitary_event_analysis.py

except changed to `UnicodeError`

* Update test_unitary_event_analysis.py

added decorator to skip certain configuration in the unittest

Author: Vahidrostami

elephant/test/test_unitary_event_analysis.py

@@ -11,6 +11,16 @@
 import types
 import elephant.unitary_event_analysis as ue
 import neo
+import sys
+import os
+
+from distutils.version import StrictVersion
+
+
+def _check_for_incompatibilty():
+    smaller_version = StrictVersion(np.__version__) < '1.10.0'
+    return sys.version_info >= (3, 0) and smaller_version
+
 
 class UETestCase(unittest.TestCase):
 
@@ -338,6 +348,149 @@ def test_jointJ_window_analysis(self):
             UE_dic['indices']['trial26'],expected_indecis_tril26))
         self.assertTrue(np.allclose(
             UE_dic['indices']['trial4'],expected_indecis_tril4))
+        
+    @staticmethod    
+    def load_gdf2Neo(fname, trigger, t_pre, t_post):
+        """
+        load and convert the gdf file to Neo format by
+        cutting and aligning around a given trigger
+        # codes for trigger events (extracted from a
+        # documentation of an old file after
+        # contacting Dr. Alexa Rihle)
+        # 700 : ST (correct) 701, 702, 703, 704*
+        # 500 : ST (error =5) 501, 502, 503, 504*
+        # 1000: ST (if no selec) 1001,1002,1003,1004*
+        # 11  : PS 111, 112, 113, 114
+        # 12  : RS 121, 122, 123, 124
+        # 13  : RT 131, 132, 133, 134
+        # 14  : MT 141, 142, 143, 144
+        # 15  : ES 151, 152, 153, 154
+        # 16  : ES 161, 162, 163, 164
+        # 17  : ES 171, 172, 173, 174
+        # 19  : RW 191, 192, 193, 194
+        # 20  : ET 201, 202, 203, 204
+        """
+        data = np.loadtxt(fname)
+    
+        if trigger == 'PS_4':
+            trigger_code = 114
+        if trigger == 'RS_4':
+            trigger_code = 124
+        if trigger == 'RS':
+            trigger_code = 12
+        if trigger == 'ES':
+            trigger_code = 15
+        # specify units
+        units_id = np.unique(data[:, 0][data[:, 0] < 7])
+        # indecies of the trigger
+        sel_tr_idx = np.where(data[:, 0] == trigger_code)[0]
+        # cutting the data by aligning on the trigger
+        data_tr = []
+        for id_tmp in units_id:
+            data_sel_units = []
+            for i_cnt, i in enumerate(sel_tr_idx):
+                start_tmp = data[i][1] - t_pre.magnitude
+                stop_tmp = data[i][1] + t_post.magnitude
+                sel_data_tmp = np.array(
+                    data[np.where((data[:, 1] <= stop_tmp) &
+                                     (data[:, 1] >= start_tmp))])
+                sp_units_tmp = sel_data_tmp[:, 1][
+                    np.where(sel_data_tmp[:, 0] == id_tmp)[0]]
+                if len(sp_units_tmp) > 0:
+                    aligned_time = sp_units_tmp - start_tmp
+                    data_sel_units.append(neo.SpikeTrain(
+                        aligned_time * pq.ms, t_start=0 * pq.ms,
+                        t_stop=t_pre + t_post))
+                else:
+                    data_sel_units.append(neo.SpikeTrain(
+                        [] * pq.ms, t_start=0 * pq.ms,
+                        t_stop=t_pre + t_post))
+            data_tr.append(data_sel_units)
+        data_tr.reverse()
+        spiketrain = np.vstack([i for i in data_tr]).T
+        return spiketrain
+
+    # test if the result of newly implemented Unitary Events in
+    # Elephant is consistent with the result of
+    # Riehle et al 1997 Science
+    # (see Rostami et al (2016) [Re] Science, 3(1):1-17)
+    @unittest.skipIf(_check_for_incompatibilty(),
+                     'Incompatible package versions')
+    def test_Riehle_et_al_97_UE(self):      
+        from neo.rawio.tests.tools import (download_test_file,
+                                           create_local_temp_dir,
+                                           make_all_directories)
+        from neo.test.iotest.tools import (cleanup_test_file)
+        url = [
+            "https://raw.githubusercontent.com/ReScience-Archives/" +
+            "Rostami-Ito-Denker-Gruen-2017/master/data",
+            "https://raw.githubusercontent.com/ReScience-Archives/" +
+            "Rostami-Ito-Denker-Gruen-2017/master/data"]
+        shortname = "unitary_event_analysis_test_data"
+        local_test_dir = create_local_temp_dir(
+            shortname, os.environ.get("ELEPHANT_TEST_FILE_DIR"))
+        files_to_download = ["extracted_data.npy", "winny131_23.gdf"]
+        make_all_directories(files_to_download,
+                             local_test_dir)
+        for f_cnt, f in enumerate(files_to_download):
+            download_test_file(f, local_test_dir, url[f_cnt])
+
+        # load spike data of figure 2 of Riehle et al 1997
+        sys.path.append(local_test_dir)
+        file_name = '/winny131_23.gdf'
+        trigger = 'RS_4'
+        t_pre = 1799 * pq.ms
+        t_post = 300 * pq.ms
+        spiketrain = self.load_gdf2Neo(local_test_dir + file_name,
+                                       trigger, t_pre, t_post)
+
+        # calculating UE ...
+        winsize = 100 * pq.ms
+        binsize = 5 * pq.ms
+        winstep = 5 * pq.ms
+        pattern_hash = [3]
+        method = 'analytic_TrialAverage'
+        t_start = spiketrain[0][0].t_start
+        t_stop = spiketrain[0][0].t_stop
+        t_winpos = ue._winpos(t_start, t_stop, winsize, winstep)
+        significance_level = 0.05
+
+        UE = ue.jointJ_window_analysis(
+            spiketrain, binsize, winsize, winstep,
+            pattern_hash, method=method)
+        # load extracted data from figure 2 of Riehle et al 1997
+        try:
+            extracted_data = np.load(
+                local_test_dir + '/extracted_data.npy').item()
+        except UnicodeError:
+            extracted_data = np.load(
+                local_test_dir + '/extracted_data.npy', encoding='latin1').item()
+        Js_sig = ue.jointJ(significance_level)
+        sig_idx_win = np.where(UE['Js'] >= Js_sig)[0]
+        diff_UE_rep = []
+        y_cnt = 0
+        for tr in range(len(spiketrain)):
+            x_idx = np.sort(
+                np.unique(UE['indices']['trial' + str(tr)],
+                          return_index=True)[1])
+            x = UE['indices']['trial' + str(tr)][x_idx]
+            if len(x) > 0:
+                # choose only the significant coincidences
+                xx = []
+                for j in sig_idx_win:
+                    xx = np.append(xx, x[np.where(
+                        (x * binsize >= t_winpos[j]) &
+                        (x * binsize < t_winpos[j] + winsize))])
+                x_tmp = np.unique(xx) * binsize.magnitude
+                if len(x_tmp) > 0:
+                    ue_trial = np.sort(extracted_data['ue'][y_cnt])
+                    diff_UE_rep = np.append(
+                        diff_UE_rep, x_tmp - ue_trial)
+                    y_cnt += +1
+        np.testing.assert_array_less(np.abs(diff_UE_rep), 0.3)
+        cleanup_test_file('dir', local_test_dir)
+
+        
 def suite():
     suite = unittest.makeSuite(UETestCase, 'test')
     return suite