+[docs]
+def main(*args):
""" Parse arguments and estimate model
"""
@@ -1620,6 +1687,13 @@
+normative_parallel — Predictive Clinical Neuroscience Toolkit 0.20 documentation
-
-
-
+
+
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
@@ -62,13 +58,24 @@
Source code for normative
exec(func + '(' + all_args + ')')
+[docs]
+def entrypoint():
+ main(sys.argv[1:])
+
+
+
# For running from the command line:
if __name__ == "__main__":
main(sys.argv[1:])
diff --git a/doc/build/html/_modules/normative_parallel.html b/doc/build/html/_modules/normative_parallel.html
index 1cb6ed98..bcd334f9 100644
--- a/doc/build/html/_modules/normative_parallel.html
+++ b/doc/build/html/_modules/normative_parallel.html
@@ -1,28 +1,24 @@
+
+
-
+
Tutorials
-
-
- Gaussian Process Regression -
- Hierarchical Bayesian Regression +
- DEMO ON NORMATIVE MODELING +
- Predictive Clinical Neuroscience Toolkit +
- Hierarchical Bayesian Regression Normative Modelling and Transfer onto unseen site.
- Braincharts: transfer -
- Bayesian Linear Regression +
- Predictive Clinical Neuroscience Toolkit +
- The Normative Modeling Framework for Computational Psychiatry Protocol +
- Data Preparation +
- Algorithm & Modeling +
- Evaluation & Interpretation
- Visualization of normative modeling outputs
- Post-hoc analysis on normative modeling outputs -
- Predictive modeling using deviation scores +
- Load Data +
- Create Train/Test Splits +
- Principal Component Regression (BBS) +
- Connectome Predictive Modelling +
- Lasso (Linear Regression + L1 Regularization) +
- Ridge (Linear Regression + L2 Regularization) +
- Elastic Net (Linear Regression + L1/L2 Regularization)
Other Useful Stuff
-
@@ -135,7 +142,8 @@
- Gaussian Process Regression -
- Hierarchical Bayesian Regression +
- DEMO ON NORMATIVE MODELING +
- Predictive Clinical Neuroscience Toolkit +
- Hierarchical Bayesian Regression Normative Modelling and Transfer onto unseen site.
- Braincharts: transfer -
- Bayesian Linear Regression +
- Predictive Clinical Neuroscience Toolkit +
- The Normative Modeling Framework for Computational Psychiatry Protocol +
- Data Preparation +
- Algorithm & Modeling +
- Evaluation & Interpretation
- Visualization of normative modeling outputs
- Post-hoc analysis on normative modeling outputs -
- Predictive modeling using deviation scores +
- Load Data +
- Create Train/Test Splits +
- Principal Component Regression (BBS) +
- Connectome Predictive Modelling +
- Lasso (Linear Regression + L1 Regularization) +
- Ridge (Linear Regression + L2 Regularization) +
- Elastic Net (Linear Regression + L1/L2 Regularization)
- Gaussian Process Regression -
- Hierarchical Bayesian Regression +
- DEMO ON NORMATIVE MODELING +
- Predictive Clinical Neuroscience Toolkit +
- Hierarchical Bayesian Regression Normative Modelling and Transfer onto unseen site.
- Braincharts: transfer -
- Bayesian Linear Regression +
- Predictive Clinical Neuroscience Toolkit +
- The Normative Modeling Framework for Computational Psychiatry Protocol +
- Data Preparation +
- Algorithm & Modeling +
- Evaluation & Interpretation
- Visualization of normative modeling outputs
- Post-hoc analysis on normative modeling outputs -
- Predictive modeling using deviation scores +
- Load Data +
- Create Train/Test Splits +
- Principal Component Regression (BBS) +
- Connectome Predictive Modelling +
- Lasso (Linear Regression + L1 Regularization) +
- Ridge (Linear Regression + L2 Regularization) +
- Elastic Net (Linear Regression + L1/L2 Regularization)
Source code for normative_parallel
import time
import numpy as np
import pandas as pd
-from subprocess import call, check_output
+from datetime import datetime
+from subprocess import run, check_output
try:
@@ -157,7 +165,9 @@ Source code for normative_parallel
PICKLE_PROTOCOL = configs.PICKLE_PROTOCOL
-[docs]def execute_nm(processing_dir,
+
+[docs]
+def execute_nm(processing_dir,
python_path,
job_name,
covfile_path,
@@ -192,6 +202,7 @@ Source code for normative_parallel
:param testrespfile_path: Full path to a .txt file that contains all test features
:param log_path: Path for saving log files
:param binary: If True uses binary format for response file otherwise it is text
+ :param cluster_spec: 'torque' for PBS Torque and 'slurm' for Slurm clusters.
:param interactive: If False (default) the user should manually
rerun the failed jobs or collect the results.
If 'auto' the job status are checked until all
@@ -211,10 +222,11 @@ Source code for normative_parallel
cv_folds = kwargs.get('cv_folds', None)
testcovfile_path = kwargs.get('testcovfile_path', None)
testrespfile_path = kwargs.get('testrespfile_path', None)
- outputsuffix = kwargs.get('outputsuffix', 'estimate')
+ outputsuffix = kwargs.get('outputsuffix', '_estimate')
cluster_spec = kwargs.pop('cluster_spec', 'torque')
log_path = kwargs.get('log_path', None)
binary = kwargs.pop('binary', False)
+ cores = kwargs.pop('n_cores_per_batch','1')
split_nm(processing_dir,
respfile_path,
@@ -234,6 +246,8 @@ Source code for normative_parallel
kwargs.update({'batch_size': str(batch_size)})
job_ids = []
+ start_time = datetime.now().strftime("%Y-%m-%dT%H:%M:%S")
+
for n in range(1, number_of_batches+1):
kwargs.update({'job_id': str(n)})
if testrespfile_path is not None:
@@ -266,9 +280,10 @@ Source code for normative_parallel
job_id = qsub_nm(job_path=batch_job_path,
log_path=log_path,
memory=memory,
- duration=duration)
+ duration=duration,
+ cores=cores)
job_ids.append(job_id)
- elif cluster_spec == 'sbatch':
+ elif cluster_spec == 'slurm':
# update the response file
kwargs.update({'testrespfile_path':
batch_testrespfile_path})
@@ -282,8 +297,10 @@ Source code for normative_parallel
memory=memory,
duration=duration,
**kwargs)
- sbatch_nm(job_path=batch_job_path,
- log_path=log_path)
+
+ job_id = sbatch_nm(job_path=batch_job_path)
+ job_ids.append(job_id)
+
elif cluster_spec == 'new':
# this part requires addition in different envioronment [
sbatchwrap_nm(processing_dir=batch_processing_dir,
@@ -310,9 +327,10 @@ Source code for normative_parallel
job_id = qsub_nm(job_path=batch_job_path,
log_path=log_path,
memory=memory,
- duration=duration)
+ duration=duration,
+ cores=cores)
job_ids.append(job_id)
- elif cluster_spec == 'sbatch':
+ elif cluster_spec == 'slurm':
sbatchwrap_nm(batch_processing_dir,
python_path,
normative_path,
@@ -323,8 +341,9 @@ Source code for normative_parallel
memory=memory,
duration=duration,
**kwargs)
- sbatch_nm(job_path=batch_job_path,
- log_path=log_path)
+
+ job_id = sbatch_nm(job_path=batch_job_path)
+ job_ids.append(job_id)
elif cluster_spec == 'new':
# this part requires addition in different envioronment [
bashwrap_nm(processing_dir=batch_processing_dir, func=func,
@@ -352,9 +371,10 @@ Source code for normative_parallel
job_id = qsub_nm(job_path=batch_job_path,
log_path=log_path,
memory=memory,
- duration=duration)
+ duration=duration,
+ cores=cores)
job_ids.append(job_id)
- elif cluster_spec == 'sbatch':
+ elif cluster_spec == 'slurm':
sbatchwrap_nm(batch_processing_dir,
python_path,
normative_path,
@@ -365,8 +385,10 @@ Source code for normative_parallel
memory=memory,
duration=duration,
**kwargs)
- sbatch_nm(job_path=batch_job_path,
- log_path=log_path)
+
+ job_id = sbatch_nm(job_path=batch_job_path)
+ job_ids.append(job_id)
+
elif cluster_spec == 'new':
# this part requires addition in different envioronment [
bashwrap_nm(processing_dir=batch_processing_dir, func=func,
@@ -376,7 +398,7 @@ Source code for normative_parallel
if interactive:
- check_jobs(job_ids, delay=60)
+ check_jobs(job_ids, cluster_spec, start_time, delay=60)
success = False
while (not success):
@@ -393,16 +415,33 @@ Source code for normative_parallel
if interactive == 'query':
response = yes_or_no('Rerun the failed jobs?')
if response:
- rerun_nm(processing_dir, log_path=log_path, memory=memory,
- duration=duration, binary=binary,
- interactive=interactive)
+ if cluster_spec == 'torque':
+ rerun_nm(processing_dir, log_path=log_path, memory=memory,
+ duration=duration, binary=binary,
+ interactive=interactive, cores=cores)
+ elif cluster_spec == 'slurm':
+ sbatchrerun_nm(processing_dir,
+ memory=memory,
+ duration=duration,
+ binary=binary,
+ log_path=log_path,
+ interactive=interactive)
+
else:
success = True
else:
print('Reruning the failed jobs ...')
- rerun_nm(processing_dir, log_path=log_path, memory=memory,
- duration=duration, binary=binary,
- interactive=interactive)
+ if cluster_spec == 'torque':
+ rerun_nm(processing_dir, log_path=log_path, memory=memory,
+ duration=duration, binary=binary,
+ interactive=interactive, cores=cores)
+ elif cluster_spec == 'slurm':
+ sbatchrerun_nm(processing_dir,
+ memory=memory,
+ duration=duration,
+ binary=binary,
+ log_path=log_path,
+ interactive=interactive)
if interactive == 'query':
response = yes_or_no('Collect the results?')
@@ -425,10 +464,13 @@ Source code for normative_parallel
outputsuffix=outputsuffix)
+
"""routines that are environment independent"""
-[docs]def split_nm(processing_dir,
+
+[docs]
+def split_nm(processing_dir,
respfile_path,
batch_size,
binary,
@@ -543,7 +585,10 @@ Source code for normative_parallel
protocol=PICKLE_PROTOCOL)
-[docs]def collect_nm(processing_dir,
+
+
+[docs]
+def collect_nm(processing_dir,
job_name,
func='estimate',
collect=False,
@@ -581,19 +626,15 @@ Source code for normative_parallel
batch_fail = []
if (func != 'fit' and func != 'extend' and func != 'merge' and func != 'tune'):
- file_example = []
# TODO: Collect_nm only depends on yhat, thus does not work when no
# prediction is made (when test cov is not specified).
- for batch in batches:
- if file_example == []:
- file_example = glob.glob(batch + 'yhat' + outputsuffix
+ files = glob.glob(processing_dir + 'batch_*/' + 'yhat' + outputsuffix
+ file_extentions)
- else:
- break
- if binary is False:
- file_example = fileio.load(file_example[0])
+ if len(files) > 0:
+ file_example = fileio.load(files[0])
else:
- file_example = pd.read_pickle(file_example[0])
+ raise ValueError(f"Missing output files (yhats at: {processing_dir + 'batch_*/' + 'yhat' + outputsuffix + file_extentions}")
+
numsubjects = file_example.shape[0]
try:
# doesn't exist if size=1, and txt file
@@ -947,7 +988,10 @@ Source code for normative_parallel
return False
-[docs]def delete_nm(processing_dir,
+
+
+[docs]
+def delete_nm(processing_dir,
binary=False):
'''This function deletes all processing for normative modelling and just keeps the combined output.
@@ -971,11 +1015,14 @@ Source code for normative_parallel
os.remove(processing_dir + 'failed_batches' + file_extentions)
+
# all routines below are envronment dependent and require adaptation in novel
# environments -> copy those routines and adapt them in accrodance with your
# environment
-[docs]def bashwrap_nm(processing_dir,
+
+[docs]
+def bashwrap_nm(processing_dir,
python_path,
normative_path,
job_name,
@@ -1036,7 +1083,7 @@ Source code for normative_parallel
job_call = [python_path + ' ' + normative_path + ' -c ' +
covfile_path + ' -f ' + func]
else:
- raise ValueError("""For 'estimate' function either testcov or cvfold
+ raise ValueError("""For 'estimate' function either testrespfile_path or cvfold
must be specified.""")
# add algorithm-specific parameters
@@ -1065,10 +1112,14 @@ Source code for normative_parallel
os.chmod(processing_dir + job_name, 0o770)
-[docs]def qsub_nm(job_path,
+
+
+[docs]
+def qsub_nm(job_path,
log_path,
memory,
- duration):
+ duration,
+ cores):
'''This function submits a job.sh scipt to the torque custer using the qsub command.
Basic usage::
@@ -1088,10 +1139,10 @@ Source code for normative_parallel
# created qsub command
if log_path is None:
qsub_call = ['echo ' + job_path + ' | qsub -N ' + job_path + ' -l ' +
- 'procs=1' + ',mem=' + memory + ',walltime=' + duration]
+ 'nodes=1:ppn='+ cores + ',mem=' + memory + ',walltime=' + duration]
else:
qsub_call = ['echo ' + job_path + ' | qsub -N ' + job_path +
- ' -l ' + 'procs=1' + ',mem=' + memory + ',walltime=' +
+ ' -l ' + 'nodes=1:ppn='+ cores + ',mem=' + memory + ',walltime=' +
duration + ' -o ' + log_path + ' -e ' + log_path]
# submits job to cluster
@@ -1102,10 +1153,15 @@ Source code for normative_parallel
return job_id
-[docs]def rerun_nm(processing_dir,
+
+
+[docs]
+def rerun_nm(processing_dir,
log_path,
memory,
duration,
+ cluster_spec,
+ cores,
binary=False,
interactive=False):
'''This function reruns all failed batched in processing_dir after collect_nm has identified the failed batches.
@@ -1133,7 +1189,8 @@ Source code for normative_parallel
job_id = qsub_nm(job_path=jobpath,
log_path=log_path,
memory=memory,
- duration=duration)
+ duration=duration,
+ cores=cores)
job_ids.append(job_id)
else:
file_extentions = '.txt'
@@ -1146,17 +1203,21 @@ Source code for normative_parallel
job_id = qsub_nm(job_path=jobpath,
log_path=log_path,
memory=memory,
- duration=duration)
+ duration=duration,
+ cores=cores)
job_ids.append(job_id)
if interactive:
- check_jobs(job_ids, delay=60)
+ check_jobs(job_ids, cluster_spec, delay=60)
+
# COPY the rotines above here and aadapt those to your cluster
# bashwarp_nm; qsub_nm; rerun_nm
-[docs]def sbatchwrap_nm(processing_dir,
+
+[docs]
+def sbatchwrap_nm(processing_dir,
python_path,
normative_path,
job_name,
@@ -1164,6 +1225,7 @@ Source code for normative_parallel
respfile_path,
memory,
duration,
+ log_path,
func='estimate',
**kwargs):
'''This function wraps normative modelling into a bash script to run it
@@ -1202,14 +1264,15 @@ Source code for normative_parallel
output_changedir = ['cd ' + processing_dir + '\n']
sbatch_init = '#!/bin/bash\n'
- sbatch_jobname = '#SBATCH --job-name=' + processing_dir + '\n'
- sbatch_account = '#SBATCH --account=p33_norment\n'
+ sbatch_jobname = '#SBATCH --job-name=' + job_name + '\n'
sbatch_nodes = '#SBATCH --nodes=1\n'
sbatch_tasks = '#SBATCH --ntasks=1\n'
sbatch_time = '#SBATCH --time=' + str(duration) + '\n'
sbatch_memory = '#SBATCH --mem-per-cpu=' + str(memory) + '\n'
- sbatch_module = 'module purge\n'
- sbatch_anaconda = 'module load anaconda3\n'
+ sbatch_log_out = '#SBATCH -o ' + log_path + '%x_%j.out' + '\n'
+ sbatch_log_error = '#SBATCH -e ' + log_path + '%x_%j.err' + '\n'
+ # sbatch_module = 'module purge\n'
+ # sbatch_anaconda = 'module load anaconda3\n'
sbatch_exit = 'set -o errexit\n'
# echo -n "This script is running on "
@@ -1217,12 +1280,13 @@ Source code for normative_parallel
bash_environment = [sbatch_init +
sbatch_jobname +
- sbatch_account +
sbatch_nodes +
sbatch_tasks +
sbatch_time +
- sbatch_module +
- sbatch_anaconda]
+ sbatch_memory +
+ sbatch_log_out +
+ sbatch_log_error
+ ]
# creates call of function for normative modelling
if (testrespfile_path is not None) and (testcovfile_path is not None):
@@ -1239,7 +1303,7 @@ Source code for normative_parallel
job_call = [python_path + ' ' + normative_path + ' -c ' +
covfile_path + ' -f ' + func]
else:
- raise ValueError("""For 'estimate' function either testcov or cvfold
+ raise ValueError("""For 'estimate' function either testrespfile_path or cv_folds
must be specified.""")
# add algorithm-specific parameters
@@ -1268,19 +1332,20 @@ Source code for normative_parallel
os.chmod(processing_dir + job_name, 0o770)
-[docs]def sbatch_nm(job_path,
- log_path):
+
+
+[docs]
+def sbatch_nm(job_path):
'''This function submits a job.sh scipt to the torque custer using the qsub
command.
Basic usage::
- sbatch_nm(job_path, log_path)
+ sbatch_nm(job_path)
:param job_path: Full path to the job.sh file
- :param log_path: The logs are currently stored in the working dir
- :outputs: Submission of the job to the (torque) cluster.
+ :outputs: Submission of the job to the slurm cluster.
written by (primarily) T Wolfers, (adapted) S Rutherford.
'''
@@ -1289,15 +1354,22 @@ Source code for normative_parallel
sbatch_call = ['sbatch ' + job_path]
# submits job to cluster
- call(sbatch_call, shell=True)
+ job_id = check_output(sbatch_call, shell=True).decode(
+ sys.stdout.encoding).replace("\n", "")
+
+ return job_id
-[docs]def sbatchrerun_nm(processing_dir,
+
+
+[docs]
+def sbatchrerun_nm(processing_dir,
memory,
duration,
new_memory=False,
new_duration=False,
binary=False,
+ interactive=False,
**kwargs):
'''This function reruns all failed batched in processing_dir after collect_nm has identified he failed batches.
@@ -1315,7 +1387,12 @@ Source code for normative_parallel
written by (primarily) T Wolfers, (adapted) S Rutherford.
'''
- log_path = kwargs.pop('log_path', None)
+
+ # log_path = kwargs.pop('log_path', None)
+
+ job_ids = []
+
+ start_time = datetime.now().strftime("%Y-%m-%dT%H:%M:%S")
if binary:
file_extentions = '.pkl'
@@ -1329,11 +1406,12 @@ Source code for normative_parallel
with fileinput.FileInput(jobpath, inplace=True) as file:
for line in file:
print(line.replace(duration, new_duration), end='')
- if new_memory != False:
- with fileinput.FileInput(jobpath, inplace=True) as file:
- for line in file:
- print(line.replace(memory, new_memory), end='')
- sbatch_nm(jobpath, log_path)
+ if new_memory != False:
+ with fileinput.FileInput(jobpath, inplace=True) as file:
+ for line in file:
+ print(line.replace(memory, new_memory), end='')
+ job_id = sbatch_nm(jobpath)
+ job_ids.append(job_id)
else:
file_extentions = '.txt'
@@ -1347,73 +1425,124 @@ Source code for normative_parallel
with fileinput.FileInput(jobpath, inplace=True) as file:
for line in file:
print(line.replace(duration, new_duration), end='')
- if new_memory != False:
- with fileinput.FileInput(jobpath, inplace=True) as file:
- for line in file:
- print(line.replace(memory, new_memory), end='')
- sbatch_nm(jobpath,
- log_path)
+ if new_memory != False:
+ with fileinput.FileInput(jobpath, inplace=True) as file:
+ for line in file:
+ print(line.replace(memory, new_memory), end='')
+ job_id = sbatch_nm(jobpath)
+ job_ids.append(job_id)
+
+ if interactive:
+ check_jobs(job_ids, cluster_spec='slurm',
+ start_time=start_time, delay=60)
+
-[docs]def retrieve_jobs():
+
+[docs]
+def retrieve_jobs(cluster_spec, start_time=None):
"""
A utility function to retrieve task status from the outputs of qstat.
+ :param cluster_spec: type of cluster, either 'torque' or 'slurm'.
+
:return: a dictionary of jobs.
"""
- output = check_output('qstat', shell=True).decode(sys.stdout.encoding)
- output = output.split('\n')
- jobs = dict()
- for line in output[2:-1]:
- (Job_ID, Job_Name, User, Wall_Time, Status, Queue) = line.split()
- jobs[Job_ID] = dict()
- jobs[Job_ID]['name'] = Job_Name
- jobs[Job_ID]['walltime'] = Wall_Time
- jobs[Job_ID]['status'] = Status
+ if cluster_spec == 'torque':
+
+ output = check_output('qstat', shell=True).decode(sys.stdout.encoding)
+ output = output.split('\n')
+ jobs = dict()
+ for line in output[2:-1]:
+ (Job_ID, Job_Name, User, Wall_Time, Status, Queue) = line.split()
+ jobs[Job_ID] = dict()
+ jobs[Job_ID]['name'] = Job_Name
+ jobs[Job_ID]['walltime'] = Wall_Time
+ jobs[Job_ID]['status'] = Status
+
+ elif cluster_spec == 'slurm':
+
+ end_time = datetime.now().strftime("%Y-%m-%dT%H:%M:%S")
+ cmd = ['sacct', '-n', '-X', '--parsable2', '--noheader',
+ '-S', start_time, '-E', end_time, '--format=JobName,State']
+ jobs = run(cmd, capture_output=True, text=True)
return jobs
-[docs]def check_job_status(jobs):
+
+
+[docs]
+def check_job_status(jobs, cluster_spec, start_time=None):
"""
A utility function to count the tasks with different status.
:param jobs: List of job ids.
- :return: returns the number of taks athat are queued, running, completed etc
+ :param cluster_spec: type of cluster, either 'torque' or 'slurm'.
+ :return returns the number of taks athat are queued, running, completed etc
"""
- running_jobs = retrieve_jobs()
+ running_jobs = retrieve_jobs(cluster_spec, start_time)
r = 0
c = 0
q = 0
u = 0
- for job in jobs:
- try:
- if running_jobs[job]['status'] == 'C':
+
+ if cluster_spec == 'torque':
+
+ for job in jobs:
+ try:
+ if running_jobs[job]['status'] == 'C':
+ c += 1
+ elif running_jobs[job]['status'] == 'Q':
+ q += 1
+ elif running_jobs[job]['status'] == 'R':
+ r += 1
+ else:
+ u += 1
+ except: # probably meanwhile the job is finished.
c += 1
- elif running_jobs[job]['status'] == 'Q':
- q += 1
- elif running_jobs[job]['status'] == 'R':
- r += 1
- else:
- u += 1
- except: # probably meanwhile the job is finished.
- c += 1
- continue
+ continue
+
+ print('Total Jobs:%d, Queued:%d, Running:%d, Completed:%d, Unknown:%d'
+ % (len(jobs), q, r, c, u))
+
+ elif cluster_spec == 'slurm':
+
+ lines = running_jobs.stdout.strip().split('\n')
+
+ for line in lines:
+ if line:
+ parts = line.split('|')
+ if len(parts) >= 2:
+ job_name, state = parts[0], parts[1]
+ if state == 'PENDING':
+ q += 1
+ elif state == 'RUNNING':
+ r += 1
+ elif state == 'COMPLETED':
+ c += 1
+ elif state == 'FAILED':
+ u += 1
+
+ print('Total Jobs:%d, Pending:%d, Running:%d, Completed:%d, Failed:%d'
+ % (len(jobs), q, r, c, u))
- print('Total Jobs:%d, Queued:%d, Running:%d, Completed:%d, Unknown:%d'
- % (len(jobs), q, r, c, u))
return q, r, c, u
-[docs]def check_jobs(jobs, delay=60):
+
+
+[docs]
+def check_jobs(jobs, cluster_spec, start_time=None, delay=60):
"""
A utility function for chacking the status of submitted jobs.
:param jobs: list of job ids.
+ :param cluster_spec: type of cluster, either 'torque' or 'slurm'.
:param delay: the delay (in sec) between two consequative checks, defaults to 60.
"""
@@ -1421,11 +1550,12 @@ Source code for normative_parallel
n = len(jobs)
while (True):
- q, r, c, u = check_job_status(jobs)
+ q, r, c, u = check_job_status(jobs, cluster_spec, start_time)
if c == n:
print('All jobs are completed!')
break
time.sleep(delay)
+
diff --git a/doc/build/html/_modules/rfa.html b/doc/build/html/_modules/rfa.html
index 986a44c9..0ea213e2 100644
--- a/doc/build/html/_modules/rfa.html
+++ b/doc/build/html/_modules/rfa.html
@@ -1,28 +1,24 @@
+
+
-
+
rfa — Predictive Clinical Neuroscience Toolkit 0.20 documentation
-
-
-
+
+
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
@@ -62,13 +58,24 @@
Tutorials
-
Other Useful Stuff
@@ -110,7 +117,9 @@ Source code for rfa
import torch
-[docs]class GPRRFA:
+
+[docs]
+class GPRRFA:
"""Random Feature Approximation for Gaussian Process Regression
Estimation and prediction of Bayesian linear regression models
@@ -191,11 +200,16 @@ Source code for rfa
return X, y, hyp
-[docs] def get_n_params(self, X):
+
+[docs]
+ def get_n_params(self, X):
return X.shape[1] + 2
-[docs] def post(self, hyp, X, y):
+
+
+[docs]
+ def post(self, hyp, X, y):
""" Generic function to compute posterior distribution.
This function will save the posterior mean and precision matrix as
@@ -240,7 +254,10 @@ Source code for rfa
if hasattr(self, '_iterations'):
self._iterations += 1
-[docs] def loglik(self, hyp, X, y):
+
+
+[docs]
+ def loglik(self, hyp, X, y):
""" Function to compute compute log (marginal) likelihood """
X, y, hyp = self._numpy2torch(X, y, hyp)
@@ -275,14 +292,20 @@ Source code for rfa
self.nlZ = nlZ
return nlZ
-[docs] def dloglik(self, hyp, X, y):
+
+
+[docs]
+ def dloglik(self, hyp, X, y):
""" Function to compute derivatives """
print("derivatives not available")
return
-[docs] def estimate(self, hyp0, X, y, optimizer='lbfgs'):
+
+
+[docs]
+ def estimate(self, hyp0, X, y, optimizer='lbfgs'):
""" Function to estimate the model """
if type(hyp0) is torch.Tensor:
@@ -324,7 +347,10 @@ Source code for rfa
return self.hyp.detach().numpy()
-[docs] def predict(self, hyp, X, y, Xs):
+
+
+[docs]
+ def predict(self, hyp, X, y, Xs):
""" Function to make predictions from the model """
X, y, hyp = self._numpy2torch(X, y, hyp)
@@ -347,7 +373,9 @@ Source code for rfa
torch.sum(Phis * torch.t(torch.solve(torch.t(Phis), self.A)[0]), 1)
# return output as numpy arrays
- return ys.detach().numpy().squeeze(), s2.detach().numpy().squeeze()
+ return ys.detach().numpy().squeeze(), s2.detach().numpy().squeeze()
+
+
diff --git a/doc/build/html/_modules/trendsurf.html b/doc/build/html/_modules/trendsurf.html
index d01ca0b3..1470dd83 100644
--- a/doc/build/html/_modules/trendsurf.html
+++ b/doc/build/html/_modules/trendsurf.html
@@ -1,28 +1,24 @@
+
+
-
+
trendsurf — Predictive Clinical Neuroscience Toolkit 0.20 documentation
-
-
-
+
+
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
@@ -62,13 +58,24 @@
Tutorials
-
Other Useful Stuff
@@ -135,7 +142,9 @@ Source code for trendsurf
from model.bayesreg import BLR
-[docs]def load_data(datafile, maskfile=None):
+
+[docs]
+def load_data(datafile, maskfile=None):
"""
Load data from disk
@@ -179,7 +188,10 @@ Source code for trendsurf
return dat, world, mask
-[docs]def create_basis(X, basis, mask):
+
+
+[docs]
+def create_basis(X, basis, mask):
"""
Create a basis set
@@ -223,7 +235,10 @@ Source code for trendsurf
return Phi
-[docs]def write_nii(data, filename, examplenii, mask):
+
+
+[docs]
+def write_nii(data, filename, examplenii, mask):
"""
Write data to nifti file
@@ -251,7 +266,10 @@ Source code for trendsurf
nib.save(array_img, filename)
-[docs]def get_args(*args):
+
+
+[docs]
+def get_args(*args):
"""
Parse command line arguments
@@ -293,7 +311,10 @@ Source code for trendsurf
return filename, maskfile, basis, args.a, args.o
-[docs]def estimate(filename, maskfile, basis, ard=False, outputall=False,
+
+
+[docs]
+def estimate(filename, maskfile, basis, ard=False, outputall=False,
saveoutput=True, **kwargs):
""" Estimate a trend surface model
@@ -406,13 +427,17 @@ Source code for trendsurf
return out
-[docs]def main(*args):
+
+
+[docs]
+def main(*args):
np.seterr(invalid='ignore')
filename, maskfile, basis, ard, outputall = get_args(args)
estimate(filename, maskfile, basis, ard, outputall)
+
# For running from the command line:
if __name__ == "__main__":
main(sys.argv[1:])
diff --git a/doc/build/html/_sources/pages/BLR_normativemodel_protocol.rst.txt b/doc/build/html/_sources/pages/BLR_normativemodel_protocol.rst.txt
index 479b5163..4c23523b 100644
--- a/doc/build/html/_sources/pages/BLR_normativemodel_protocol.rst.txt
+++ b/doc/build/html/_sources/pages/BLR_normativemodel_protocol.rst.txt
@@ -1,27 +1,22 @@
-.. title:: BLR tutorial
+`Predictive Clinical Neuroscience Toolkit `__
+======================================================================================
-Bayesian Linear Regression
+The Normative Modeling Framework for Computational Psychiatry Protocol
======================================================================
-The Normative Modeling Framework for Computational Psychiatry. Nature Protocols. https://www.nature.com/articles/s41596-022-00696-5.
+Using Bayesian Linear Regression and Multi-Site Cortical Thickness Data
+-----------------------------------------------------------------------
Created by `Saige Rutherford `__
-Using Multi-Site Cortical Thickness Data
-
-.. image:: https://colab.research.google.com/assets/colab-badge.svg
- :target: https://colab.research.google.com/github/predictive-clinical-neuroscience/PCNtoolkit-demo/blob/main/tutorials/BLR_protocol/BLR_normativemodel_protocol.ipynb
-
-
-.. figure:: ./blr_fig2.png
- :height: 400px
- :align: center
-
Data Preparation
----------------------------------------------
+================
Install necessary libraries & grab data files
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+---------------------------------------------
+
+Step 1.
+~~~~~~~
Begin by cloning the GitHub repository using the following commands.
This repository contains the necessary code and example data. Then
@@ -33,33 +28,28 @@ your computer).
! git clone https://github.com/predictive-clinical-neuroscience/PCNtoolkit-demo.git
-
-.. parsed-literal::
-
- Cloning into 'PCNtoolkit-demo'...
- remote: Enumerating objects: 855, done.[K
- remote: Counting objects: 100% (855/855), done.[K
- remote: Compressing objects: 100% (737/737), done.[K
- remote: Total 855 (delta 278), reused 601 (delta 101), pack-reused 0[K
- Receiving objects: 100% (855/855), 18.07 MiB | 13.53 MiB/s, done.
- Resolving deltas: 100% (278/278), done.
-
-
.. code:: ipython3
import os
# set this path to the git cloned PCNtoolkit-demo repository --> Uncomment whichever line you need for either running on your own computer or on Google Colab.
- #os.chdir('/Users/PCNtoolkit-demo/tutorials/BLR_protocol') # if running on your own computer, use this line (change the path to match where you cloned the repository)
- os.chdir('/content/PCNtoolkit-demo/tutorials/BLR_protocol') # if running on Google Colab, use this line
+ #wdir = '/PCNtoolkit-demo' # if running on your own computer, use this line (change the path to match where you cloned the repository)
+ wdir ='/content/PCNtoolkit-demo' # if running on Google Colab, use this line
+
+ os.chdir(os.path.join(wdir,'tutorials','BLR_protocol'))
+
.. code:: ipython3
+ ! pip install nutpie
+ ! pip install pcntoolkit
! pip install -r requirements.txt
-
Prepare covariate data
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+----------------------
+
+Step 2.
+~~~~~~~
The data set (downloaded in Step 1) includes a multi-site dataset from
the `Human Connectome Project Young Adult
@@ -88,8 +78,8 @@ depending on the research question.
.. code:: ipython3
# if running in Google colab, remove the "data/" folder from the path
- hcp = pd.read_csv('/content/PCNtoolkit-demo/data/HCP1200_age_gender.csv')
- ixi = pd.read_csv('/content/PCNtoolkit-demo/data/IXI_age_gender.csv')
+ hcp = pd.read_csv(os.path.join(wdir,'data','HCP1200_age_gender.csv'))
+ ixi = pd.read_csv(os.path.join(wdir,'data','IXI_age_gender.csv'))
.. code:: ipython3
@@ -98,8 +88,8 @@ depending on the research question.
.. parsed-literal::
- /usr/local/lib/python3.7/dist-packages/pandas/core/reshape/merge.py:1218: UserWarning: You are merging on int and float columns where the float values are not equal to their int representation
- UserWarning,
+ :1: UserWarning: You are merging on int and float columns where the float values are not equal to their int representation.
+ cov = pd.merge(hcp, ixi, on=["participant_id", "age", "sex", "site"], how='outer')
.. code:: ipython3
@@ -108,29 +98,348 @@ depending on the research question.
.. code:: ipython3
- sns.displot(cov, x="age", hue="site", multiple="stack", height=6)
+ sns.displot(cov, x="age", hue="site", multiple="stack", height=6);
+.. image:: BLR_normativemodel_protocol_files/BLR_normativemodel_protocol_16_0.png
-.. parsed-literal::
-
+.. code:: ipython3
+ cov.groupby(['site']).describe()
-.. image:: BLR_normativemodel_protocol_files/BLR_normativemodel_protocol_15_1.png
+.. raw:: html
+
+
+
+
+
+
+
+
+ age
+ sex
+
+
+ count
+ mean
+ std
+ min
+ 25%
+ 50%
+ 75%
+ max
+ count
+ mean
+ std
+ min
+ 25%
+ 50%
+ 75%
+ max
+
+
+ site
+
+
+
+
+ hcp
+ 1206.0
+ 28.837479
+ 3.690534
+ 22.000000
+ 26.000000
+ 29.00000
+ 32.000000
+ 37.00000
+ 1206.0
+ 1.543947
+ 0.498272
+ 1.0
+ 1.0
+ 2.0
+ 2.0
+ 2.0
+
+
+ ixi
+ 590.0
+ 49.476531
+ 16.720864
+ 19.980835
+ 34.027721
+ 50.61191
+ 63.413415
+ 86.31896
+ 590.0
+ 1.555932
+ 0.497283
+ 1.0
+ 1.0
+ 2.0
+ 2.0
+ 2.0
+
+
+
+
+
+
-.. code:: ipython3
- cov.groupby(['site']).describe()
+Preprare brain data
+-------------------
-Prepare brain data
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Step 3.
+~~~~~~~
Next, format and combine the MRI data using the following commands. The
example data contains cortical thickness maps estimated by running
@@ -139,7 +448,7 @@ was reduced by using ROIs from the Desikan-Killiany atlas. Including the
Euler number as a covariate is also recommended, as this is a proxy
metric for data quality. The `Euler
number `__ from
-each subjects recon-all output folder was extracted into a text file
+each subject’s recon-all output folder was extracted into a text file
and is merged into the cortical thickness data frame. The Euler number
is site-specific, thus, to use the same exclusion threshold across sites
it is important to center the site by subtracting the site median from
@@ -147,7 +456,7 @@ all subjects at a site. Then take the square root and multiply by
negative one and exclude any subjects with a square root above 10.
Here is some psuedo-code (run from a terminal in the folder that has all
-subjects recon-all output folders) that was used to extract these ROIs:
+subject’s recon-all output folders) that was used to extract these ROIs:
``export SUBJECTS_DIR=/path/to/study/freesurfer_data/``
@@ -157,14 +466,14 @@ subjects recon-all output folders) that was used to extract these ROIs:
.. code:: ipython3
- hcpya = pd.read_csv('/content/PCNtoolkit-demo/data/HCP1200_aparc_thickness.csv')
- ixi = pd.read_csv('/content/PCNtoolkit-demo/data/IXI_aparc_thickness.csv')
+ hcpya = pd.read_csv(os.path.join(wdir,'data','HCP1200_aparc_thickness.csv'))
+ ixi = pd.read_csv(os.path.join(wdir,'data','IXI_aparc_thickness.csv'))
.. code:: ipython3
brain_all = pd.merge(ixi, hcpya, how='outer')
-We extracted the euler number from each subjects recon-all output
+We extracted the euler number from each subject’s recon-all output
folder into a text file and we now need to format and combine these into
our brain dataframe.
@@ -173,12 +482,12 @@ recon-all.log for each subject. Run this from the terminal in the folder
where your subjects recon-all output folders are located. This assumes
that all of your subject IDs start with “sub-” prefix.
-:literal:`for i in sub-*; do if [[ -e ${i}/scripts/recon-all.log ]]; then cat ${i}/scripts/recon-all.log | grep -A 1 "Computing euler" > temp_log; lh_en=$(cat temp_log | head -2 | tail -1 | awk -F '=' '{print $2}' | awk -F ',' '{print $1}'); rh_en=$(cat temp_log | head -2 | tail -1 | awk -F '=' '{print $3}'); echo "${i}, ${lh_en}, ${rh_en}" >> euler.csv; echo ${i}; fi; done`
+:literal:`for i in sub-\*; do if [[ -e ${i}/scripts/recon-all.log ]]; then cat ${i}/scripts/recon-all.log | grep -A 1 "Computing euler" > temp_log; lh_en=`cat temp_log | head -2 | tail -1 | awk -F '=' '{print $2}' | awk -F ',' '{print $1}'\`; rh_en=`cat temp_log | head -2 | tail -1 | awk -F '=' '{print $3}'\`; echo "${i}, ${lh_en}, ${rh_en}" >> euler.csv; echo ${i}; fi; done`
.. code:: ipython3
- hcp_euler = pd.read_csv('/content/PCNtoolkit-demo/data/hcp-ya_euler.csv')
- ixi_euler = pd.read_csv('/content/PCNtoolkit-demo/data/ixi_euler.csv')
+ hcp_euler = pd.read_csv(os.path.join(wdir,'data','hcp-ya_euler.csv'))
+ ixi_euler = pd.read_csv(os.path.join(wdir,'data','ixi_euler.csv'))
.. code:: ipython3
@@ -221,7 +530,7 @@ inclusion is not too strict or too lenient.
.. code:: ipython3
- df_euler.groupby(by='site').median()
+ df_euler.groupby(by='site')[['lh_euler', 'rh_euler', 'avg_euler']].median()
@@ -229,9 +538,8 @@ inclusion is not too strict or too lenient.
.. raw:: html
-
-
-
+
+
-
+
+
+
+
+
+