init commit

LICENSE.txt

@@ -0,0 +1,32 @@
+The Clear BSD License
+
+Copyright (c) 2022, authors of DeepSIF
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted (subject to the limitations in the disclaimer
+below) provided that the following conditions are met:
+
+     * Redistributions of source code must retain the above copyright notice,
+     this list of conditions and the following disclaimer.
+
+     * Redistributions in binary form must reproduce the above copyright
+     notice, this list of conditions and the following disclaimer in the
+     documentation and/or other materials provided with the distribution.
+
+     * Neither the name of the copyright holder nor the names of its
+     contributors may be used to endorse or promote products derived from this
+     software without specific prior written permission.
+
+NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
+THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
+CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
+IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.

README.md

@@ -1 +1,95 @@
-Placeholder directory for: Deep Neural Networks Constrained by Neural Mass Models Improve Electrophysiological Source Imaging of Spatio-temporal Brain Dynamics
+# DeepSIF: Deep Learning based Source Imaging Framework
+
+
+DeepSIF is an EEG/MEG source imaging framework aiming at providing an estimation of the location, size, and temporal activity of the brain activities from scalp EEG/MEG recordings. There are three components: training data generation (```forward/```), neural network training (```main.py```), and model evaluation (```eval_sim.py```,```eval_recal.py```), as detailed below. The codes are provided as a service to the scientific community, and should be used at users’ own risks.
+
+
+This work was supported in part by the National Institutes of Health grants NS096761, EB021027, AT009263, MH114233, EB029354, and NS124564, awarded to Dr. Bin He, Carnegie Mellon University. Additional data in 20 human epilepsy patients tested in this work can be found at
+
+
+https://figshare.com/s/580622eaf17108da49d7.
+
+
+Please cite the following publication if you are using any part of the codes:
+
+Sun R, Sohrabpour A, Worrell GA, He B: “Deep Neural Networks Constrained by Neural Mass Models Improve Electrophysiological Source Imaging of Spatio-temporal Brain Dynamics.” Proceedings of the National Academy of Sciences, 2022.
+
+
+
+## Train Data Generation
+#### The Virtual Brain Simulation
+```bash
+python generate_tvb_data.py --a_start 0 --a_end 10
+```
+The simulation for each region can also run in parallel. (Require multiprocessing installed.)
+
+#### Process Raw TVB Data and Prepare Training/Testing Dataset 
+Run in Matlab
+```matlab
+process_raw_nmm
+generate_sythetic_source
+```
+The output of ```generate_sythetic_source``` can be used as input for ```loaders.SpikeEEGBuild``` or ```loaders.SpikeEEGBuildEval```
+
+## Training
+
+After sythetic training dataset is generated, ```main.py``` can be used to train a DeepSIF model. ```network.py``` contains the architecture used
+ in the paper.  ```loaders.py``` provides two ways to load the dataset. If the data is already saved in seperate input/output files
+ , ```SpikeEEGLoad``` can be used. If training data is generated on the run, ```SpikeEEGBuild``` can be used to generate different types of
+  training data. To train a model, use 
+
+```bash
+python main.py --model_id 1
+```
+**Parameters:**
+
+    '--save', type=int, default=True, help='save each epoch or not'
+    '--workers', default=0, type=int, help='number of data loading workers'
+    '--batch_size', default=64, type=int, help='batch size'
+    '--device', default='cuda:0', type=str, help='device running the code'
+    '--arch', default='TemporalInverseNet', type=str, help='network achitecture class'
+    '--dat', default='SpikeEEGBuild', type=str, help='data loader class'
+    '--train', default='test_sample_source2.mat', type=str, help='train dataset name or directory'
+    '--test', default='test_sample_source2.mat', type=str, help='test dataset name or directory'
+    '--model_id', default=75, type=int, help='model id'
+    '--lr', default=3e-4, type=float, help='learning rate'
+    '--resume', default='1', type=str, help='epoch id to resume'
+    '--epoch', default=20, type=int, help='total number of epoch'
+    '--fwd', default='leadfield_75_20k.mat', type=str, help='forward matrix to use'
+    '--rnn_layer', default=3, type=int, help='number of rnn layer'
+    '--info', default='', type=str, help='other information regarding this model'
+
+## Evaluation
+
+#### Simulation :
+After a model is trained, ```eval_sim.py``` can be used to evaluate the trained model in simulations under different conditions. Some examples are:
+```bash
+python eval_sim.py --model_id 75
+```
+Additional tests: narrow-band input
+```bash
+python eval_sim.py --model_id 75 --lfreq 1 --hfreq 3
+```
+Additional tests: different noise type
+```bash
+python eval_sim.py --model_id 75 --snr_rsn_ratio 0.5
+```
+Additional tests: different head / conductivity value / electrode locations
+```bash
+python eval_sim.py --model_id 75 --fwd <the forward matrix file>
+```
+#### Real data :
+Or use real data as the model input as shown in ```eval_real.py```:
+```bash
+python eval_real.py
+```
+Default subject folder : VEP
+<!---Results visualization examples are in tutorials/-->
+
+## Dependencies
+* [Python >= 3.8.3](https://www.python.org/downloads/)
+* [PyTorch>=1.6.0](https://pytorch.org/)
+* [tvb](https://www.thevirtualbrain.org/tvb/zwei)
+* [mne](https://mne.tools/stable/index.html)
+* [h5py](https://www.h5py.org/)
+* [numpy](https://numpy.org/)

anatomy/README.md

@@ -0,0 +1,38 @@
+### This folder contains miscellaneous files related to the forward process
+#### For TVB simulation
+* connectivity_76.zip
+* connectivity_998.zip
+
+Connectivity file used for the TVB simulation, downloaded from
+https://github.com/the-virtual-brain/tvb-data
+
+#### Template headmodel:
+- **fs_cortex_20k.mat**: fsaverage5 cortex
+    - pos, tri, ori: vertices positions, triangulations, and dipole orientations; resolution 20484
+    - left_ind, left_tri: index for left hemisphere, triangulations in the left hemisphere
+    - right_ind, right_tri: index for right hemisphere, triangulations in the right hemisphere
+
+- **fs_cortex_20k_inflated.mat**: inflated fsaverage5 cortex
+    - pos, tri: vertices positions, triangulations; resolution 20484
+    - posl, tril: vertices positions and triangulations in the left hemisphere
+    - posr, trir: vertices positions and triangulations in the right hemisphere
+
+- **fs_cortex_20k_region_mapping.mat** : map fsaverage5 to 994 regions
+    - rm: region mapping id, size 1*20484
+    - nbs: neighbours for each region
+
+- **leadfield_75_20k.mat** : leadfield matrix for fsaverage5, 75 channels
+    - fwd: size 75*994
+
+- **dis_matrix_fs_20k.mat** : distance between source centres
+    - raw_dis_matrix: size 994*994
+
+- **electrode_75.mat** : 75 EEG channels in EEGLAB format
+    - eloc75
+
+- **fsaverage5/**: contains the files for the raw freesurfer output, for plotting in mne    
+
+#### Simulations:
+- **realistic_noise.mat** : resting data extracted from 75 channel EEG recordings
+    - data: num_examples * num_time* num_channel; 4*500*75
+    - npower: the power for each channel; 4*75 

eval_real.py

@@ -0,0 +1,80 @@
+import argparse
+import os
+import time
+from scipy.io import loadmat, savemat
+import numpy as np
+import glob
+
+import torch
+import network
+
+
+def main():
+    start_time = time.time()
+    # parse the input
+    parser = argparse.ArgumentParser(description='DeepSIF Model')
+    parser.add_argument('--device', default='cpu', type=str, help='device running the code')
+    parser.add_argument('--model_id', type=int, default=64, help='model id')
+    parser.add_argument('--resume', default='', type=str, help='epoch id to resume')
+    parser.add_argument('--info', default='', type=str, help='other information regarding this model')
+    args = parser.parse_args()
+
+    # ======================= PREPARE PARAMETERS =====================================================================================================
+    use_cuda = (False) and torch.cuda.is_available()  # Only use GPU during training
+    device = torch.device(args.device if use_cuda else "cpu")
+    result_root = 'model_result/{}_the_model'.format(args.model_id)
+    if not os.path.exists(result_root):
+        print("ERROR: No model {}".format(args.model_id))
+        return
+
+    # =============================== LOAD MODEL =====================================================================================================
+    if args.resume:
+        fn = fn = os.path.join(result_root, 'epoch_' + args.resume)
+    else:
+        fn = os.path.join(result_root, 'model_best.pth.tar')
+    print("=> Load checkpoint", fn)
+    if os.path.isfile(fn):
+        print("=> Found checkpoint '{}'".format(fn))
+        checkpoint = torch.load(fn, map_location=torch.device('cpu'))
+        best_result = checkpoint['best_result']
+        net = network.__dict__[checkpoint['arch']](*checkpoint['attribute_list']).to(device)  # redefine the weights architecture
+        net.load_state_dict(checkpoint['state_dict'], strict=False)
+        print("=> Loaded checkpoint {}, current results: {}".format(fn, best_result))
+    else:
+        print("ERROR: no checkpoint found")
+        return
+
+    print('Number of parameters:', net.count_parameters())
+    print('Prepare time:', time.time() - start_time)
+
+    # =============================== EVALUATION =====================================================================================================
+    net.eval()
+    subject_list = ['VEP']
+    for pii in subject_list:
+        folder_name = 'source/{}'.format(pii)
+        start_time = time.time()
+        flist = glob.glob(folder_name + '/data*.mat')
+        if len(flist) == 0:
+            print('WARNING: NO FILE IN FOLDER {}.'.format(folder_name))
+            continue
+        flist = sorted(flist, key=lambda name: int(os.path.basename(name)[4:-4]))  # sort file based on nature number
+        test_data = []
+        for i in flist:
+            data = loadmat(i)['data']
+            # data = data - np.mean(data, 0, keepdims=True)
+            # data = data - np.mean(data, 1, keepdims=True)
+            data = data / np.max(np.abs(data[:]))
+            test_data.append(data)
+
+        data = torch.from_numpy(np.array(test_data)).to(device, torch.float)
+        out = net(data)['last']
+        # calculate the loss
+        all_out = out.detach().cpu().numpy()
+        # visualize the result in Matlab
+        savemat(folder_name + '/rnn_test_{}_{}.mat'.format(args.model_id, fn[-8:]), {'all_out': all_out})
+        print('Save output as:', folder_name + '/rnn_test_{}_{}.mat'.format(args.model_id, fn[-8:]))
+    print('Total run time:', time.time() - start_time)
+
+if __name__ == '__main__':
+    main()
+