Include the code for data processing in the same tutorial

Author: htwangtw

content/gcn_decoding.md

@@ -68,65 +68,118 @@ print(X.shape)
 So we have 1452 time points, with one cognitive annotations each, and for each time point we have recordings of fMRI activity across 675 voxels. We can also see that the cognitive annotations span 9 different categories.
 
 ```{code-cell} ipython3
-proc_path = os.path.join(data_dir, 'haxby_proc/')
-concat_path = os.path.join(data_dir, 'haxby_concat/')
 conn_path = os.path.join(data_dir, 'haxby_connectomes/')
-split_path = os.path.join(data_dir, 'haxby_split_win/')
-```
+if not os.path.exists(conn_path):
+    os.makedirs(conn_path)
 
-```{code-cell} ipython3
 import glob
 import nilearn.connectome
 import numpy as np
 import warnings
 warnings.filterwarnings(action='once')
-# Estimating connectomes
+# Estimating connectomes and save for pytorch to load
 corr_measure = nilearn.connectome.ConnectivityMeasure(kind="correlation")
 conn = corr_measure.fit_transform([X])[0]
-
-sub_no = 4
 np.save(os.path.join(conn_path, 'conn_subj{}.npy'.format(sub_no)), conn)
+```
 
-conn_files = sorted(glob.glob(conn_path + '/*.npy')) 
+```{code-cell} ipython3
+concat_path = os.path.join(data_dir, 'haxby_concat/')
+if not os.path.exists(concat_path):
+    os.makedirs(concat_path)
+
+concat_bold_files = []
+for i in range(0,len(y)):
+    label = y[i]
+    concat_bold_files = X[i:i+1]
+    concat_file_name = concat_path + '{}_concat_fMRI.npy'.format(label)
+    
+    if os.path.isfile(concat_file_name):
+        concat_file = np.load(concat_file_name, allow_pickle = True)
+        concat_file = np.concatenate((concat_file, concat_bold_files), axis = 0)
+        np.save(concat_file_name, concat_file)
+    else:
+        np.save(concat_file_name, concat_bold_files)
+
+
+
+import pandas as pd
+split_path = os.path.join(data_dir, 'haxby_split_win/')
+if not os.path.exists(split_path):
+    os.makedirs(split_path)
 
+dic_labels = {'rest':0,'face':1,'chair':2,'scissors':3,'shoe':4,'scrambledpix':5,'house':6,'cat':7,'bottle':8}
+label_df = pd.DataFrame(columns=['label', 'filename'])
+processed_bold_files = sorted(glob.glob(concat_path + '/*.npy'))
+window_length = 1
+out_file = os.path.join(split_path, '{}_{:04d}.npy')
+out_csv = os.path.join(split_path, 'labels.csv')
+
+for proc_bold in processed_bold_files:
+    
+    ts_data = np.load(proc_bold)
+    ts_duration = len(ts_data)
+
+    ts_filename = os.path.basename(proc_bold)
+    ts_label = ts_filename.split('_', 1)[0]
+
+    valid_label = dic_labels[ts_label]
+    
+    # Split the timeseries
+    rem = ts_duration % window_length
+    n_splits = int(np.floor(ts_duration / window_length))
+
+    ts_data = ts_data[:(ts_duration-rem), :]   
+    
+    for j, split_ts in enumerate(np.split(ts_data, n_splits)):
+        ts_output_file_name = out_file.format(ts_filename, j)
+
+        split_ts = np.swapaxes(split_ts, 0, 1)
+        np.save(ts_output_file_name, split_ts)
+        curr_label = {'label': valid_label, 'filename': os.path.basename(ts_output_file_name)}
+        label_df = label_df.append(curr_label, ignore_index=True)
+    
+label_df.to_csv(out_csv, index=False)  
+```
+
+```{code-cell} ipython3
 # split dataset
 import sys
-sys.path.append(os.path.join(".."))
 sys.path.append('../src')
-
-import gcn_windows_dataset
+from gcn_windows_dataset import TimeWindowsDataset
 
 random_seed = 0
 
-train_dataset = gcn_windows_dataset.TimeWindowsDataset(
-    data_dir=split_path
-    , partition="train"
-    , random_seed=random_seed
-    , pin_memory=True
-    , normalize=True,
-    shuffle = True)
-
-valid_dataset = gcn_windows_dataset.TimeWindowsDataset(
-    data_dir=split_path
-    , partition="valid"
-    , random_seed=random_seed
-    , pin_memory=True
-    , normalize=True, 
-    shuffle = True)
-
-test_dataset = gcn_windows_dataset.TimeWindowsDataset(
-    data_dir=split_path
-    , partition="test"
-    , random_seed=random_seed
-    , pin_memory=True
-    , normalize=True, 
-    shuffle = True)
+train_dataset = TimeWindowsDataset(
+    data_dir=split_path, 
+    partition="train", 
+    random_seed=random_seed, 
+    pin_memory=True, 
+    normalize=True,
+    shuffle=True)
+
+valid_dataset = TimeWindowsDataset(
+    data_dir=split_path, 
+    partition="valid", 
+    random_seed=random_seed, 
+    pin_memory=True, 
+    normalize=True,
+    shuffle=True)
+
+test_dataset = TimeWindowsDataset(
+    data_dir=split_path, 
+    partition="test", 
+    random_seed=random_seed, 
+    pin_memory=True, 
+    normalize=True,
+    shuffle=True)
 
 print("train dataset: {}".format(train_dataset))
 print("valid dataset: {}".format(valid_dataset))
 print("test dataset: {}".format(test_dataset))
+```
 
-
+```{code-cell} ipython3
 import torch
 
 torch.manual_seed(random_seed)
@@ -136,10 +189,6 @@ test_generator = torch.utils.data.DataLoader(test_dataset, batch_size=16, shuffl
 train_features, train_labels = next(iter(train_generator))
 print(f"Feature batch shape: {train_features.size()}; mean {torch.mean(train_features)}")
 print(f"Labels batch shape: {train_labels.size()}; mean {torch.mean(torch.Tensor.float(train_labels))}")
-
-connectomes = []
-for conn_file in conn_files:
-      connectomes += [np.load(conn_file)]
 ```
 
 ## Building brain graphs
@@ -152,11 +201,14 @@ Each node is only connected to *k* other neighbours, which is __8 nodes__ with t
 For more details you please check out __*src/graph_construction.py*__ script.
 
 ```{code-cell} ipython3
-import graph_construction
+conn_files = sorted(glob.glob(conn_path + '/*.npy')) 
+connectomes = []
+for conn_file in conn_files:
+      connectomes += [np.load(conn_file)]
 
+from graph_construction import make_group_graph
 
-graph = graph_construction.make_group_graph(connectomes, self_loops=False, 
-                                            k=8, symmetric=True)
+graph = make_group_graph(connectomes, self_loops=False, k=8, symmetric=True)
 ```
 
 ## Running model
@@ -172,11 +224,10 @@ In this example we will continue with __*window_length = 1*__, which means each
 TR is cycle time between corresponding points in fMRI.
 
 ```{code-cell} ipython3
-import gcn_model
+from gcn_model import GCN
 
 window_length = 1
-gcn = gcn_model.GCN(graph.edge_index, graph.edge_attr, 
-                           n_timepoints=window_length)
+gcn = GCN(graph.edge_index, graph.edge_attr, n_timepoints=window_length)
 gcn
 ```
 

src/graph_construction.py

@@ -3,7 +3,7 @@
 import torch_geometric as tg
 
 
-def make_undirected(mat):
+def _make_undirected(mat):
     """
     Takes an input adjacency matrix and makes it undirected (symmetric).
 
@@ -21,7 +21,7 @@ def make_undirected(mat):
     return sym
 
 
-def knn_graph_quantile(mat, self_loops=False, k=8, symmetric=True):
+def _knn_graph_quantile(mat, self_loops=False, k=8, symmetric=True):
     """
     Takes an input correlation matrix and returns a k-Nearest
     Neighbour weighted undirected adjacency matrix.
@@ -46,7 +46,7 @@ def knn_graph_quantile(mat, self_loops=False, k=8, symmetric=True):
     if not self_loops:
         np.fill_diagonal(adj, 0)
     if symmetric:
-        adj = make_undirected(adj)
+        adj = _make_undirected(adj)
     return adj
 
 
@@ -75,7 +75,7 @@ def make_group_graph(connectomes, k=8, self_loops=False, symmetric=True):
     # Group average connectome and nndirected 8 k-NN graph
     avg_conn = np.array(connectomes).mean(axis=0)
     avg_conn = np.round(avg_conn, 6)
-    avg_conn_k = knn_graph_quantile(avg_conn, k=k, self_loops=self_loops, symmetric=symmetric)
+    avg_conn_k = _knn_graph_quantile(avg_conn, k=k, self_loops=self_loops, symmetric=symmetric)
 
     # Format matrix into graph for torch_geometric
     adj_sparse = tg.utils.dense_to_sparse(torch.from_numpy(avg_conn_k))