Refact: checked the static network toolbox paper example worked.

examples/toolbox_paper/ctf_rest/static_networks.py

@@ -33,8 +33,9 @@ def mat2vec(mat):
 
 # Load source data
 data = Data(
-    "/well/woolrich/projects/toolbox_paper/ctf_rest/training_data/networks",
+    "training_data/networks",
     sampling_frequency=250,
+    n_jobs=4,
 )
 
 # Calculate static power spectra
@@ -44,6 +45,7 @@ def mat2vec(mat):
     sampling_frequency=data.sampling_frequency,
     standardize=True,
     calc_coh=True,
+    n_jobs=4,
 )
 
 # Calculate power maps for different frequency bands
@@ -61,7 +63,11 @@ def mat2vec(mat):
     mean_pow_map,
     mask_file="MNI152_T1_8mm_brain.nii.gz",
     parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz",
-    plot_kwargs={"views": ["lateral"], "vmax": mean_pow_map.max()},
+    plot_kwargs={
+        "views": ["lateral"],
+        "symmetric_cbar": True,
+        "vmax": mean_pow_map.max(),  # puts all plots on the same scale
+    },
     filename=f"{output_dir}/pow_.png",
 )
 
@@ -82,14 +88,18 @@ def mat2vec(mat):
     mean_coh_net,
     parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz",
     threshold=0.95,
-    plot_kwargs={"edge_cmap": "Reds", "edge_vmin": 0, "edge_vmax": mean_coh_net.max()},
+    plot_kwargs={
+        "edge_cmap": "Reds",
+        "edge_vmin": 0,
+        "edge_vmax": mean_coh_net.max(),
+    },
     filename=f"{output_dir}/coh_.png",
 )
 
 # Calculate AEC networks for different frequency bands
 aec = []
 for l, h in [(1, 4), (4, 7), (7, 13), (13, 30)]:
-    data.filter(low_freq=l, high_freq=h)
+    data.filter(low_freq=l, high_freq=h, use_raw=True)
     data.amplitude_envelope()
     data.standardize()
     aec.append(static.functional_connectivity(data.time_series()))
@@ -104,7 +114,11 @@ def mat2vec(mat):
     mean_aec,
     parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz",
     threshold=0.95,
-    plot_kwargs={"edge_cmap": "Reds", "edge_vmin": 0, "edge_vmax": mean_aec.max()},
+    plot_kwargs={
+        "edge_cmap": "Reds",
+        "edge_vmin": 0,
+        "edge_vmax": mean_aec.max(),
+    },
     filename=f"{output_dir}/aec_.png",
 )
 
@@ -120,27 +134,30 @@ def mat2vec(mat):
 pow_diff, pvalues = statistics.group_diff_max_stat_perm(
     pow_map, assignments, n_perm=1000, n_jobs=4
 )  # pow_diff is group 1 (old) minus group 2 (young)
-pvalues[pvalues > 0.05] = 0.2
 
 # Plot significant power map differences
 power.save(
     pow_diff,
     mask_file="MNI152_T1_8mm_brain.nii.gz",
     parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz",
-    plot_kwargs={"views": ["lateral"]},
+    plot_kwargs={
+        "views": ["lateral"],
+        "symmetric_cbar": True,
+    },
     filename=f"{output_dir}/pow_diff_.png",
 )
 power.save(
     pvalues,
     mask_file="MNI152_T1_8mm_brain.nii.gz",
     parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz",
-    asymmetric_data=True,
-        plot_kwargs={
+    plot_kwargs={
         "views": ["lateral"],
         "cmap": "hot_r",
         "bg_on_data": 1,
         "darkness": 0.6,
-        "alpha": 1
+        "alpha": 1,
+        "vmin": 0,
+        "vmax": 0.1,
     },
     filename=f"{output_dir}/pow_diff_pvalues_.png",
 )

osl_dynamics/analysis/power.py

@@ -94,6 +94,7 @@ def variance_from_spectra(
     psd,
     components=None,
     frequency_range=None,
+    method="mean",
 ):
     """Calculates variance from power spectra.
 
@@ -109,6 +110,10 @@ def variance_from_spectra(
         Spectral components. Shape must be (n_components, n_freq).
     frequency_range : list, optional
         Frequency range in Hz to integrate the PSD over. Default is full range.
+    method : str
+        Should take the sum of the PSD over the frequency range
+        (:code:`method="sum"`) or take the average value of the PSD
+        (:code:`method="mean"`).
 
     Returns
     -------
@@ -164,6 +169,9 @@ def variance_from_spectra(
             "If frequency_range is passed, frequenices must also be passed."
         )
 
+    if method not in ["mean", "sum"]:
+        raise ValueError("method should be 'mean' or 'sum'.")
+
     # Number of spectral components
     if components is None:
         n_components = 1
@@ -195,10 +203,16 @@ def variance_from_spectra(
         else:
             # Integrate over the given frequency range
             if frequency_range is None:
-                p = np.sum(psd, axis=-1)
+                if method == "sum":
+                    p = np.sum(psd, axis=-1)
+                else:
+                    p = np.mean(psd, axis=-1)
             else:
                 [min_arg, max_arg] = get_frequency_args_range(f, frequency_range)
-                p = np.sum(psd[..., min_arg:max_arg], axis=-1)
+                if method == "sum":
+                    p = np.sum(psd[..., min_arg:max_arg], axis=-1)
+                else:
+                    p = np.mean(psd[..., min_arg:max_arg], axis=-1)
 
         p = p.reshape(n_components, n_modes, n_channels)
         var.append(p)

osl_dynamics/analysis/spectral.py

@@ -1117,8 +1117,7 @@ def get_frequency_args_range(frequencies, frequency_range):
     f_max_arg = np.argwhere(frequencies <= frequency_range[1])[-1, 0]
     if f_max_arg <= f_min_arg:
         raise ValueError("Cannot select requested frequency range.")
-    args_range = [f_min_arg, f_max_arg + 1]
-    return args_range
+    return [f_min_arg, f_max_arg]
 
 
 def _welch_spectrogram(