From d0919572e8b9bd1e369f399eebbc01f66a983915 Mon Sep 17 00:00:00 2001 From: John Wang Date: Mon, 9 Jun 2025 14:23:52 +0000 Subject: [PATCH 1/2] example for analysis using model --- examples/eg__analysis_model.py | 178 +++++++++++++++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 examples/eg__analysis_model.py diff --git a/examples/eg__analysis_model.py b/examples/eg__analysis_model.py new file mode 100644 index 00000000..0838892e --- /dev/null +++ b/examples/eg__analysis_model.py @@ -0,0 +1,178 @@ +""" +.. _ex-modelanalysis: + +======================================================== +Modelling TMS-EEG evoked responses +======================================================== + +This example shows the analysis from model: + +1. model parameters +2. networks +3. neural states + +""" +# %% +# First we must import the necessary packages required for the example: + +# System-based packages +import os +import sys +sys.path.append('..') + + +# Whobpyt modules taken from the whobpyt package +import whobpyt +from whobpyt.datatypes import Parameter as par, Timeseries +from whobpyt.models.jansen_rit import JansenRitModel,JansenRitParams +from whobpyt.run import ModelFitting +from whobpyt.optimization.custom_cost_JR import CostsJR +from whobpyt.datasets.fetchers import fetch_egtmseeg + +# Python Packages used for processing and displaying given analytical data (supported for .mat and Google Drive files) +import numpy as np +import pandas as pd +import scipy.io +import gdown +import pickle +import warnings +warnings.filterwarnings('ignore') +import matplotlib.pyplot as plt # Plotting library (For Visualization) +import seaborn as sns + +import mne # Neuroimaging package + + + + + + +# %% +# load in a previously completed model fitting results object +full_run_fname = os.path.join(data_dir, 'Subject_1_low_voltage_fittingresults_stim_exp.pkl') +F = pickle.load(open(full_run_fname, 'rb')) + + +### get labels for Yeo 200 +url = 'https://raw.githubusercontent.com/ThomasYeoLab/CBIG/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/Parcellations/MNI/Centroid_coordinates/Schaefer2018_200Parcels_7Networks_order_FSLMNI152_2mm.Centroid_RAS.csv' +atlas = pd.read_csv(url) +labels = atlas['ROI Name'] + +# get networks +nets = [label.split('_')[2] for label in labels] +net_names = np.unique(np.array(nets)) + +# %% +# 1. model parameters +# ----------------------------------- + +# %% +# Plots of parameter values over Training (check if converges) +fig, axs = plt.subplots(2,2, figsize = (12,8)) +paras = ['c1', 'c2', 'c3', 'c4'] +for i in range(len(paras)): + axs[i//2,i%2].plot(F.trainingStats.fit_params[paras[i]]) + axs[i//2, i%2].set_title(paras[i]) +plt.title("Select Variables Changing Over Training Epochs") + +# %% +# Plots of parameter values over Training (prior vs post) +fig, axs = plt.subplots(2,2, figsize = (12,8)) +paras = ['c1', 'c2', 'c3', 'c4'] +for i in range(len(paras)): + axs[i//2,i%2].hist(F.trainingStats.fit_params[paras[i]][:500], label='prior') + axs[i//2,i%2].hist(F.trainingStats.fit_params[paras[i]][-500:], label='post') + axs[i//2, i%2].set_title(paras[i]) +plt.title("Prior vs Post") + +# %% +# 2. Networks +# ----------------------------------- +fig, axs = plt.subplots(1,3, figsize = (12,8)) +networks_frommodels = ['p2p', 'p2e', 'p2i'] +sns.heatmap(F.model.w_p2p.detach().numpy(), cmap = 'bwr', center=0, ax=axs[0]) +axs[0].set_title(networks_frommodels[0]) +sns.heatmap(F.model.w_p2p.detach().numpy(), cmap = 'bwr', center=0, ax=axs[1]) +axs[1].set_title(networks_frommodels[1]) +sns.heatmap(F.model.w_p2p.detach().numpy(), cmap = 'bwr', center=0, ax=axs[2]) +axs[2].set_title(networks_frommodels[2]) + + +# %% +# 3. Neural states +# ----------------------------------- + +#### plot E response on each networks +fig, ax = plt.subplots(2,4, figsize=(12,10), sharey= True) +t = np.linspace(-0.1,0.3, 400) + +for i, net in enumerate(net_names): + mask = np.array(nets) == net + ax[i//4, i%4].plot(t, F.lastRec['E'].npTS()[mask,:].mean(0).T) + ax[i//4, i%4].set_title(net) +plt.suptitle('Test: E') +plt.show() + +### plot I response at each networks +fig, ax = plt.subplots(2,4, figsize=(12,10), sharey= True) +t = np.linspace(-0.1,0.3, 400) + +for i, net in enumerate(net_names): + mask = np.array(nets) == net + ax[i//4, i%4].plot(t, F.lastRec['I'].npTS()[mask,:].mean(0).T) + ax[i//4, i%4].set_title(net) +plt.suptitle('Test: I') +plt.show() + +### plot P response at each networks +fig, ax = plt.subplots(2,4, figsize=(12,10), sharey= True) +t = np.linspace(-0.1,0.3, 400) + +for i, net in enumerate(net_names): + mask = np.array(nets) == net + ax[i//4, i%4].plot(t, F.lastRec['P'].npTS()[mask,:].mean(0).T) + ax[i//4, i%4].set_title(net) +plt.suptitle('Test: P') +plt.show() + + +### plot phase of E at each network +j = complex(0,1) +fig, ax = plt.subplots(2,4, figsize=(12,10), sharey= True) +t = np.linspace(-0.1,0.3, 400) + +phase = np.angle(F.lastRec['E'].npTS()+j*F.lastRec['Ev'].npTS()) +for i, net in enumerate(net_names): + mask = np.array(nets) == net + ax[i//4, i%4].plot(t, phase[mask,:].mean(0).T) + ax[i//4, i%4].set_title(net) +plt.suptitle('Test: phase E') +plt.show() + +### plot I phase at each network +j = complex(0,1) +fig, ax = plt.subplots(2,4, figsize=(12,10), sharey= True) +t = np.linspace(-0.1,0.3, 400) + +phase = np.angle(F.lastRec['I'].npTS()+j*F.lastRec['Iv'].npTS()) +for i, net in enumerate(net_names): + mask = np.array(nets) == net + ax[i//4, i%4].plot(t, phase[mask,:].mean(0).T) + ax[i//4, i%4].set_title(net) +plt.suptitle('Test: phase I') +plt.show() + +### plot P phase at each network + +j = complex(0,1) +fig, ax = plt.subplots(2,4, figsize=(12,10), sharey= True) +t = np.linspace(-0.1,0.3, 400) + +phase = np.angle(F.lastRec['P'].npTS()+j*F.lastRec['Pv'].npTS()) +for i, net in enumerate(net_names): + mask = np.array(nets) == net + ax[i//4, i%4].plot(t, phase[mask,:].mean(0).T) + ax[i//4, i%4].set_title(net) +plt.suptitle('Test: phase P') +plt.show() + From 3ddf2a242f194f32628ded1aaf5e31639abd0885 Mon Sep 17 00:00:00 2001 From: John Wang Date: Mon, 14 Jul 2025 16:36:42 +0000 Subject: [PATCH 2/2] move delay out of step loop --- whobpyt/models/jansen_rit/jansen_rit.py | 26 ++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/whobpyt/models/jansen_rit/jansen_rit.py b/whobpyt/models/jansen_rit/jansen_rit.py index 67029554..09161dd2 100644 --- a/whobpyt/models/jansen_rit/jansen_rit.py +++ b/whobpyt/models/jansen_rit/jansen_rit.py @@ -437,27 +437,27 @@ def forward(self, external, hx, hE): # Run through the number of specified sample points for this window for i_window in range(self.TRs_per_window): + # Collect the delayed inputs: + + # i) index the history of E + Ed = pttranspose(hE.clone().gather(1,self.delays), 0, 1) + + # ii) multiply the past states by the connectivity weights matrix, and sum over rows + LEd_p2e = ptsum(w_n_f * Ed, 1) + LEd_p2i = -ptsum(w_n_b * Ed, 1) + LEd_p2p = ptsum(w_n_l * Ed, 1) + # iii) reshape for next step + LEd_p2e = ptreshape(LEd_p2e, (n_nodes, 1)) + LEd_p2i = ptreshape(LEd_p2i, (n_nodes, 1)) + LEd_p2p = ptreshape(LEd_p2p, (n_nodes, 1)) # For each sample point, run the model by solving the differential # equations for a defined number of integration steps, # and keep only the final activity state within this set of steps for step_i in range(self.steps_per_TR): - # Collect the delayed inputs: - - # i) index the history of E - Ed = pttranspose(hE.clone().gather(1,self.delays), 0, 1) - - # ii) multiply the past states by the connectivity weights matrix, and sum over rows - LEd_p2e = ptsum(w_n_f * Ed, 1) - LEd_p2i = -ptsum(w_n_b * Ed, 1) - LEd_p2p = ptsum(w_n_l * Ed, 1) - # iii) reshape for next step - LEd_p2e = ptreshape(LEd_p2e, (n_nodes, 1)) - LEd_p2i = ptreshape(LEd_p2i, (n_nodes, 1)) - LEd_p2p = ptreshape(LEd_p2p, (n_nodes, 1)) # iv) if specified, add the laplacian component (self-connections from diagonals) if self.use_laplacian: