From 563c9c0ba09ac11302b8abb539e78eb0cb07c966 Mon Sep 17 00:00:00 2001 From: cgohil8 <53237863+cgohil8@users.noreply.github.com> Date: Wed, 28 Jun 2023 09:09:45 +0100 Subject: [PATCH] Remove: OSL_HMM/HMM_MAR and updated examples. (#158) --- examples/meg/uk_meg_notts/dynemo_fixed-cov.py | 92 --------- examples/meg/uk_meg_notts/hmm_tinda.py | 32 +-- examples/meg/uk_meg_notts/sage_fixed-cov.py | 77 -------- examples/plotting/activity-maps_means.py | 29 --- examples/plotting/activity-maps_tde-cov.py | 47 ----- examples/plotting/activity-maps_vars.py | 30 --- examples/plotting/activity-maps_workbench.py | 35 ---- ...nnectivity-maps_corr-abs.py => aec_abs.py} | 35 +++- ...nnectivity-maps_corr-sep.py => aec_sep.py} | 36 ++-- ...ectra.py => coherence_networks_spectra.py} | 13 +- ... coherence_networks_spectra_glassbrain.py} | 12 +- ....py => coherence_networks_spectra_nnmf.py} | 12 +- .../plotting/connectivity-maps_tde-cov.py | 59 ------ examples/plotting/mean_activity_maps.py | 38 ++++ .../plotting/mean_activity_maps_workbench.py | 44 +++++ ...-maps_spectra.py => power_maps_spectra.py} | 16 +- examples/plotting/power_maps_tde_cov.py | 53 +++++ osl_dynamics/data/__init__.py | 3 +- osl_dynamics/data/osl.py | 187 ------------------ 19 files changed, 220 insertions(+), 630 deletions(-) delete mode 100644 examples/meg/uk_meg_notts/dynemo_fixed-cov.py delete mode 100644 examples/meg/uk_meg_notts/sage_fixed-cov.py delete mode 100644 examples/plotting/activity-maps_means.py delete mode 100644 examples/plotting/activity-maps_tde-cov.py delete mode 100644 examples/plotting/activity-maps_vars.py delete mode 100644 examples/plotting/activity-maps_workbench.py rename examples/plotting/{connectivity-maps_corr-abs.py => aec_abs.py} (56%) rename examples/plotting/{connectivity-maps_corr-sep.py => aec_sep.py} (57%) rename examples/plotting/{connectivity-maps_spectra.py => coherence_networks_spectra.py} (78%) rename examples/plotting/{connectivity-maps_spectra_glassbrain.py => coherence_networks_spectra_glassbrain.py} (87%) rename examples/plotting/{connectivity-maps_spectra-nnmf.py => coherence_networks_spectra_nnmf.py} (84%) delete mode 100644 examples/plotting/connectivity-maps_tde-cov.py create mode 100644 examples/plotting/mean_activity_maps.py create mode 100644 examples/plotting/mean_activity_maps_workbench.py rename examples/plotting/{activity-maps_spectra.py => power_maps_spectra.py} (79%) create mode 100644 examples/plotting/power_maps_tde_cov.py delete mode 100644 osl_dynamics/data/osl.py diff --git a/examples/meg/uk_meg_notts/dynemo_fixed-cov.py b/examples/meg/uk_meg_notts/dynemo_fixed-cov.py deleted file mode 100644 index 9ea4c94a..00000000 --- a/examples/meg/uk_meg_notts/dynemo_fixed-cov.py +++ /dev/null @@ -1,92 +0,0 @@ -"""Example script for running inference on resting-state MEG data for one subject. - -- The data is stored on the BMRC cluster: /well/woolrich/projects/uk_meg_notts -- Uses the final covariances inferred by an HMM fit from OSL for the covariance of each - mode. -- Covariances are NOT trainable. -- Achieves a dice coefficient of ~0.94 (when compared to the OSL HMM mode time course). -""" - -print("Setting up") -from osl_dynamics.data import OSL_HMM, Data, processing -from osl_dynamics.inference import metrics, modes, tf_ops -from osl_dynamics.models.dynemo import Config, Model - -# GPU settings -tf_ops.gpu_growth() - -# Settings -config = Config( - n_modes=6, - n_channels=80, - sequence_length=400, - inference_n_units=64, - inference_normalization="layer", - model_n_units=64, - model_normalization="layer", - theta_normalization="layer", - learn_alpha_temperature=True, - initial_alpha_temperature=1.0, - learn_means=False, - learn_covariances=False, - do_kl_annealing=True, - kl_annealing_curve="tanh", - kl_annealing_sharpness=10, - n_kl_annealing_epochs=100, - batch_size=64, - learning_rate=0.01, - n_epochs=200, -) - -# Read MEG data -print("Reading MEG data") -prepared_data = Data( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/prepared_data/subject1.mat", - sampling_frequency=250, - n_embeddings=15, -) - -# Prepare dataset -training_dataset = prepared_data.dataset( - config.sequence_length, - config.batch_size, - shuffle=True, -) -prediction_dataset = prepared_data.dataset( - config.sequence_length, - config.batch_size, - shuffle=False, -) - -# Initialise covariances with the final HMM covariances -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/results/Subj1-1_K-6/hmm.mat" -) -config.initial_covariances = hmm.covariances - -# Build model -model = Model(config) -model.summary() - -print("Training model") -history = model.fit( - training_dataset, - epochs=config.n_epochs, - save_best_after=config.n_kl_annealing_epochs, - save_filepath="tmp/model", -) - -# Free energy = Log Likelihood - KL Divergence -free_energy = model.free_energy(prediction_dataset) -print(f"Free energy: {free_energy}") - -# Inferred mode mixing factors and mode time courses -alpha = model.get_alpha(prediction_dataset) -inf_stc = modes.argmax_time_courses(alpha) -hmm_stc = processing.trim_time_series( - time_series=hmm.mode_time_course(), - sequence_length=config.sequence_length, -) - -# Dice coefficient -print("Dice coefficient:", metrics.dice_coefficient(hmm_stc, inf_stc)) diff --git a/examples/meg/uk_meg_notts/hmm_tinda.py b/examples/meg/uk_meg_notts/hmm_tinda.py index ed345964..b5407c81 100644 --- a/examples/meg/uk_meg_notts/hmm_tinda.py +++ b/examples/meg/uk_meg_notts/hmm_tinda.py @@ -3,28 +3,38 @@ """ import os +import pickle import numpy as np import matplotlib.pyplot as plt from scipy.stats import ttest_rel -from osl_dynamics.inference.modes import argmax_time_courses +from osl_dynamics.inference import modes from osl_dynamics.analysis.tinda import tinda, plot_cycle, optimise_sequence -from osl_dynamics.data import HMM_MAR -# Some settings +# Settings os.makedirs("figures", exist_ok=True) # Directory for plots do_bonferroni_correction = True -# Load precomputed HMM -print("Loading HMM") -hmm = HMM_MAR( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/results/Subj1-55_K-12/hmm.mat" -) +# We will download example data hosted on osf.io/by2tc. +# Note, osfclient must be installed. This can be installed with pip: +# +# pip install osfclient + +def get_data(name, output_dir): + if os.path.exists(output_dir): + print(f"{output_dir} already downloaded. Skipping..") + return + os.system(f"osf -p by2tc fetch trained_models/{name}.zip") + os.system(f"unzip -o {name}.zip -d {output_dir}") + os.remove(f"{name}.zip") + print(f"Data downloaded to: {output_dir}") + +get_data("tde-hmm_notts_rest_55_subj", output_dir="notts_tde_hmm") -# Get state time courses -gamma = hmm.gamma() -stc = argmax_time_courses(gamma) +# Get the state time course +alpha = pickle.load(open("notts_tde_hmm/alpha.pkl", "rb")) +stc = modes.argmax_time_courses(alpha) # Run TINDA on every subject print("Running TINDA") diff --git a/examples/meg/uk_meg_notts/sage_fixed-cov.py b/examples/meg/uk_meg_notts/sage_fixed-cov.py deleted file mode 100644 index 402f130e..00000000 --- a/examples/meg/uk_meg_notts/sage_fixed-cov.py +++ /dev/null @@ -1,77 +0,0 @@ -"""Example script for running inference on resting-state MEG data for one subject. - -- The data is stored on the BMRC cluster: /well/woolrich/projects/uk_meg_notts -- Uses the final covariances inferred by an HMM fit from OSL for the covariance of each - mode. -- Covariances are NOT trainable. -- Achieves a dice coefficient of ~0.94 (when compared to the OSL HMM mode time course). -""" - -print("Setting up") -from osl_dynamics.data import OSL_HMM, Data, processing -from osl_dynamics.inference import metrics, modes, tf_ops -from osl_dynamics.models.sage import Config, Model - -# GPU settings -tf_ops.gpu_growth() - -# Settings -config = Config( - n_modes=6, - n_channels=80, - sequence_length=400, - inference_n_units=64, - inference_normalization="layer", - model_n_units=64, - model_normalization="layer", - discriminator_n_units=16, - discriminator_normalization="layer", - learn_means=False, - learn_covariances=False, - batch_size=64, - learning_rate=0.01, - n_epochs=200, -) - -# Read MEG data -print("Reading MEG data") -prepared_data = Data( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/prepared_data/subject1.mat", - sampling_frequency=250, - n_embeddings=15, -) - -# Prepare dataset -training_dataset = prepared_data.dataset( - config.sequence_length, - config.batch_size, - shuffle=True, -) -prediction_dataset = prepared_data.dataset( - config.sequence_length, - config.batch_size, - shuffle=False, -) - -# Initialise covariances with the final HMM covariances -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/results/Subj1-1_K-6/hmm.mat" -) -config.initial_covariances = hmm.covariances - -# Build model -model = Model(config) - -print("Training model") -history = model.fit(training_dataset) - -# Inferred mode mixing factors and mode time courses -alpha = model.get_alpha(prediction_dataset) -inf_stc = modes.argmax_time_courses(alpha) -hmm_stc = processing.trim_time_series( - time_series=hmm.mode_time_course(), - sequence_length=config.sequence_length, -) - -# Dice coefficient -print("Dice coefficient:", metrics.dice_coefficient(hmm_stc, inf_stc)) diff --git a/examples/plotting/activity-maps_means.py b/examples/plotting/activity-maps_means.py deleted file mode 100644 index c1b0dfef..00000000 --- a/examples/plotting/activity-maps_means.py +++ /dev/null @@ -1,29 +0,0 @@ -"""Example script for plotting the activity maps from the state/mode means. - -In this example we plot the state means from an HMM fit, but this can be -easily substituted with mode means from DyNeMo. -""" - -from osl_dynamics.analysis import power -from osl_dynamics.data import OSL_HMM - -# Load an HMM trained on amplitude envelope data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/summer21/results/Subj1-55_K-8_zeroMean-0/hmm.mat" -) - -# Get the inferred state means -means = hmm.means - -# Files used to source reconstruct the data the HMM was trained on -mask_file = "MNI152_T1_8mm_brain.nii.gz" -parcellation_file = "fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz" - -# Mean activity maps -power.save( - means, - filename="maps_.png", - mask_file=mask_file, - parcellation_file=parcellation_file, - subtract_mean=True, -) diff --git a/examples/plotting/activity-maps_tde-cov.py b/examples/plotting/activity-maps_tde-cov.py deleted file mode 100644 index 2ccf14e2..00000000 --- a/examples/plotting/activity-maps_tde-cov.py +++ /dev/null @@ -1,47 +0,0 @@ -"""Example code for plotting power maps from the state/mode covariances inferred -on time-delay embedded training data. - -In this example we use the inferred matrix directly. -""" - -from osl_dynamics.analysis import power, modes -from osl_dynamics.data import OSL_HMM, rw - -# Load an HMM trained on time-delay embedded data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/results/Subj1-10_K-6/hmm.mat" -) - -# Get the inferred covariances -# These covariances describe the training data in the time-delay embedded/PCA space -# We must convert this matrix into the source space by reversing the PCA and -# time-delay embedding -covs = hmm.covariances - -# This is the number of embeddings and PCA components used to prepare the training data -n_embeddings = 15 -pca_components = rw.loadmat( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/prepared_data/pca_components.mat" -) - -# Files used to create the source reconstructed data -mask_file = "MNI152_T1_8mm_brain.nii.gz" -parcellation_file = ( - "fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz" -) - -# Convert from the time-delay embedded/PCA space to the original source space -power_map = modes.raw_covariances( - covs, - n_embeddings=n_embeddings, - pca_components=pca_components, -) - -# Save the power maps -power.save( - power_map=power_map, - filename="maps_.png", - mask_file=mask_file, - parcellation_file=parcellation_file, - subtract_mean=True, -) diff --git a/examples/plotting/activity-maps_vars.py b/examples/plotting/activity-maps_vars.py deleted file mode 100644 index 4b7d2ffd..00000000 --- a/examples/plotting/activity-maps_vars.py +++ /dev/null @@ -1,30 +0,0 @@ -"""Example script for plotting the activity maps from the state/mode covariances. - -In this example our measure of activity is the variance (i.e. power) which we -extract from the diagonal of the covariance matrix directly. This approach is -typically used when we train with zero mean. -""" - -from osl_dynamics.analysis import power -from osl_dynamics.data import OSL_HMM - -# Load an HMM trained on amplitude envelope data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/summer21/results/Subj1-55_K-8_zeroMean-1/hmm.mat" -) - -# Get the inferred state covariances -covs = hmm.covariances - -# Files used during source reconstruction -mask_file = "MNI152_T1_8mm_brain.nii.gz" -parcellation_file = "fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz" - -# Plot power maps (this function will automatically extract the diagonal) -power.save( - covs, - filename="maps_.png", - mask_file=mask_file, - parcellation_file=parcellation_file, - subtract_mean=True, -) diff --git a/examples/plotting/activity-maps_workbench.py b/examples/plotting/activity-maps_workbench.py deleted file mode 100644 index 93aee19c..00000000 --- a/examples/plotting/activity-maps_workbench.py +++ /dev/null @@ -1,35 +0,0 @@ -"""Example code for saving maps as nii files and plotting with workbench. - -""" - -from osl_dynamics.analysis import power, workbench -from osl_dynamics.data import OSL_HMM - -# Load an HMM trained on amplitude envelope data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/summer21/results/Subj1-55_K-6_zeroMean-0/hmm.mat" -) - -# Mean activity maps (state means) -means = hmm.means - -# Source reconstruction files used to create the data the HMM was trained on -mask_file = "MNI152_T1_8mm_brain.nii.gz" -parcellation_file = "fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz" - -# Save nii file -power.save( - power_map=means, - filename="maps.nii.gz", - mask_file=mask_file, - parcellation_file=parcellation_file, - subtract_mean=True, -) - -# Setup workbench by specifying the path it is located -workbench.setup("/well/woolrich/projects/software/workbench/bin_rh_linux64") - -# Use workbench to save the maps in the nii file as images -# "tmp" is a directory used to store temporary files created by workbench -# this directory can be safely deleted after the maps have been saved -workbench.render("maps.nii.gz", "tmp", gui=False, image_name="maps_.png") diff --git a/examples/plotting/connectivity-maps_corr-abs.py b/examples/plotting/aec_abs.py similarity index 56% rename from examples/plotting/connectivity-maps_corr-abs.py rename to examples/plotting/aec_abs.py index a001c9fa..2e458999 100644 --- a/examples/plotting/connectivity-maps_corr-abs.py +++ b/examples/plotting/aec_abs.py @@ -7,31 +7,46 @@ a DyNeMo fit. """ +import os +import numpy as np + from osl_dynamics.array_ops import cov2corr from osl_dynamics.analysis import connectivity -from osl_dynamics.data import OSL_HMM -# Load an HMM trained on amplitude envelope data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/summer21/results/Subj1-55_K-8_zeroMean-1/hmm.mat" -) +# We will download example data hosted on osf.io/by2tc. +# Note, osfclient must be installed. This can be installed with pip: +# +# pip install osfclient + +def get_data(name, output_dir): + if os.path.exists(output_dir): + print(f"{output_dir} already downloaded. Skipping..") + return + os.system(f"osf -p by2tc fetch trained_models/{name}.zip") + os.system(f"unzip -o {name}.zip -d {output_dir}") + os.remove(f"{name}.zip") + print(f"Data downloaded to: {output_dir}") + +get_data("ae-hmm_notts_rest_55_subj", output_dir="notts_ae_hmm") -# Get the inferred state covariances -covariances = hmm.covariances +# Load the inferred state covariances +covariances = np.load("notts_ae_hmm/covs.npy") # Convert the covariance matrix to a correlation matrix correlations = cov2corr(covariances) +# Take the absolute value +conn_map = abs(correlations) + # We have many options for how to threshold the maps: # - Plot the top X % of connections. # - Plot the top X % of connections relative to the means across states/modes. # - Use a GMM to calculate a threshold (X) for us. # # Here we will just plot the top 2% of the absolute values -conn_map = abs(correlations) connectivity.save( - connectivity_map=conn_map, + conn_map, threshold=0.98, - filename="corr_.png", parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz", + filename="corr_.png", ) diff --git a/examples/plotting/connectivity-maps_corr-sep.py b/examples/plotting/aec_sep.py similarity index 57% rename from examples/plotting/connectivity-maps_corr-sep.py rename to examples/plotting/aec_sep.py index ff0beb17..02e11662 100644 --- a/examples/plotting/connectivity-maps_corr-sep.py +++ b/examples/plotting/aec_sep.py @@ -7,26 +7,34 @@ a DyNeMo fit. """ +import os import numpy as np from osl_dynamics.array_ops import cov2corr from osl_dynamics.analysis import connectivity -from osl_dynamics.data import OSL_HMM -# Load an HMM trained on amplitude envelope data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/summer21/results/Subj1-55_K-8_zeroMean-1/hmm.mat" -) +# We will download example data hosted on osf.io/by2tc. +# Note, osfclient must be installed. This can be installed with pip: +# +# pip install osfclient + +def get_data(name, output_dir): + if os.path.exists(output_dir): + print(f"{output_dir} already downloaded. Skipping..") + return + os.system(f"osf -p by2tc fetch trained_models/{name}.zip") + os.system(f"unzip -o {name}.zip -d {output_dir}") + os.remove(f"{name}.zip") + print(f"Data downloaded to: {output_dir}") -# Get the inferred state covariances -covariances = hmm.covariances +get_data("ae-hmm_notts_rest_55_subj", output_dir="notts_ae_hmm") + +# Load the inferred state covariances +covariances = np.load("notts_ae_hmm/covs.npy") # Convert the covariance matrix to a correlation matrix correlations = cov2corr(covariances) -# Files used to source reconstruct the data the HMM was trained on -parcellation_file = "fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz" - # We have many options for how to threshold the maps: # - Plot the top X % of connections. # - Plot the top X % of connections relative to the means across states/modes. @@ -39,12 +47,12 @@ ) pos_conn_map, neg_conn_map = connectivity.separate_edges(conn_map) connectivity.save( - connectivity_map=pos_conn_map, + pos_conn_map, + parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz", filename="pos_corr_.png", - parcellation_file=parcellation_file, ) connectivity.save( - connectivity_map=neg_conn_map, + neg_conn_map, + parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz", filename="neg_corr_.png", - parcellation_file=parcellation_file, ) diff --git a/examples/plotting/connectivity-maps_spectra.py b/examples/plotting/coherence_networks_spectra.py similarity index 78% rename from examples/plotting/connectivity-maps_spectra.py rename to examples/plotting/coherence_networks_spectra.py index cea1484b..46604f9f 100644 --- a/examples/plotting/connectivity-maps_spectra.py +++ b/examples/plotting/coherence_networks_spectra.py @@ -4,8 +4,8 @@ The spectra can be calculate with a multitaper (in the case of a state time course) or regression (in the case of a mode time course). -See examples/analysis/multitaper_spectra.py for how to calculate a multitaper -and examples/analysis/regression_spectra.py for how to calculate a regression. +See examples/minimal/multitaper_spectra.py for how to calculate a multitaper +and examples/minimal/regression_spectra.py for how to calculate regression spectra. In this script we assume this has been done and we have the group-average spectra files: f.npy and coh.npy. (The other file: psd.npy is not needed for this script.) @@ -15,11 +15,6 @@ from osl_dynamics.analysis import connectivity -# Source reconstruction files used to create the training data -parcellation_file = ( - "fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz" -) - # Load subject-specific state/mode spectra f = np.load("f.npy") coh = np.load("coh.npy") @@ -44,7 +39,7 @@ # Plot connectivity maps connectivity.save( - connectivity_map=conn_map, + conn_map, + parcellation_file="fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz", filename="coh_.png", - parcellation_file=parcellation_file, ) diff --git a/examples/plotting/connectivity-maps_spectra_glassbrain.py b/examples/plotting/coherence_networks_spectra_glassbrain.py similarity index 87% rename from examples/plotting/connectivity-maps_spectra_glassbrain.py rename to examples/plotting/coherence_networks_spectra_glassbrain.py index 2b359f7a..03dac90d 100644 --- a/examples/plotting/connectivity-maps_spectra_glassbrain.py +++ b/examples/plotting/coherence_networks_spectra_glassbrain.py @@ -5,7 +5,7 @@ course) or regression (in the case of a mode time course). See examples/analysis/multitaper_spectra.py for how to calculate a multitaper -and examples/analysis/regression_spectra.py for how to calculate a regression. +and examples/analysis/regression_spectra.py for how to calculate regression spectra. In this script we assume this has been done and we have the subject-specific spectra: - f.npy, the frequency axis. @@ -19,10 +19,6 @@ from osl_dynamics.analysis import connectivity -# Source reconstruction files used to create the training data -parcellation_file = ( - "fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz" -) # Load subject-specific spectra f = np.load("f.npy") @@ -48,8 +44,8 @@ # Plot connectivity maps connectivity.save( - connectivity_map=conn_map, - filename="coh_.html", - parcellation_file=parcellation_file, + conn_map, + parcellation_file="fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz", glassbrain=True, + filename="coh_.html", ) diff --git a/examples/plotting/connectivity-maps_spectra-nnmf.py b/examples/plotting/coherence_networks_spectra_nnmf.py similarity index 84% rename from examples/plotting/connectivity-maps_spectra-nnmf.py rename to examples/plotting/coherence_networks_spectra_nnmf.py index 7cdd61ff..674fdefd 100644 --- a/examples/plotting/connectivity-maps_spectra-nnmf.py +++ b/examples/plotting/coherence_networks_spectra_nnmf.py @@ -5,7 +5,6 @@ course) or regression (in the case of a mode time course). See examples/minimal/multitaper_spectra.py for how to calculate a multitaper -and examples/minimal/regression_spectra.py for how to calculate a regression. """ import numpy as np @@ -13,11 +12,6 @@ from osl_dynamics.analysis import connectivity, spectral from osl_dynamics.utils import plotting -# Source reconstruction files used to create the training data -parcellation_file = ( - "fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz" -) - # Load subject-specific state/mode spectra f = np.load("f.npy") coh = np.load("coh.npy") @@ -56,8 +50,8 @@ # Plot connectivity maps connectivity.save( - connectivity_map=conn_map, - filename="coh_.png", - parcellation_file=parcellation_file, + conn_map, + parcellation_file="fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz", component=0, # only plot the first spectral component + filename="coh_.png", ) diff --git a/examples/plotting/connectivity-maps_tde-cov.py b/examples/plotting/connectivity-maps_tde-cov.py deleted file mode 100644 index 175e6465..00000000 --- a/examples/plotting/connectivity-maps_tde-cov.py +++ /dev/null @@ -1,59 +0,0 @@ -"""Example code for plotting connectivity based on a correlation matrix. - -In this example we use the inferred state/mode covariance matrices from -training on time-embedded/PCA data to calculate the state/mode correlation -matrices. - -In this script we use an HMM fit, but this can be easily substituted for -a DyNeMo fit. -""" - -from osl_dynamics.array_ops import cov2corr -from osl_dynamics.analysis import connectivity, modes -from osl_dynamics.data import OSL_HMM, rw - -# Load an HMM trained on time-delay embedded data -hmm = OSL_HMM( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/results/Subj1-10_K-6/hmm.mat" -) - -# Get the inferred covariances -# These covariances describe the training data in the time-delay embedded/PCA space -# We must convert this matrix into the source space by reversing the PCA and -# time-delay embedding -covs = hmm.covariances - -# This is the number of embeddings and PCA components used to prepare the training data -n_embeddings = 15 -pca_components = rw.loadmat( - "/well/woolrich/projects/uk_meg_notts/eo/natcomms18/prepared_data/pca_components.mat" -) - -# Files used to create the source reconstructed data -parcellation_file = ( - "fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz" -) - -# Convert from the time-delay embedded/PCA space to the original source space -covs = modes.raw_covariances( - covs, - n_embeddings=n_embeddings, - pca_components=pca_components, - zero_lag=False, # Should we just consider the zero-lag covariance? -) - -# Convert from covariance to correlation -corrs = cov2corr(covs) - -# We have many options for how to threshold the maps: -# - Plot the top X % of connections. -# - Plot the top X % of connections relative to the means across states/modes. -# - Use a GMM to calculate a threshold (X) for us. -# -# Here we will just plot the top 2% of the absolute values -connectivity.save( - connectivity_map=abs(corrs), - filename="corr_.png", - threshold=0.98, - parcellation_file=parcellation_file, -) diff --git a/examples/plotting/mean_activity_maps.py b/examples/plotting/mean_activity_maps.py new file mode 100644 index 00000000..0d1a3438 --- /dev/null +++ b/examples/plotting/mean_activity_maps.py @@ -0,0 +1,38 @@ +"""Example script for plotting the activity maps from the state/mode means. + +In this example we plot the state means from an HMM fit, but this can be +easily substituted with mode means from DyNeMo. +""" + +import os +import numpy as np + +from osl_dynamics.analysis import power + +# We will download example data hosted on osf.io/by2tc. +# Note, osfclient must be installed. This can be installed with pip: +# +# pip install osfclient + +def get_data(name, output_dir): + if os.path.exists(output_dir): + print(f"{output_dir} already downloaded. Skipping..") + return + os.system(f"osf -p by2tc fetch trained_models/{name}.zip") + os.system(f"unzip -o {name}.zip -d {output_dir}") + os.remove(f"{name}.zip") + print(f"Data downloaded to: {output_dir}") + +get_data("ae-hmm_notts_rest_55_subj", output_dir="notts_ae_hmm") + +# Load the inferred state means +means = np.load("notts_ae_hmm/means.npy") + +# Mean activity maps +power.save( + means, + mask_file="MNI152_T1_8mm_brain.nii.gz", + parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz", + subtract_mean=True, + filename="maps_.png", +) diff --git a/examples/plotting/mean_activity_maps_workbench.py b/examples/plotting/mean_activity_maps_workbench.py new file mode 100644 index 00000000..5f60e4e1 --- /dev/null +++ b/examples/plotting/mean_activity_maps_workbench.py @@ -0,0 +1,44 @@ +"""Example code for saving maps as nii files and plotting with workbench. + +""" + +import os +import numpy as np + +from osl_dynamics.analysis import power, workbench + +# We will download example data hosted on osf.io/by2tc. +# Note, osfclient must be installed. This can be installed with pip: +# +# pip install osfclient + +def get_data(name, output_dir): + if os.path.exists(output_dir): + print(f"{output_dir} already downloaded. Skipping..") + return + os.system(f"osf -p by2tc fetch trained_models/{name}.zip") + os.system(f"unzip -o {name}.zip -d {output_dir}") + os.remove(f"{name}.zip") + print(f"Data downloaded to: {output_dir}") + +get_data("ae-hmm_notts_rest_55_subj", output_dir="notts_ae_hmm") + +# Load mean activity maps (state means) +means = np.load("notts_ae_hmm/means.npy") + +# Save nii file +power.save( + means, + mask_file="MNI152_T1_8mm_brain.nii.gz", + parcellation_file="fmri_d100_parcellation_with_PCC_reduced_2mm_ss5mm_ds8mm.nii.gz", + subtract_mean=True, + filename="maps.nii.gz", +) + +# Setup workbench by specifying the path it is located +workbench.setup("/well/woolrich/projects/software/workbench/bin_rh_linux64") + +# Use workbench to save the maps in the nii file as images +# "tmp" is a directory used to store temporary files created by workbench +# this directory can be safely deleted after the maps have been saved +workbench.render("maps.nii.gz", "tmp", gui=False, image_name="maps_.png") diff --git a/examples/plotting/activity-maps_spectra.py b/examples/plotting/power_maps_spectra.py similarity index 79% rename from examples/plotting/activity-maps_spectra.py rename to examples/plotting/power_maps_spectra.py index f4a596ed..0080c436 100644 --- a/examples/plotting/activity-maps_spectra.py +++ b/examples/plotting/power_maps_spectra.py @@ -5,7 +5,7 @@ course) or regression (in the case of a mode time course). See examples/minimal/multitaper_spectra.py for how to calculate a multitaper -and examples/minimal/regression_spectra.py for how to calculate a regression. +and examples/minimal/regression_spectra.py for how to calculate regression spectra. In this script we assume this has been done and we have the subject-specific spectra files: @@ -28,21 +28,15 @@ # Calculate the group average gpsd = np.average(psd, axis=0, weights=w) -# Source reconstruction files used to create the training data -mask_file = "MNI152_T1_8mm_brain.nii.gz" -parcellation_file = ( - "fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz" -) - # Calculate power maps from the spectra # (frequency_range is an optional argument) power_map = power.variance_from_spectra(f, gpsd, frequency_range=[1, 30]) # Save the power maps as images power.save( - power_map=power_map, - filename="maps_.png", - mask_file=mask_file, - parcellation_file=parcellation_file, + power_map, + mask_file="MNI152_T1_8mm_brain.nii.gz", + parcellation_file="fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz", subtract_mean=True, + filename="maps_.png", ) diff --git a/examples/plotting/power_maps_tde_cov.py b/examples/plotting/power_maps_tde_cov.py new file mode 100644 index 00000000..d898012a --- /dev/null +++ b/examples/plotting/power_maps_tde_cov.py @@ -0,0 +1,53 @@ +"""Example code for plotting power maps from the state/mode covariances inferred +on time-delay embedded training data. + +""" + +import os +import numpy as np + +from osl_dynamics.analysis import power, modes +from osl_dynamics.data import rw + +# We will download example data hosted on osf.io/by2tc. +# Note, osfclient must be installed. This can be installed with pip: +# +# pip install osfclient + +def get_data(name, output_dir): + if os.path.exists(output_dir): + print(f"{output_dir} already downloaded. Skipping..") + return + os.system(f"osf -p by2tc fetch trained_models/{name}.zip") + os.system(f"unzip -o {name}.zip -d {output_dir}") + os.remove(f"{name}.zip") + print(f"Data downloaded to: {output_dir}") + +get_data("tde-hmm_notts_rest_55_subj", output_dir="notts_tde_hmm") + +# Get the inferred covariances +# +# These covariances describe the training data in the time-delay embedded/PCA space +# We must convert this matrix into the source space by reversing the PCA and +# time-delay embedding +covs = np.load("notts_tde_hmm/covs.npy") + +# This is the number of embeddings and PCA components used to prepare the training data +n_embeddings = 15 +pca_components = rw.loadmat("notts_tde_hmm/pca_components.mat") + +# Convert from the time-delay embedded/PCA space to the original source space +power_map = modes.raw_covariances( + covs, + n_embeddings=n_embeddings, + pca_components=pca_components, +) + +# Save the power maps +power.save( + power_map, + mask_file="MNI152_T1_8mm_brain.nii.gz", + parcellation_file="fmri_d100_parcellation_with_3PCC_ips_reduced_2mm_ss5mm_ds8mm_adj.nii.gz", + subtract_mean=True, + filename="maps_.png", +) diff --git a/osl_dynamics/data/__init__.py b/osl_dynamics/data/__init__.py index 11b93523..ff52445b 100644 --- a/osl_dynamics/data/__init__.py +++ b/osl_dynamics/data/__init__.py @@ -3,6 +3,5 @@ """ from osl_dynamics.data.base import Data -from osl_dynamics.data.osl import OSL_HMM, HMM_MAR -__all__ = ["Data", "OSL_HMM", "HMM_MAR"] +__all__ = ["Data"] diff --git a/osl_dynamics/data/osl.py b/osl_dynamics/data/osl.py deleted file mode 100644 index 16755f84..00000000 --- a/osl_dynamics/data/osl.py +++ /dev/null @@ -1,187 +0,0 @@ -"""Class to read HMM-MAR objects from MATLAB. - -""" - -import numpy as np -from osl_dynamics import array_ops -from osl_dynamics.data.processing import trim_time_series -from osl_dynamics.data.rw import loadmat - - -def OSL_HMM(filename): - """Wrapper for osl_dynamics.data.osl.HMM_MAR. - - Parameters - ---------- - filename : str - The location of the HMM-MAR object saved as a mat7.3 file. - """ - return HMM_MAR(filename) - - -class HMM_MAR: - """Imports and encapsulates HMM-MAR objects as a python class. - - Parameters - ---------- - filename : str - The location of the HMM-MAR object saved as a mat7.3 file. - """ - - def __init__(self, filename): - self.filename = filename - hmm = loadmat(filename) - self.hmm = hmm["hmm"] if "hmm" in hmm else hmm - - self.state = self.hmm["state"] - self.mode = self.hmm["state"] - self.k = int(self.hmm["K"]) - self.p = self.hmm["P"] - self.dir_2d_alpha = self.hmm["Dir2d_alpha"] - self.pi = self.hmm["Pi"] - self.dir_alpha = self.hmm["Dir_alpha"] - self.prior = self.hmm["prior"] - self.train = self.hmm["train"] - - # State probabilities - self.Gamma = None - for gamma in ["gamma", "Gamma"]: - if gamma in self.hmm: - self.Gamma = self.hmm[gamma].astype(np.float32) - elif gamma in hmm: - self.Gamma = hmm[gamma].astype(np.float32) - - # State time course - if self.Gamma is not None: - vpath = self.Gamma.argmax(axis=1) - self.vpath = array_ops.get_one_hot(vpath).astype(np.float32) - else: - self.vpath = None - - # State means - self.means = np.array([state["W"]["Mu_W"] for state in self.state]) - - # State covariances - self.covariances = np.array( - [ - state["Omega"]["Gam_rate"] - / (state["Omega"]["Gam_shape"] - len(state["Omega"]["Gam_rate"]) - 1) - for state in self.state - ] - ) - - # Transition probability matrix - self.trans_prob = self.p - - # Discontinuities in the training data which indicate the number of data - # points for different subjects - if "T" in self.hmm: - self.discontinuities = [np.squeeze(T).astype(int) for T in self.hmm["T"]] - elif self.Gamma is not None: - # Assume gamma has no discontinuities - self.discontinuities = [self.Gamma.shape[0]] - else: - self.discontinuities = None - - def __str__(self): - return f"OSL HMM object from file {self.filename}" - - def gamma(self, concatenate=False, pad=None): - """State probabilities for each subject. - - Parameters - ---------- - concatenate : bool - Should we concatenate the gammas for each subjects? - pad : int - Pad the gamma for each subject with zeros to replace the data points lost - by performing n_embeddings. Default is no padding. - - Returns - ------- - gamma : np.ndarray or list - State probabilities. - """ - if self.Gamma is None: - return None - - if pad is None: - if concatenate or len(self.discontinuities) == 1: - return self.Gamma - else: - return np.split(self.Gamma, np.cumsum(self.discontinuities[:-1])) - else: - padded_gamma = [ - np.pad(gamma, [[pad, pad], [0, 0]]) - for gamma in np.split(self.Gamma, np.cumsum(self.discontinuities[:-1])) - ] - if concatenate: - return np.concatenate(padded_gamma) - else: - return padded_gamma - - def trimmed_gamma(self, sequence_length, concatenate=False): - """Trimmed state probabilities for each subject. - - Data points that would be lost due to separating the time series into - sequences are removed from each subject. - - Parameters - ---------- - sequence_length : int - Sequence length. - concatenate : bool - Should we concatenate the gammas for each subject? - - Returns - ------- - gamma : np.ndarray or list - State probabilities. - """ - return trim_time_series( - self.gamma(), sequence_length=sequence_length, concatenate=concatenate - ) - - def state_time_course(self, concatenate=False, pad=None): - """State time course for each subject. - - Parameters - ---------- - concatenate : bool - Should we concatenate the state time course for each subjects? - pad : int - Pad the state time course for each subject with zeros to replace the data - points lost by performing n_embeddings. Default is no padding. - - Returns - ------- - stc : np.ndarray or list - State time course. - """ - if self.vpath is None: - return None - - if pad is None: - if concatenate or len(self.discontinuities) == 1: - return self.vpath - else: - return np.split(self.vpath, np.cumsum(self.discontinuities[:-1])) - else: - padded_stc = [ - np.pad(self.vpath, [[pad, pad], [0, 0]]) - for stc in np.split(self.vpath, np.cumsum(self.discontinuities[:-1])) - ] - if concatenate: - return np.concatenate(padded_stc) - else: - return padded_stc - - def mode_time_course(self, *args, **kwargs): - """Wrapper for the state_time_course method. - - Returns - ------- - mtc : np.ndarray or list - Mode time course. - """ - return self.state_time_course(*args, **kwargs)