Skip to content
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

Gtm tortuosity #84

Closed
wants to merge 3 commits into from
Closed

Gtm tortuosity #84

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fc2d6ff
add GTM parameters in set_tortuous_parameter
matteofrigo Oct 28, 2019
e86da9a
Merge branch 'master' of github.com:AthenaEPI/dmipy into gtm-tortuosity
matteofrigo Oct 28, 2019
0748eaa
change behaviour of S0 response assignment in GTM
matteofrigo Oct 31, 2019
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
26 changes: 23 additions & 3 deletions dmipy/core/modeling_framework.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -558,7 +558,8 @@ def _add_fixed_parameter_array(self, parameter_name, parameter_array):
def set_tortuous_parameter(self, lambda_perp_parameter_name,
lambda_par_parameter_name,
volume_fraction_intra_parameter_name,
volume_fraction_extra_parameter_name):
volume_fraction_extra_parameter_name,
S0_responses=None):
Copy link
Collaborator

Choose a reason for hiding this comment

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

Again (see below) I don't think we should be able to give S0 responses twice for any MC-model - it should be unique in the definition.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Is there a shortcut for knowing which of the S0s in self.S0_tissue_responses list corresponds to the intra and which one to the extra?

Copy link
Collaborator

Choose a reason for hiding this comment

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

I didn't even thing about this when I first coded this... I especially kept the parameter names as input because I didn't want to have to hardcode the ordering or definitions of models to get the tortuosity regardless.

Actually we do need to find a solution here - what you were doing with adding the S0_responses as a parameter may be the easiest way to make this work after all... Either this or we need to define a nice way to just point to the models that are intra and extra, and then make the function just find which parameters and S0s are associated with this...

"""
Allows the user to set a tortuosity constraint on the perpendicular
diffusivity of the extra-axonal compartment, which depends on the
Expand All @@ -581,6 +582,8 @@ def set_tortuous_parameter(self, lambda_perp_parameter_name,
volume_fraction_extra_parameter_name: string
name of the extra-axonal volume fraction parameter, see
self.parameter_names.
S0_responses : list, (optional)
s0 response of the tissues associated to IC and EC.
"""
params = [lambda_perp_parameter_name, lambda_par_parameter_name,
volume_fraction_intra_parameter_name,
Expand All @@ -597,7 +600,8 @@ def set_tortuous_parameter(self, lambda_perp_parameter_name,
self.parameter_links.append([model, name, T1_tortuosity, [
self._parameter_map[lambda_par_parameter_name],
self._parameter_map[volume_fraction_intra_parameter_name],
self._parameter_map[volume_fraction_extra_parameter_name]]
self._parameter_map[volume_fraction_extra_parameter_name],
S0_responses]
])
del self.parameter_ranges[lambda_perp_parameter_name]
del self.parameter_cardinality[lambda_perp_parameter_name]
Expand Down Expand Up @@ -1212,8 +1216,24 @@ def fit(self, acquisition_scheme, data,
start = time()
mt_fractions = np.empty(
np.r_[N_voxels, self.N_models], dtype=float)

Copy link
Collaborator

Choose a reason for hiding this comment

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

You need to add this also for MC-SM models.

S0_responses = []
Copy link
Collaborator

Choose a reason for hiding this comment

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

We know what will be the length of S0_responses, so it's cleaner to initiate it just with S0_responses = np.zeros(self.N_models) or something.

for i, (m, s) in enumerate(zip(self.models,
matteofrigo marked this conversation as resolved.
Show resolved Hide resolved
self.S0_tissue_responses)):
if s is None:
# If the S0 of the tissue is None, the model is composite
# and the corresponding S0 is treated as the weighted
# average of the S0 tissues modelled by the "submodels".
s0 = 0.0
for k, mm in enumerate(m.models):
Copy link
Collaborator

Choose a reason for hiding this comment

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

Nicely done, but while this may work here at the top level to the second level - just for this project - it means that the top level now has the logic to infer what is the S0 of a lower level distributed model.

In my opinion this is not so clean, as the we should be able to ask the distributed model what is it's S0 based on its signal volume fractions. This has the added benefit of it being to (potentially) generalize to more than 2 layers - we want to have the code as modular as possible.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It makes sense.
If I understand the flow correctly, this could be done by moving the definition of S0_responses to the __init__ method, where self.S0_tissue_responses will be defined inferring the different S0s from the "submodels".
Am I right?

Copy link
Collaborator

Choose a reason for hiding this comment

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

What I mean is that we should be able to say S0 = m.S0_response(volume_fractions) to any submodel, instead of having to manually sum up the s0 values in line 1229.

Hold on, I how is line 1229 even getting the 'partial_volume' values? mm is just a DistributedModel instance right? not a dict? I don't thing this works like this, or am I wrong?

s0 += (mm['partial_volume_{}'.format(k)] *
mm.S0_tissue_responses[k])
S0_responses.append(s0)
else:
S0_responses.append(s)

fit_func = MultiTissueConvexOptimizer(
acquisition_scheme, self, self.S0_tissue_responses)
acquisition_scheme, self, S0_responses)
for idx, pos in enumerate(zip(*mask_pos)):
voxel_S = data_[pos]
parameters = fitted_parameters_lin[idx]
Expand Down
10 changes: 8 additions & 2 deletions dmipy/distributions/distribute_models.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,8 @@ def set_fixed_parameter(self, parameter_name, value):

def set_tortuous_parameter(self, lambda_perp,
lambda_par,
volume_fraction_intra):
volume_fraction_intra,
S0_responses=None):
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think the S0_responses should not be given as a parameter of the set_tortuous_parameter function, but should be given to the SD1_WatsonDistributedModel as a second argument, and before a property as self.S0_responses, just as in the upper MC-models.

You could, potentially, add a parameter here like use_S0_debiasing_if_possible or something to just be able to choose here.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

For the moment I would opt for self.S0_tissue_responses for consistency with other functions.
Do you like set_tortuous_parameter(..., use_tissue_S0=True) for triggering it?

Copy link
Collaborator

Choose a reason for hiding this comment

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

ah of course keep things consistent. I just meant that I didn't like the multiple places to give S0_tissue_responses, which would make it possible to do things differently in different places as the same time.

perfect on the option.

Copy link
Collaborator

Choose a reason for hiding this comment

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

See my comment above, we do need to find a way to nicely define which model is associated with what S0 response for intra or extra.

Perhaps we can make the signature even simpler and do def set_tortuous_parameter(model_intra, model_extra, use_tissue_S0=True), and accept the model definitions themselves as input parameters, like

stick=C1Stick()
zepp = G2Zeppelin()
watson = SD1WatsonDistributedModel(models = [stick, sepp], S0s=[1, 2])
watson.set_tortuous_parameter(model_intra=stick, model_extra=zepp, use_tissue_S0=True)

Internally we can just check for the appropriate model with isinstance and get the names for the parameters and associated S0s... what do you think?

"""
Allows the user to set a tortuosity constraint on the perpendicular
diffusivity of the extra-axonal compartment, which depends on the
Expand All @@ -295,6 +296,8 @@ def set_tortuous_parameter(self, lambda_perp,
volume_fraction_intra: string
name of the intra-axonal volume fraction parameter, see
self.parameter_names.
S0_responses: list
list of the s0 responses of the tissues associated to IC and EC.
"""
params = [lambda_perp, lambda_par, volume_fraction_intra]
for param in params:
Expand All @@ -308,7 +311,10 @@ def set_tortuous_parameter(self, lambda_perp,
model, name = self._parameter_map[lambda_perp]
self.parameter_links.append([model, name, T1_tortuosity, [
self._parameter_map[lambda_par],
self._parameter_map[volume_fraction_intra]]
self._parameter_map[volume_fraction_intra],
None,
S0_responses
]
])
del self.parameter_ranges[lambda_perp]
del self.parameter_scales[lambda_perp]
Expand Down
20 changes: 13 additions & 7 deletions dmipy/utils/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ def rotation_matrix_100_to_theta_phi_psi(theta, phi, psi):
return np.dot(R_100_to_theta_phi, R_around_100)


def T1_tortuosity(lambda_par, vf_intra, vf_extra=None):
def T1_tortuosity(lambda_par, vf_intra, vf_extra=None, S0_responses=None):
"""Tortuosity model for perpendicular extra-axonal diffusivity [1, 2, 3].
If vf_extra=None, then vf_intra must be a nested volume fraction, in the
sense that E_bundle = vf_intra * E_intra + (1 - vf_intra) * E_extra, with
Expand All @@ -186,6 +186,8 @@ def T1_tortuosity(lambda_par, vf_intra, vf_extra=None):
intra-axonal volume fraction [0, 1].
vf_extra : float, (optional)
extra-axonal volume fraction [0, 1].
S0_responses : list, (optional)
Copy link
Collaborator

Choose a reason for hiding this comment

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

indicate the length of the list

s0 response of the tissues associated to IC and EC.

Returns
-------
Expand All @@ -204,12 +206,16 @@ def T1_tortuosity(lambda_par, vf_intra, vf_extra=None):
.. [3] Szafer et al. "Theoretical model for water diffusion in tissues."
Magnetic resonance in medicine 33.5 (1995): 697-712.
"""
if vf_extra is None:
lambda_perp = (1 - vf_intra) * lambda_par
else:
fraction_intra = vf_intra / (vf_intra + vf_extra)
lambda_perp = (1 - fraction_intra) * lambda_par
return lambda_perp
fraction_intra = vf_intra
fraction_extra = 1.0 - fraction_intra if vf_extra is None else vf_extra

if S0_responses is not None:
Copy link
Collaborator

Choose a reason for hiding this comment

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

add detailed comment either here or in the doctring of the tortuosity function here - add a preliminary citation to your abstract if you want.

fraction_intra /= S0_responses[0]
fraction_extra /= S0_responses[1]

fraction_intra /= fraction_extra + fraction_intra

return (1. - fraction_intra) * lambda_par


def parameter_equality(param):
Expand Down