Merge pull request #274 from rsagroup/demo_boot

Demo updates

Author: HeikoSchuett

demos/demo_bootstrap.ipynb

@@ -144,17 +144,19 @@ def eval_dual_bootstrap(
     if use_correction and n_cv > 1:
         # we essentially project from the two points for 1 repetition and
         # for n_cv repetitions to infinitely many cv repetitions
-        evals_mean = np.mean(np.mean(evaluations[eval_ok], -1), -1)
-        evals_1 = np.mean(evaluations[eval_ok], -2)
-        noise_ceil_mean = np.mean(noise_ceil[:, eval_ok], -1)
-        noise_ceil_1 = noise_ceil[:, eval_ok]
-        var_mean = np.cov(
-            np.concatenate([evals_mean.T, noise_ceil_mean]))
-        var_1 = []
-        for i in range(n_cv):
-            var_1.append(np.cov(np.concatenate([
-                evals_1[:, :, i].T, noise_ceil_1[:, :, i]])))
-        var_1 = np.mean(np.array(var_1), axis=0)
+        evals_nonan = np.mean(np.mean(evaluations[eval_ok], -2), -2)
+        evals_1 = np.mean(evaluations[eval_ok], -3)
+        noise_ceil_nonan = np.mean(
+            noise_ceil[:, eval_ok], -2).transpose([1, 0, 2])
+        noise_ceil_1 = noise_ceil[:, eval_ok].transpose([1, 0, 2, 3])
+        matrix = np.concatenate([evals_nonan, noise_ceil_nonan], 1)
+        matrix -= np.mean(matrix, 0, keepdims=True)
+        var_mean = np.einsum('ijk,ilk->kjl', matrix, matrix) \
+            / (matrix.shape[0] - 1)
+        matrix_1 = np.concatenate([evals_1, noise_ceil_1], 1)
+        matrix_1 -= np.mean(matrix_1, 0, keepdims=True)
+        var_1 = np.einsum('ijmk,ilmk->kjl', matrix_1, matrix_1) \
+            / (matrix_1.shape[0] - 1) / matrix_1.shape[2]
         # this is the main formula for the correction:
         variances = (n_cv * var_mean - var_1) / (n_cv - 1)
     else:

demos/demo_dissimilarities.ipynb

@@ -102,14 +102,14 @@ def summary(self, test_type='t-test'):
         p_noise = self.test_noise(test_type=test_type)
         # header of the results table
         summary += 'Model' + (' ' * (name_length - 5))
-        summary += '|  Eval \u00B1 SEM   |'
+        summary += '|   Eval \u00B1 SEM   |'
         summary += ' p (against 0) |'
         summary += ' p (against NC) |\n'
-        summary += '-' * (name_length + 50)
+        summary += '-' * (name_length + 51)
         summary += '\n'
         for i, m in enumerate(self.models):
             summary += m.name + (' ' * (name_length - len(m.name)))
-            summary += f'| {means[i]:5.3f} \u00B1 {sems[i]:4.3f} |'
+            summary += f'| {means[i]: 5.3f} \u00B1 {sems[i]:4.3f} |'
             if p_zero[i] < 0.001:
                 summary += '      < 0.001  |'
             else:

rsatoolbox/inference/evaluate.py

@@ -296,7 +296,9 @@ def show_rdm_panel(
     rdmat = rdm.get_matrices()[0, :, :]
     if np.any(nanmask):
         rdmat[nanmask] = np.nan
-    image = ax.imshow(rdmat, cmap=cmap, vmin=vmin, vmax=vmax)
+    image = ax.imshow(
+        rdmat, cmap=cmap, vmin=vmin, vmax=vmax,
+        interpolation='none')
     ax.set_xlim(-0.5, rdm.n_cond - 0.5)
     ax.set_ylim(rdm.n_cond - 0.5, -0.5)
     ax.xaxis.set_ticks(gridlines)

rsatoolbox/inference/result.py

@@ -12,7 +12,6 @@ class TestDemos(unittest.TestCase):
     def test_example_dataset(self):
         import numpy as np
         import matplotlib.pyplot as plt
-        import rsatoolbox
         import rsatoolbox.data as rsd  # abbreviation to deal with dataset
 
         # import the measurements for the dataset
@@ -101,8 +100,6 @@ def test_example_dataset(self):
     def test_example_dissimilarities(self):
         # relevant imports
         import numpy as np
-        from scipy import io
-        import rsatoolbox
         import rsatoolbox.data as rsd  # abbreviation to deal with dataset
         import rsatoolbox.rdm as rsr
         # create a dataset object
@@ -122,11 +119,12 @@ def test_example_dissimilarities(self):
                            channel_descriptors=chn_des)
         # calculate an RDM
         RDM_euc = rsr.calc_rdm(data)
-        RDM_corr = rsr.calc_rdm(data,
-                                method='correlation', descriptor='conds')
+        _ = rsr.calc_rdm(
+            data,
+            method='correlation', descriptor='conds')
         # create an RDM object
         rdm_des = {'RDM': np.array(['RDM_1'])}
-        RDM_euc2 = rsr.RDMs(
+        _ = rsr.RDMs(
             RDM_euc.dissimilarities,
             dissimilarity_measure=RDM_euc.dissimilarity_measure,
             descriptors=RDM_euc.descriptors,
@@ -136,6 +134,49 @@ def test_example_dissimilarities(self):
         dist_matrix = RDM_euc.get_matrices()
         print(dist_matrix)
 
+    def test_dual_boot(self):
+        import numpy as np
+        from scipy import io
+        import rsatoolbox
+        import os
+        path = os.path.dirname(os.path.abspath(__file__))
+        matlab_data = io.matlab.loadmat(
+            os.path.join(path, '..', 'demos',
+                         'rdms_inferring', 'modelRDMs_A2020.mat'))
+        matlab_data = matlab_data['modelRDMs']
+        n_models = len(matlab_data[0])
+        model_names = [matlab_data[0][i][0][0] for i in range(n_models)]
+        measurement_model = [matlab_data[0][i][1][0] for i in range(n_models)]
+        rdms_array = np.array([matlab_data[0][i][3][0] for i in range(n_models)])
+        model_rdms = rsatoolbox.rdm.RDMs(
+            rdms_array,
+            rdm_descriptors={
+                'brain_computational_model': model_names,
+                'measurement_model': measurement_model},
+            dissimilarity_measure='Euclidean'
+        )
+        model_names = [matlab_data[0][i][0][0] for i in range(n_models)]
+        matlab_data = io.matlab.loadmat(
+            os.path.join(path, '..', 'demos',
+                         'rdms_inferring', 'noisyModelRDMs_A2020.mat'))
+        rdms_matlab = matlab_data['noisyModelRDMs']
+        rdms_matrix = rdms_matlab.squeeze().astype('float')
+        i_rep = 2  # np.random.randint(len(repr_names))
+        i_noise = 1  # np.random.randint(len(noise_std))
+        i_fwhm = 0  # np.random.randint(len(fwhms))
+        rdms_data = rsatoolbox.rdm.RDMs(rdms_matrix[:, i_rep, i_fwhm, i_noise, :].transpose())
+        models_flex = []
+        for i_model in np.unique(model_names):
+            models_flex.append(
+                rsatoolbox.model.ModelSelect(
+                    i_model,
+                    model_rdms.subset('brain_computational_model', i_model)
+                    )
+                )
+        results_3_full = rsatoolbox.inference.eval_dual_bootstrap(
+            models_flex, rdms_data, k_pattern=4, k_rdm=2, method='corr', N=5)
+        print(results_3_full)
+
     def test_exercise_all(self):
         import numpy as np
         from scipy import io
@@ -161,8 +202,9 @@ def test_exercise_all(self):
 
         conv1_rdms = model_rdms.subset('brain_computational_model', 'conv1')
         plt.figure(figsize=(10, 10))
-        rsatoolbox.vis.show_rdm(conv1_rdms,
-                           rdm_descriptor='measurement_model')
+        rsatoolbox.vis.show_rdm(
+            conv1_rdms,
+            rdm_descriptor='measurement_model')
 
         conv1_rdms = model_rdms.subset('brain_computational_model', 'conv1')
         print(conv1_rdms)
@@ -229,8 +271,10 @@ def test_exercise_all(self):
 
         models_flex = []
         for i_model in np.unique(model_names):
-            models_flex.append(rsatoolbox.model.ModelSelect(i_model,
-                model_rdms.subset('brain_computational_model', i_model)))
+            models_flex.append(
+                rsatoolbox.model.ModelSelect(
+                    i_model,
+                    model_rdms.subset('brain_computational_model', i_model)))
 
         print('created the following models:')
         for i in range(len(models_flex)):
@@ -332,9 +376,10 @@ def test_temporal_rsa(self):
             '%0.0f ms' % (np.round(x * 1000, 2))
             for x in rdms_data_binned.rdm_descriptors['time']]
 
-        rsatoolbox.vis.show_rdm(rdms_data_binned,
-                           pattern_descriptor='conds',
-                           rdm_descriptor='time_formatted')
+        rsatoolbox.vis.show_rdm(
+            rdms_data_binned,
+            pattern_descriptor='conds',
+            rdm_descriptor='time_formatted')
         from rsatoolbox.rdm import get_categorical_rdm
         rdms_model_in = get_categorical_rdm(['%d' % x for x in range(4)])
         rdms_model_lr = get_categorical_rdm(['l', 'r', 'l', 'r'])
@@ -347,8 +392,10 @@ def test_temporal_rsa(self):
         model_rdms.rdm_descriptors['model_names'] = model_names
         model_rdms.pattern_descriptors['cond_names'] = cond_names
         plt.figure(figsize=(10, 10))
-        rsatoolbox.vis.show_rdm(model_rdms, rdm_descriptor='model_names',
-                           pattern_descriptor='cond_names')
+        rsatoolbox.vis.show_rdm(
+            model_rdms,
+            rdm_descriptor='model_names',
+            pattern_descriptor='cond_names')
         from rsatoolbox.rdm import compare
         r = []
         for mod in model_rdms:
@@ -375,16 +422,18 @@ def test_demo_rdm_scatterplot(self):
         matlab_data = loadmat(
             os.path.join(path, '..', 'demos', '92imageData/92_brainRDMs.mat'))['RDMs']
         n_rdms = len(matlab_data[0])
-        rdms_ = RDMs(np.array([matlab_data[0][i][0][0] for i in range(n_rdms)]),
-                  pattern_descriptors=condition_descriptors,
-                  rdm_descriptors={'name': np.array([f"RDM{i}"
-                                                     for i in range(4)])}
-                  )
+        rdms_ = RDMs(
+            np.array([matlab_data[0][i][0][0] for i in range(n_rdms)]),
+            pattern_descriptors=condition_descriptors,
+            rdm_descriptors={
+                'name': np.array([f"RDM{i}" for i in range(4)])}
+            )
 
         rdms_a = concat([rdms_[0], rdms_[1]])
         rdms_b = concat([rdms_[2], rdms_[3]])
 
-        rdm_comparison_scatterplot((rdms_a, rdms_b),
+        rdm_comparison_scatterplot(
+            (rdms_a, rdms_b),
             show_marginal_distributions=True,
             show_identity_line=True,
             show_legend=False,
@@ -397,5 +446,6 @@ def test_demo_rdm_scatterplot(self):
             }
             )
 
+
 if __name__ == '__main__':
     unittest.main()

rsatoolbox/vis/rdm_plot.py

@@ -338,10 +338,10 @@ def test_result_print(self):
         self.assertEqual(
             string,
             'Results for running fixed evaluation for corr on 2 models:\n\n'
-            + 'Model |  Eval ± SEM   | p (against 0) | p (against NC) |\n'
-            + '--------------------------------------------------------\n'
-            + 't1    | 0.873 ± 0.100 |        0.006  |         0.496  |\n'
-            + 't2    | 0.120 ± 0.141 |        0.243  |         0.037  |\n\n'
+            + 'Model |   Eval ± SEM   | p (against 0) | p (against NC) |\n'
+            + '---------------------------------------------------------\n'
+            + 't1    |  0.873 ± 0.100 |        0.006  |         0.496  |\n'
+            + 't2    |  0.120 ± 0.141 |        0.243  |         0.037  |\n\n'
             + 'p-values are based on uncorrected t-tests')