Skip to content

Random permutation of FCMA #217

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
xoltar merged 10 commits into from
Apr 28, 2017
Merged

Random permutation of FCMA #217

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
81972f7
add random permutation to FCMA
yidawang Apr 27, 2017
507bbf3
addressing review comments from Bryn
yidawang Apr 27, 2017
5f04e95
addressing review comments from Mihai
yidawang Apr 27, 2017
cec3445
further changes according to the comments
yidawang Apr 27, 2017
ec340ec
addressing review comments
yidawang Apr 27, 2017
0e28424
Merge remote-tracking branch 'upstream/master' into fcma-random-permu…
yidawang Apr 27, 2017
e90671a
Merge remote-tracking branch 'upstream/master' into fcma-random-permu…
yidawang Apr 28, 2017
11383d1
addressing Bryn's comments
yidawang Apr 28, 2017
ec78ccd
updating docstrings
yidawang Apr 28, 2017
e8d60a8
fix building error issue
yidawang Apr 28, 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
108 changes: 98 additions & 10 deletions brainiak/fcma/preprocessing.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,21 @@
import logging
from scipy.stats.mstats import zscore
from mpi4py import MPI
from enum import Enum

from ..image import mask_images, multimask_images


logger = logging.getLogger(__name__)

__all__ = [
"RandomType",
"prepare_fcma_data",
"generate_epochs_info",
"prepare_mvpa_data",
"prepare_searchlight_mvpa_data",
]


def _separate_epochs(activity_data, epoch_list):
""" create data epoch by epoch
Expand All @@ -37,9 +46,9 @@ def _separate_epochs(activity_data, epoch_list):

Parameters
----------
activity\_data: list of 2D array in shape [nVoxels, nTRs]
activity_data: list of 2D array in shape [nVoxels, nTRs]
the masked activity data organized in voxel*TR formats of all subjects
epoch\_list: list of 3D array in shape [condition, nEpochs, nTRs]
epoch_list: list of 3D array in shape [condition, nEpochs, nTRs]
specification of epochs and conditions
assuming all subjects have the same number of epochs
len(epoch_list) equals the number of subjects
Expand Down Expand Up @@ -83,8 +92,69 @@ def _separate_epochs(activity_data, epoch_list):
return raw_data, labels


def _randomize_single_subject(data, seed=None):
"""Randomly permute the voxels of the subject.

The subject is organized as Voxel x TR,
this method shuffles the voxel dimension in place.

Parameters
----------
data: 2D array in shape [nVoxels, nTRs]
Activity data to be shuffled.
seed: Optional[int]
Seed for random state used implicitly for shuffling.

Returns
-------
None.
"""
if seed is not None:
np.random.seed(seed)
np.random.shuffle(data)


def _randomize_subject_list(data_list, random):
"""Randomly permute the voxels of a subject list.

The method shuffles the subject one by one in place according to
the random type. If RandomType.NORANDOM, return the original list.

Parameters
----------
data_list: list of 2D array in shape [nVxels, nTRs]
Activity data list to be shuffled.
random: RandomType
Randomization type.

Returns
-------
None.
"""
if random == RandomType.REPRODUCIBLE:
for i in range(len(data_list)):
_randomize_single_subject(data_list[i], seed=i)
elif random == RandomType.UNREPRODUCIBLE:
for data in data_list:
_randomize_single_subject(data)


class RandomType(Enum):
"""Define the random types as enumeration

NORANDOM means do not randomize the data;
REPRODUCIBLE means randomize the data with a fixed seed so that the
permutation holds between different runs;
UNREPRODUCIBLE means truly randomize the data which returns different
results in different runs.
"""
NORANDOM = 0
REPRODUCIBLE = 1
UNREPRODUCIBLE = 2


def prepare_fcma_data(images, conditions, mask1, mask2=None,
comm=MPI.COMM_WORLD):
random=RandomType.NORANDOM, comm=MPI.COMM_WORLD):
"""Prepare data for correlation-based computation and analysis.

Generate epochs of interests, then broadcast to all workers.
Expand All @@ -102,6 +172,9 @@ def prepare_fcma_data(images, conditions, mask1, mask2=None,
If it is not specified, the method will assign None to the returning
variable raw_data2 and the self-correlation on raw_data1 will be
computed
random: Optional[RandomType]
Randomize the data within subject or not.
Default NORANDOM
comm: MPI.Comm
MPI communicator to use for MPI operations.

Expand All @@ -122,14 +195,17 @@ def prepare_fcma_data(images, conditions, mask1, mask2=None,
raw_data1 = []
raw_data2 = []
if rank == 0:
logger.info('start to apply masks and separate epochs')
if mask2 is not None:
masks = (mask1, mask2)
activity_data1, activity_data2 = zip(*multimask_images(images,
masks,
np.float32))
_randomize_subject_list(activity_data2, random)
raw_data2, _ = _separate_epochs(activity_data2, conditions)
else:
activity_data1 = list(mask_images(images, mask1, np.float32))
_randomize_subject_list(activity_data1, random)

@xoltar xoltar Apr 27, 2017

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

_randomize_single_subject modifies data in place (and then returns the same array). _randomize_subject_list returns a new list, containing old arrays that have been modified in place, which is why these lines work even though you're ignoring the return value. I would say either write them both to work in-place, or write them both to return new values. Mixing the two makes for confusion. If you decide you want to have _randomize_single_subject return a new array instead of working in-place, consider numpy.random.permutation instead of shuffle.

@yidawang yidawang Apr 27, 2017

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

Thanks for the suggestion! Somehow I ignored the return value here although I should put it as activity_data1 = _randomize_subject_list(activity_data1, random). I would prefer working in-place, let me modify _randomize_subject_list accordingly.

raw_data1, labels = _separate_epochs(activity_data1, conditions)
time1 = time.time()
raw_data_length = len(raw_data1)
Expand Down Expand Up @@ -170,7 +246,7 @@ def generate_epochs_info(epoch_list):

Returns
-------
epoch\_info: list of tuple (label, sid, start, end).
epoch_info: list of tuple (label, sid, start, end).
label is the condition labels of the epochs;
sid is the subject id, corresponding to the index of raw_data;
start is the start TR of an epoch (inclusive);
Expand Down Expand Up @@ -250,7 +326,8 @@ def prepare_mvpa_data(images, conditions, mask):
return processed_data, labels


def prepare_searchlight_mvpa_data(images, conditions, data_type=np.float32):
def prepare_searchlight_mvpa_data(images, conditions, data_type=np.float32,
random=RandomType.NORANDOM):
""" obtain the data for activity-based voxel selection using Searchlight

Average the activity within epochs and z-scoring within subject,
Expand All @@ -266,10 +343,12 @@ def prepare_searchlight_mvpa_data(images, conditions, data_type=np.float32):
Condition specification.
data_type
Type to cast image to.
random: Optional[RandomType]
Randomize the data within subject or not.

Returns
-------
processed\_data: 4D array in shape [brain 3D + epoch]
processed_data: 4D array in shape [brain 3D + epoch]
averaged epoch by epoch processed data

labels: 1D array
Expand All @@ -290,10 +369,19 @@ def prepare_searchlight_mvpa_data(images, conditions, data_type=np.float32):
# counting the epochs per subject for z-scoring
subject_count = np.zeros(len(conditions), dtype=np.int32)

logger.info('start to apply masks and separate epochs')
for sid, f in enumerate(images):
data = f.get_data().astype(data_type)
[d1, d2, d3, d4] = data.shape
if random == RandomType.REPRODUCIBLE:
data = data.reshape((d1 * d2 * d3, d4))
_randomize_single_subject(data, seed=sid)
data = data.reshape((d1, d2, d3, d4))
elif random == RandomType.UNREPRODUCIBLE:
data = data.reshape((d1 * d2 * d3, d4))
_randomize_single_subject(data)
data = data.reshape((d1, d2, d3, d4))
if processed_data is None:
[d1, d2, d3, _] = data.shape
processed_data = np.empty([d1, d2, d3, num_epochs],
dtype=data_type)
# averaging
Expand All @@ -303,9 +391,9 @@ def prepare_searchlight_mvpa_data(images, conditions, data_type=np.float32):
processed_data[:, :, :, idx] = \
np.mean(data[:, :, :, epoch[2]:epoch[3]], axis=3)

logger.info(
'file %s is loaded and processed, with data shape %s' %
(f.get_filename(), data.shape)
logger.debug(
'file %s is loaded and processed, with data shape %s',
f.get_filename(), data.shape
)
# z-scoring
cur_epoch = 0
Expand Down
16 changes: 15 additions & 1 deletion brainiak/fcma/voxelselector.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,8 @@ def _master(self):
the length of array equals the number of voxels
"""
logger.info(
'Master starts to allocate tasks'
'Master at rank %d starts to allocate tasks',
MPI.COMM_WORLD.Get_rank()
)
results = []
comm = MPI.COMM_WORLD
Expand All @@ -181,6 +182,10 @@ def _master(self):
if current_task[1] == 0:
using_size = i
break
logger.debug(
'master starts to send a task to worker %d' %
i
)
comm.send(current_task,
dest=i,
tag=self._WORKTAG)
Expand Down Expand Up @@ -235,6 +240,10 @@ def _worker(self, clf):
-------
None
"""
logger.debug(
'worker %d is running, waiting for tasks from master at rank %d' %
(MPI.COMM_WORLD.Get_rank(), self.master_rank)
)
comm = MPI.COMM_WORLD
status = MPI.Status()
while 1:
Expand Down Expand Up @@ -264,6 +273,11 @@ def _correlation_computation(self, task):
time1 = time.time()
s = task[0]
nEpochs = len(self.raw_data)
logger.debug(
'start to compute the correlation: #epochs: %d, '
'#processed voxels: %d, #total voxels to compute against: %d' %
(nEpochs, task[1], self.num_voxels2)
)
corr = np.zeros((task[1], nEpochs, self.num_voxels2),
np.float32, order='C')
count = 0
Expand Down
11 changes: 10 additions & 1 deletion brainiak/io.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
from pathlib import Path
from typing import Callable, Iterable, List, Union

import logging
import nibabel as nib
import numpy as np

Expand All @@ -31,6 +32,8 @@

from .image import SingleConditionSpec

logger = logging.getLogger(__name__)


def load_images_from_dir(in_dir: Union[str, Path], suffix: str = "nii.gz",
) -> Iterable[SpatialImage]:
Expand Down Expand Up @@ -58,6 +61,9 @@ def load_images_from_dir(in_dir: Union[str, Path], suffix: str = "nii.gz",
in_dir = Path(in_dir)
files = sorted(in_dir.glob("*" + suffix))
for f in files:
logger.debug(
'Starting to read file %s', f
)
yield nib.load(str(f))


Expand Down Expand Up @@ -86,6 +92,9 @@ def load_images(image_paths: Iterable[Union[str, Path]]
string_path = str(image_path)
else:
string_path = image_path
logger.debug(
'Starting to read file %s', string_path
)
yield nib.load(string_path)


Expand Down Expand Up @@ -130,7 +139,7 @@ def load_labels(path: Union[str, Path]) -> List[SingleConditionSpec]:
List[SingleConditionSpec]
List of SingleConditionSpec stored in labels file.
"""
condition_specs = np.load(path)
condition_specs = np.load(str(path))
return [c.view(SingleConditionSpec) for c in condition_specs]


Expand Down
17 changes: 15 additions & 2 deletions examples/fcma/mvpa_voxel_selection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@
MPI.COMM_WORLD.Get_size()
)
data_dir = sys.argv[1]
extension = sys.argv[2]
suffix = sys.argv[2]
mask_file = sys.argv[3]
epoch_file = sys.argv[4]

Expand All @@ -53,11 +53,24 @@
'mask size: %d' %
np.sum(mask)
)
images = io.load_images_from_dir(data_dir, extension)
images = io.load_images_from_dir(data_dir, suffix=suffix)
conditions = io.load_labels(epoch_file)
data, labels = prepare_searchlight_mvpa_data(images, conditions)

# setting the random argument produces random voxel selection results
# for non-parametric statistical analysis.
# There are three random options:
# RandomType.NORANDOM is the default
# RandomType.REPRODUCIBLE permutes the voxels in the same way every run
# RandomType.UNREPRODUCIBLE permutes the voxels differently across runs
# example
#from brainiak.fcma.preprocessing import RandomType
#data, labels = prepare_searchlight_mvpa_data(images, conditions,
# random=RandomType.UNREPRODUCIBLE)

# the following line is an example to leaving a subject out
#epoch_info = [x for x in epoch_info if x[1] != 0]

num_subjs = int(sys.argv[5])
# create a Searchlight object
sl = Searchlight(sl_rad=1)
Expand Down
18 changes: 16 additions & 2 deletions examples/fcma/voxel_selection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,19 +38,33 @@
MPI.COMM_WORLD.Get_size()
)
data_dir = sys.argv[1]
extension = sys.argv[2]
suffix = sys.argv[2]
mask_file = sys.argv[3]
epoch_file = sys.argv[4]
images = io.load_images_from_dir(data_dir, extension)
images = io.load_images_from_dir(data_dir, suffix=suffix)
mask = io.load_boolean_mask(mask_file)
conditions = io.load_labels(epoch_file)
raw_data, _, labels = prepare_fcma_data(images, conditions, mask)

# setting the random argument produces random voxel selection results
# for non-parametric statistical analysis.
# There are three random options:
# RandomType.NORANDOM is the default
# RandomType.REPRODUCIBLE permutes the voxels in the same way every run
# RandomType.UNREPRODUCIBLE permutes the voxels differently across runs
# example:
# from brainiak.fcma.preprocessing import RandomType
# raw_data, _, labels = prepare_fcma_data(images, conditions, mask,
# random=RandomType.REPRODUCIBLE)

# if providing two masks, just append the second mask as the last input argument
# and specify raw_data2
# example:
# images = io.load_images_from_dir(data_dir, extension)
# mask2 = io.load_boolean_mask('face_scene/mask.nii.gz')
# raw_data, raw_data2, labels = prepare_fcma_data(images, conditions, mask,
# mask2)

epochs_per_subj = int(sys.argv[5])
num_subjs = int(sys.argv[6])
# the following line is an example to leaving a subject out
Expand Down
Loading