Skip to content

[RTM] RF: Split recon-all into coarse chunks to improve resource use estimates #506

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 13 commits into from
May 17, 2017
Merged

[RTM] RF: Split recon-all into coarse chunks to improve resource use estimates #506

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
2a4a013
RF: Split recon-all directives to adjust parallelism
effigies May 2, 2017
b1c9588
De-duplicate surface outputs
effigies May 3, 2017
cd9a0e6
Remove get_surfaces iterable node
effigies May 3, 2017
47e15ff
Use autorecon3 instead of recon-all -all to finish
effigies May 5, 2017
c9b0c8a
Update recon node names
effigies May 5, 2017
4172e59
DOC: Update changelog
effigies May 10, 2017
cc7e105
Adjust recon flags
effigies May 11, 2017
1833ab3
RF: Use iterables instead of MapNode to manage parallelism
effigies May 14, 2017
7c8772e
RF: Separate autorecon_resume and gifti_surface workflows
effigies May 14, 2017
2b60563
Add sync point, postpone hyporelabel
effigies May 14, 2017
ba3367e
RF: Switch back to MapNodes for conciseness
effigies May 16, 2017
01de20a
DOC: Describe subprocess breakdown
effigies May 16, 2017
a52a040
DOC: Update doc pins
effigies May 17, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGES.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ Next release
============

* [ENH] Improved EPI skull stripping and tissue contrast enhancements (#519)
* [ENH] Improve resource use estimates in FreeSurfer workflow (#506)

0.4.3 (10th of May 2017)
========================
Expand Down
4 changes: 2 additions & 2 deletions docs/environment.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,6 @@ dependencies:
- svgutils
- nitime
- nilearn
- git+https://github.com/nipy/nipype.git@d82a18f80126b4bc99d0b4f6d14ff3b6e075f297#egg=nipype
- git+https://github.com/poldracklab/niworkflows.git@c309563cb5763464b607c45a20ffde92293194a1#egg=niworkflows
- git+https://github.com/nipy/nipype.git@4e48d2644b19e74581904e5404569e505279835a#egg=nipype
- git+https://github.com/poldracklab/niworkflows.git@1f6ee8fe28f85b30398e33424c8b4d0b2f20992f#egg=niworkflows
- git+https://github.com/INCF/pybids.git@7205ae01fbdca8e8bfd26ac773b1134b84c8af0c#egg=pybids
21 changes: 18 additions & 3 deletions docs/workflows.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ Surface preprocessing
:mod:`fmriprep.workflows.anatomical.init_surface_recon_wf`

.. workflow::
:graph2use: colored
:graph2use: orig
:simple_form: yes

from fmriprep.workflows.anatomical import init_surface_recon_wf
Expand Down Expand Up @@ -112,10 +112,25 @@ would be processed by the following command::
The second phase imports the brainmask calculated in the `T1w/T2w preprocessing`_
sub-workflow.
The final phase resumes reconstruction, using the T2w image to assist
in finding the pial surface, if available::
in finding the pial surface, if available.
In order to utilize resources efficiently, this is broken down into six
sub-stages; the first and last are run serially, while each pair of
per-hemisphere stages are run in parallel, if possible::

$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-all -T2pial
-autorecon2-volonly
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon2-perhemi -hemi lh
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon2-perhemi -hemi rh
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon-hemi lh -T2pial \
-noparcstats -noparcstats2 -noparcstats3 -nohyporelabel -nobalabels
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon-hemi rh -T2pial \
-noparcstats -noparcstats2 -noparcstats3 -nohyporelabel -nobalabels
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon3

Reconstructed white and pial surfaces are included in the report.

Expand Down
180 changes: 131 additions & 49 deletions fmriprep/workflows/anatomical.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -307,8 +307,7 @@ def bidsinfo(in_file):
fs.ReconAll(
directive='autorecon1',
flags='-noskullstrip',
openmp=omp_nthreads,
parallel=True),
openmp=omp_nthreads),
name='autorecon1')
autorecon1.interface._can_resume = False
autorecon1.interface.num_threads = omp_nthreads
Expand Down Expand Up @@ -341,32 +340,142 @@ def inject_skullstripped(subjects_dir, subject_id, skullstripped):
output_names=['subjects_dir', 'subject_id']),
name='skull_strip_extern')

reconall = pe.Node(
fs_transform = pe.Node(
fs.Tkregister2(fsl_out='freesurfer2subT1.mat', reg_header=True),
name='fs_transform')

autorecon_resume_wf = init_autorecon_resume_wf(omp_nthreads=omp_nthreads)
gifti_surface_wf = init_gifti_surface_wf()

workflow.connect([
# Configuration
(inputnode, recon_config, [('t1w', 't1w_list'),
('t2w', 't2w_list')]),
(inputnode, bids_info, [(('t1w', fix_multi_T1w_source_name), 'in_file')]),
# Passing subjects_dir / subject_id enforces serial order
(inputnode, autorecon1, [('subjects_dir', 'subjects_dir')]),
(bids_info, autorecon1, [('subject_id', 'subject_id')]),
(autorecon1, skull_strip_extern, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(skull_strip_extern, autorecon_resume_wf, [('subjects_dir', 'inputnode.subjects_dir'),
('subject_id', 'inputnode.subject_id')]),
(autorecon_resume_wf, gifti_surface_wf, [
('outputnode.subjects_dir', 'inputnode.subjects_dir'),
('outputnode.subject_id', 'inputnode.subject_id')]),
# Reconstruction phases
(inputnode, autorecon1, [('t1w', 'T1_files')]),
(recon_config, autorecon1, [('t2w', 'T2_file'),
('hires', 'hires'),
# First run only (recon-all saves expert options)
('mris_inflate', 'mris_inflate')]),
(inputnode, skull_strip_extern, [('skullstripped_t1', 'skullstripped')]),
(recon_config, autorecon_resume_wf, [('use_T2', 'inputnode.use_T2')]),
# Construct transform from FreeSurfer conformed image to FMRIPREP
# reoriented image
(inputnode, fs_transform, [('t1w', 'target_image')]),
(autorecon1, fs_transform, [('T1', 'moving_image')]),
# Output
(autorecon_resume_wf, outputnode, [('outputnode.subjects_dir', 'subjects_dir'),
('outputnode.subject_id', 'subject_id'),
('outputnode.out_report', 'out_report')]),
(gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces')]),
(fs_transform, outputnode, [('fsl_file', 'fs_2_t1_transform')]),
])

return workflow


def init_autorecon_resume_wf(omp_nthreads, name='autorecon_resume_wf'):
workflow = pe.Workflow(name=name)

inputnode = pe.Node(
niu.IdentityInterface(
fields=['subjects_dir', 'subject_id', 'use_T2']),
name='inputnode')

outputnode = pe.Node(
niu.IdentityInterface(
fields=['subjects_dir', 'subject_id', 'out_report']),
name='outputnode')

autorecon2_vol = pe.Node(
fs.ReconAll(
directive='autorecon2-volonly',
openmp=omp_nthreads),
name='autorecon2_vol')
autorecon2_vol.interface.num_threads = omp_nthreads

autorecon2_surfs = pe.MapNode(
fs.ReconAll(
directive='autorecon2-perhemi',
openmp=omp_nthreads),
iterfield='hemi',
name='autorecon2_surfs')
autorecon2_surfs.interface.num_threads = omp_nthreads
autorecon2_surfs.inputs.hemi = ['lh', 'rh']

autorecon_surfs = pe.MapNode(
fs.ReconAll(
directive='autorecon-hemi',
flags=['-noparcstats', '-noparcstats2', '-noparcstats3',
'-nohyporelabel', '-nobalabels'],
openmp=omp_nthreads),
iterfield='hemi',
name='autorecon_surfs')
autorecon_surfs.interface.num_threads = omp_nthreads
autorecon_surfs.inputs.hemi = ['lh', 'rh']
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On a 6-core machine, I'm getting both 5-thread jobs in this mapnode running simultaneously. @oesteban I think you're more familiar with the resource manager than I am. Do you see something wrong here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

doesn't omp_nthreahds default to nthreads-1?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. But we shouldn't be able to run two 5-thread jobs on a 6-thread queue... I think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes that is correct. I wonder if this is because MultiProc is not handling MapNodes correctly.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I think the possibilities are:

  1. MapNodes are handled incorrectly by MultiProc
  2. MapNodes need to have num_threads specified differently from other nodes
  3. I've screwed something subtle up here

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be worth opening an issue at nipy/nipype and pinging Cameron.


autorecon3 = pe.Node(
ReconAllRPT(
flags='-noskullstrip',
directive='autorecon3',
openmp=omp_nthreads,
parallel=True,
out_report='reconall.svg',
generate_report=True),
name='reconall')
reconall.interface.num_threads = omp_nthreads
name='autorecon3')
autorecon3.interface.num_threads = omp_nthreads

fs_transform = pe.Node(
fs.Tkregister2(fsl_out='freesurfer2subT1.mat', reg_header=True),
name='fs_transform')
def _dedup(in_list):
vals = set(in_list)
if len(vals) > 1:
raise ValueError(
"Non-identical values can't be deduplicated:\n{!r}".format(in_list))
return vals.pop()

get_surfaces = pe.Node(nio.FreeSurferSource(), iterables=('hemi', ('lh', 'rh')),
name='get_surfaces')
workflow.connect([
(inputnode, autorecon_surfs, [('use_T2', 'use_T2')]),
(inputnode, autorecon2_vol, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(autorecon2_vol, autorecon2_surfs, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(autorecon2_surfs, autorecon_surfs, [(('subjects_dir', _dedup), 'subjects_dir'),
(('subject_id', _dedup), 'subject_id')]),
(autorecon_surfs, autorecon3, [(('subjects_dir', _dedup), 'subjects_dir'),
(('subject_id', _dedup), 'subject_id')]),
(autorecon3, outputnode, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id'),
('out_report', 'out_report')]),
])

return workflow


def init_gifti_surface_wf(name='gifti_surface_wf'):
workflow = pe.Workflow(name=name)

inputnode = pe.Node(niu.IdentityInterface(['subjects_dir', 'subject_id']), name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(['surfaces']), name='outputnode')

midthickness = pe.Node(MakeMidthickness(thickness=True, distance=0.5, out_name='midthickness'),
name='midthickness')
get_surfaces = pe.Node(nio.FreeSurferSource(), name='get_surfaces')

midthickness = pe.MapNode(
MakeMidthickness(thickness=True, distance=0.5, out_name='midthickness'),
iterfield='in_file',
name='midthickness')

save_midthickness = pe.Node(nio.DataSink(parameterization=False),
name='save_midthickness')

surface_list = pe.JoinNode(niu.Merge(4, ravel_inputs=True), name='surface_list',
joinsource='get_surfaces', joinfield=['in1', 'in2', 'in3', 'in4'],
run_without_submitting=True)
surface_list = pe.Node(niu.Merge(4, ravel_inputs=True),
name='surface_list', run_without_submitting=True)
fs_2_gii = pe.MapNode(fs.MRIsConvert(out_datatype='gii'),
iterfield='in_file', name='fs_2_gii')

Expand Down Expand Up @@ -421,37 +530,10 @@ def normalize_surfs(in_file):
name='fix_surfs')

workflow.connect([
# Configuration
(inputnode, recon_config, [('t1w', 't1w_list'),
('t2w', 't2w_list')]),
(inputnode, bids_info, [(('t1w', fix_multi_T1w_source_name), 'in_file')]),
# Passing subjects_dir / subject_id enforces serial order
(inputnode, autorecon1, [('subjects_dir', 'subjects_dir')]),
(bids_info, autorecon1, [('subject_id', 'subject_id')]),
(autorecon1, skull_strip_extern, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(skull_strip_extern, reconall, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(reconall, get_surfaces, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(reconall, save_midthickness, [('subjects_dir', 'base_directory'),
('subject_id', 'container')]),
(reconall, outputnode, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id'),
('out_report', 'out_report')]),
# Reconstruction phases
(inputnode, autorecon1, [('t1w', 'T1_files')]),
(recon_config, autorecon1, [('t2w', 'T2_file'),
('hires', 'hires'),
# First run only (recon-all saves expert options)
('mris_inflate', 'mris_inflate')]),
(inputnode, skull_strip_extern, [('skullstripped_t1', 'skullstripped')]),
(recon_config, reconall, [('use_T2', 'use_T2')]),
# Construct transform from FreeSurfer conformed image to FMRIPREP
# reoriented image
(inputnode, fs_transform, [('t1w', 'target_image')]),
(autorecon1, fs_transform, [('T1', 'moving_image')]),
(fs_transform, outputnode, [('fsl_file', 'fs_2_t1_transform')]),
(inputnode, get_surfaces, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(inputnode, save_midthickness, [('subjects_dir', 'base_directory'),
('subject_id', 'container')]),
# Generate midthickness surfaces and save to FreeSurfer derivatives
(get_surfaces, midthickness, [('smoothwm', 'in_file'),
('graymid', 'graymid')]),
Expand Down