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.

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

FEAT: add skeleton for amico implementation #89

Open
wants to merge 31 commits into
base: master
Choose a base branch
from

Conversation

matteofrigo
Copy link
Contributor

@matteofrigo matteofrigo commented Apr 30, 2020

With this PR we add the skeleton for the implementation of AMICO. Let's open the discussion @Sara04 @rutgerfick .

dmipy.core.modeling_framework

For the MultiCompartmentAMICOModel and the MultiCompartmentAMICOSphericalMeanModel we have to do the following:

  • Implement _create_forward_model_matrix
  • Choose a strategy for the solver to use. The difficulty stands in keeping the syntax consistent with the other models
  • Update documentation of fit function
  • Complete the setup of the optimization before model fitting in the fit function, taking into account that the look-up table must be built or loaded.
  • Implement the fit_func method in the fit function
  • Preprocess the x0 as required by the solver
  • Implement the multi-tissue feature correctly

dmipy.core.fitted_modeling_framework

  • I added the fitted_distribution property to the FittedMultiCompartmentAMICOModel and FittedMultiCompartmentAMICOSphericalMeanModel classes. This is intended to return the distribution of each parameter obtained from the AMICO fitting.
  • The fitted_parameter property could be built on top of the fitted_distribution.

This commit adds the mentioned function, which builds a set of N
directions uniformly dispersed on a unitary hemisphere.
Copy link
Collaborator

@rutgerfick rutgerfick left a comment

Choose a reason for hiding this comment

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

Damn that's a good skeleton @matteofrigo !

dmipy/core/fitted_modeling_framework.py Outdated Show resolved Hide resolved
dmipy/core/fitted_modeling_framework.py Show resolved Hide resolved
dmipy/core/fitted_modeling_framework.py Show resolved Hide resolved
dmipy/core/fitted_modeling_framework.py Outdated Show resolved Hide resolved
dmipy/core/fitted_modeling_framework.py Outdated Show resolved Hide resolved
dmipy/core/tests/test_amico_models.py Show resolved Hide resolved
return model, data


def test_forward_model():
Copy link
Collaborator

Choose a reason for hiding this comment

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

these tests should also include a signal fitting test were we define a distribution by hand and see if AMICO can fit that distribution.

I can already expect that the AMICO fit will be ill-posed if the distribution is densely sampled, so expect that to test the correctness the forward model will have to be quite small (very course sampling of parameters)

dmipy/core/tests/test_amico_spherical_mean_models.py Outdated Show resolved Hide resolved
dmipy/core/tests/test_amico_models.py Show resolved Hide resolved
dmipy/utils/build_sphere.py Outdated Show resolved Hide resolved
@rutgerfick
Copy link
Collaborator

As for the creation of the AMICO solver, I think we should build the skeleton the same as for https://github.com/AthenaEPI/dmipy/blob/master/dmipy/optimizers_fod/csd_cvxpy.py.

So without a lookup-table, I think in the first iteration we just build the M-matrix on-the-fly in the main call of the class.

In a next iteration we prepare the look-up table etc in the init, and then use that in the call to speed up the process.

However, since we're using multi-processing we have to consider that a lookup table is a big piece of memory and it will be painful to keep creating processes with this block of data, but we'll tackle this problem when we get there.

This commit removes the implementation of the skeleton for the spherical
mean version of AMICO.

This commit also changes how the parameters are defined in AMICO, which
now is done through two properties of the devoted class:

* self.parameter_indices tells which indices of the forward model matrix
correspond to the wanted parameter
* self.parameter_ranges tells the grid of parameters corresponding to
the distribution obtained from the fitting
# orientations. Also, the coefficients for each orientation must
# be resampled.
pass
raise NotImplementedError
Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't think the predict function should be nonimplemented in this first MR since it addresses basic functionality of the fitted model.

Since in the first version we will not even have the lookup table, we will have to calculate the forward model matrix no matter if acquisition_scheme=None or something else.

The code you wrote here appropriate for the follow-up MR with look-up table.

Copy link
Collaborator

Choose a reason for hiding this comment

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

let's define if the LUT is part of the first MR or not ;)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

We agreed on no LUT for the moment.

Anyway, we need to think this properly. If we have the parameters fitted on a predefined set of orientations, how do we "register" this to the new acquisition scheme? Do we want to use a spherical harmonics representation and resample it accordingly?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Still, it would be tricky for different shells

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 function should look something like

if acquisition_scheme is None:
    acquisition_scheme = self.scheme
E_shape = self.model.mask.shape + (acquisition_scheme.N_dwi,)
E = np.zeros(E_shape)
x, y, z = np.where(self.mask)
for x_, y_, z_ in zip(x, y, z):
    parameter_array = self.fitted_parameter_vectors[x_, y_, z_]
    M = self.forward_model(xxxxx)  # same way as we instantiated it before the fitting
    E[x_, y_, z_] = M.dot(parameter_array)
return E

raise ValueError(msg)

def _check_if_model_orientations_are_fixed(self):
# TODO: return True if the orientations are fixed, False o/wise.
Copy link
Collaborator

Choose a reason for hiding this comment

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

for the other checks I just raise the error inside these internal _check function, see the function directly above.

This allows me to keep the model.fit() function cleaner, or we can change this if you have a good argument for it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Noted. In the future I think it could make sense to make it boolean in order to check when we have to set the orientations automatically.


@property
def forward_model(self):
"""Return the forward model matrix associated to the AMICO model"""
Copy link
Collaborator

Choose a reason for hiding this comment

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

one can only generate the forward model given the directions of the model at this voxel. How are you including this 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.

It reads the directions from the model parameters, that should be set manually before calling anything that needs the forward model matrix

Copy link
Collaborator

Choose a reason for hiding this comment

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

I get that, but do you want to pre-calculate the forward models for all unique orientations in the fixed orientations??
This can quickly become a huge chunk of data.

I think this function should calculate a forward model for a single set of parameters on-the-fly, this way it can also be used in the optimizer, generating a new forward_model matrix for each incoming voxel, and the same in the fitted model predict function.

This means that the function call should be def forward_model(self, orientation_parameters), since the sampling over the other parameters would be set beforehand. What do you think about this?

every parameter. TODO: ?
maxiter : integer,
for MIX optimization, how many iterations are allowed. TODO: ?
N_sphere_samples : integer,
Copy link
Collaborator

Choose a reason for hiding this comment

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

first MR would generate M matrix on-the-fly right? so this parameter does not exist?

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 generates the matrix the first time it is asked for, then reuses the same. This is thought in such a way that we can keep track at the same time of the forward model matrix and the indices associated to each parameter.

Copy link
Collaborator

Choose a reason for hiding this comment

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

what do you mean by re-using?

matteofrigo and others added 7 commits May 17, 2020 19:10
This commit adds numerous details to the AMICO skeleton. The most
relevant changes concern the check of the fixed orientations, that now
is designed to raise an error, and the hemisphere that will be used as
default set of orientation, that now is defined from the symmetric724
of dipy.
@@ -2228,6 +2233,12 @@ def __init__(self, models, S0_tissue_responses=None, parameter_links=None):
self._prepare_parameter_links()
self._prepare_model_properties()
self._check_for_double_model_class_instances()

# QUESTION: I think we should create multi-compartment model
Copy link
Collaborator

Choose a reason for hiding this comment

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

+1

# QUESTION:
# Should we return simulated signals with b=0 values
# or only diffusion weighted, I think b=0 impacts a lot
# results when solving nnls
Copy link
Collaborator

Choose a reason for hiding this comment

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

why should it matter to nnls?

Copy link
Contributor

Choose a reason for hiding this comment

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

Because b=0 have higher values than diffusion weighted values (and less affected by noise), so it will force that sum of b=0 values of atoms is equal to 1 and in this way will prevent that some atoms are fitted to noise. Also, there is a difference if we use for example 18 b=0 values (as in HCP dataset) or just single b=0 value that is equal to the mean of all b=0 values (that would be 1 after normalization), because the higher the number of b=0 values, nnls will give more importance to fitting b=0 parts of the atoms.

Copy link
Collaborator

Choose a reason for hiding this comment

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

So you vote for keeping the b=0 measurements then, if I understand correctly?

Copy link
Collaborator

Choose a reason for hiding this comment

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

I agree with this also.

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, but maybe it is something that we can investigate

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 we can leave the default as being just sample anything that is in the acquisition scheme.
we can test the effect of having more B0s by just making a scheme with more or less b0s in it.

"""Creates forward model matrix, including parameter tessellation grid
and corresponding indices.
"""
"""
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Double docstring. Just needs to be merged.

Copy link
Collaborator

Choose a reason for hiding this comment

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

what do you mean? you just want to merge this in so you can continue on the other branch?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

No, actually that the two docstrings should be merged in a single docstring.

matteofrigo and others added 19 commits August 18, 2020 16:57
This commit adds the missing interface between AmicoCvxpyOptimizer and
MultiCompartmentAMICOModel in the MultiCompartmentAMICOModel.fit
function.

This commit also adapts the output of the fit function to the dedicated
Fitted model.

Tests are still to be developed.
This commit adds a default value to the regularization parameters in
the AMICO optimization framework, which was not present in the previous
version. The default is set to zero.
Any parameter that is set in an AMICO model, it must be set also in the
underlying MultiCompartmentModel. This commit adds this operation.
This commit removes the forward model matrix from the required inputs of
the FittedMultiCompartmentAMICOModel class. The prediction of the signal
will have to be made dependent on the underlying MC model.
This commit changes the signature of the fit_func of the
MultiCompartmentAMICOModel class, which now requires the voxel-specific
directions to be passed.

The check for orientations to be fixed a priori has been fixed to take
into account when distinct orientations have been fixed for each voxel.

Still to be added:
- Check if some non-directional parameter have been set to
  voxel-specific values (forbidden)
- Create a fitted_parameters volume of the right dimension
…add direction to the array of fitted parameters; update shape of fitted parameters
…bution); raise ValueError if number of atoms in forward model matrix is greater than max_atoms
This commit forbids the user from setting the tortuosity constraint
from the MultiCompartmentAMICOModel interface. It must be set from the
DistributeModel interface.

This is made necessary in order to allow the intra axonal fraction to
be properly sampled in the forward model matrix.

Also, minor PEP8 issues are fixed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants