From 3b8210ad869cc73c14b4ea1875d53d522dae804a Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 30 Apr 2020 15:19:24 +0200 Subject: [PATCH 01/28] FEAT: add skeleton for amico implementation --- dmipy/core/fitted_modeling_framework.py | 626 +++++++++++++++ dmipy/core/modeling_framework.py | 733 +++++++++++++++++- dmipy/core/tests/test_amico_models.py | 99 +++ .../tests/test_amico_spherical_mean_models.py | 100 +++ 4 files changed, 1557 insertions(+), 1 deletion(-) create mode 100644 dmipy/core/tests/test_amico_models.py create mode 100644 dmipy/core/tests/test_amico_spherical_mean_models.py diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index 3970f4c8..434f4380 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -897,3 +897,629 @@ def mean_squared_error(self, data): mse = np.mean((data_ - y_hat) ** 2, axis=-1) mse[~self.mask] = 0 return mse + +class FittedMultiCompartmentAMICOModel: + """ + The FittedMultiCompartmentModel instance contains information about the + original MultiCompartmentModel, the estimated S0 values, the fitting mask + and the fitted model parameters. + + Parameters + ---------- + model : MultiCompartmentAMICOModel instance, + A dmipy MultiCompartmentAMICOModel. + S0 : array of size (Ndata,) or (N_data, N_DWIs), + Array containing the estimated S0 values of the data. If data is 4D, + then S0 is 3D if there is only one TE, and the same 4D size of the data + if there are multiple TEs. + mask : array of size (N_data,), + boolean mask of voxels that were fitted. + fitted_parameters_vector : array of size (N_data, Nparameters), + fitted model parameters array. + """ + + def __init__(self, model, S0, mask, fitted_parameters_vector, + fitted_multi_tissue_fractions_vector=None): + self.model = model + self.S0 = S0 + self.mask = mask + self.fitted_parameters_vector = fitted_parameters_vector + self.fitted_multi_tissue_fractions_vector = ( + fitted_multi_tissue_fractions_vector) + + @property + def fitted_distribution(self): + """Returns the distribution of each parameter.""" + # TODO: we have to find a way to return the distribution of each + # parameter obtained from the optimization + raise NotImplementedError + + @property + def fitted_parameters(self): + "Returns the fitted parameters as a dictionary." + # TODO: this could also rely on the fitted_distribution property + return NotImplementedError + + @property + def fitted_multi_tissue_fractions(self): + "Returns the fitted multi tissue fractions as a dictionary." + return self._return_fitted_multi_tissue_fractions() + + @property + def fitted_multi_tissue_fractions_normalized(self): + "Returns the normalized fitted multi tissue fractions as a dictionary" + return self._return_fitted_multi_tissue_fractions(normalized=True) + + def _return_fitted_multi_tissue_fractions(self, normalized=False): + """ + Returns the multi-tissue estimated volume fractions. + + Parameters + ---------- + normalized: boolean, + whether or not to normalize returned multi-tissue volume fractions. + NOTE: This does not mean the unity constraint was enforced during + the estimation of the fractions - it is just a normalization after + the fact. + + Returns + ------- + mt_partial_volumes: dict, + contains the multi-tissue volume fractions by name. + NOTE: if the MC-model only consisted of 1 model, then the name will + be 'partial_volume_0', but will not have a counterpart in + self.fitted_parameters. + """ + mt_fract_vec = self.fitted_multi_tissue_fractions_vector.copy() + if normalized: + mt_fract_sum = np.sum(mt_fract_vec, axis=-1) + mt_fract_mask = mt_fract_sum > 0 + mt_fract_vec[mt_fract_mask] = ( + mt_fract_vec[mt_fract_mask] / + mt_fract_sum[mt_fract_mask][..., None]) + if self.model.N_models > 1: + fract_names = self.model.partial_volume_names + else: + fract_names = ['partial_volume_0'] + mt_partial_volumes = {} + for i, partial_volume_name in enumerate(fract_names): + mt_partial_volumes[partial_volume_name] = mt_fract_vec[..., i] + return mt_partial_volumes + + @property + def fitted_and_linked_parameters(self): + "Returns the fitted and linked parameters as a dictionary." + fitted_parameters = self.model.parameter_vector_to_parameters( + self.fitted_parameters_vector) + return self.model.add_linked_parameters_to_parameters( + fitted_parameters) + + def fod(self, vertices, visual_odi_lower_bound=0.): + """ + Returns the Fiber Orientation Distribution if it is available. + + Parameters + ---------- + vertices : array of size (Nvertices, 3), + Array of cartesian unit vectors at which to sample the FOD. + visual_odi_lower_bound : float, + gives a lower bound to the Orientation Distribution Index (ODI) of + FODs of Watson and Bingham distributions. This can be useful to + visualize FOD fields where some FODs are extremely sharp. + + Returns + ------- + fods : array of size (Ndata, Nvertices), + the FODs of the fitted model, scaled by volume fraction. + """ + if not self.model.fod_available: + msg = ('FODs not available for current model.') + raise ValueError(msg) + dataset_shape = self.fitted_parameters_vector.shape[:-1] + N_samples = len(vertices) + fods = np.zeros(np.r_[dataset_shape, N_samples]) + mask_pos = np.where(self.mask) + for pos in zip(*mask_pos): + parameters = self.model.parameter_vector_to_parameters( + self.fitted_parameters_vector[pos]) + if visual_odi_lower_bound > 0: + param_keys = parameters.keys() + for key in param_keys: + if key[-3:] == 'odi': + parameters[key] = np.clip(parameters[key], + visual_odi_lower_bound, 1) + fods[pos] = self.model(vertices, quantity='FOD', **parameters) + return fods + + def fod_sh(self, sh_order=8, basis_type=None): + """ + Returns the spherical harmonics coefficients of the Fiber Orientation + Distribution (FOD) if it is available. Uses are 724 spherical + tessellation to do the spherical harmonics transform. + + Parameters + ---------- + sh_order : integer, + the maximum spherical harmonics order of the coefficient expansion. + basis_type : string, + type of spherical harmonics basis to use for the expansion, see + sh_to_sf_matrix for more info. + + Returns + ------- + fods_sh : array of size (Ndata, Ncoefficients), + spherical harmonics coefficients of the FODs, scaled by volume + fraction. + """ + if not self.model.fod_available: + msg = ('FODs not available for current model.') + raise ValueError(msg) + sphere = get_sphere(name='repulsion724') + vertices = sphere.vertices + _, inv_sh_matrix = sh_to_sf_matrix( + sphere, sh_order, basis_type=basis_type, return_inv=True) + fods_sf = self.fod(vertices) + + dataset_shape = self.fitted_parameters_vector.shape[:-1] + number_coef_used = int((sh_order + 2) * (sh_order + 1) // 2) + fods_sh = np.zeros(np.r_[dataset_shape, number_coef_used]) + mask_pos = np.where(self.mask) + for pos in zip(*mask_pos): + fods_sh[pos] = np.dot(inv_sh_matrix.T, fods_sf[pos]) + return fods_sh + + def peaks_spherical(self): + "Returns the peak angles of the model." + mu_params = [] + for name, card in self.model.parameter_cardinality.items(): + if name[-2:] == 'mu' and card == 2: + mu_params.append(self.fitted_parameters[name]) + if len(mu_params) == 0: + msg = ('peaks not available for current model.') + raise ValueError(msg) + if len(mu_params) == 1: + return np.expand_dims(mu_params[0], axis=-2) + return np.concatenate([mu[..., None] for mu in mu_params], axis=-1) + + def peaks_cartesian(self): + "Returns the cartesian peak unit vectors of the model." + peaks_spherical = self.peaks_spherical() + peaks_cartesian = unitsphere2cart_Nd(peaks_spherical) + return peaks_cartesian + + def predict(self, acquisition_scheme=None, S0=None, mask=None): + """ + simulates the dMRI signal of the fitted MultiCompartmentModel for the + estimated model parameters. If no acquisition_scheme is given, then + the same acquisition_scheme that was used for the fitting is used. If + no S0 is given then it is assumed to be the estimated one. If no mask + is given then all voxels are assumed to have been fitted. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy. + S0 : None or float, + Signal intensity without diffusion sensitization. If None, uses + estimated SO from fitting process. If float, uses that value. + mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), + mask of voxels to simulate data at. + + Returns + ------- + predicted_signal : array of size (Ndata, N_DWIS), + predicted DWIs for the given model parameters and acquisition + scheme. + """ + if acquisition_scheme is None: + acquisition_scheme = self.model.scheme + self.model._check_model_params_with_acquisition_params( + acquisition_scheme) + + dataset_shape = self.fitted_parameters_vector.shape[:-1] + if S0 is None: + S0 = self.S0 + elif isinstance(S0, float): + S0 = np.ones(dataset_shape) * S0 + if mask is None: + mask = self.mask + + N_samples = len(acquisition_scheme.bvalues) + + predicted_signal = np.zeros(np.r_[dataset_shape, N_samples]) + mask_pos = np.where(mask) + for pos in zip(*mask_pos): + parameters = self.model.parameter_vector_to_parameters( + self.fitted_parameters_vector[pos]) + predicted_signal[pos] = self.model( + acquisition_scheme, **parameters) * S0[pos] + return predicted_signal + + def R2_coefficient_of_determination(self, data): + "Calculates the R-squared of the model fit." + data_ = data / self.S0[..., None] + + y_hat = self.predict(S0=1.) + y_bar = np.mean(data_, axis=-1) + SStot = np.sum((data_ - y_bar[..., None]) ** 2, axis=-1) + SSres = np.sum((data_ - y_hat) ** 2, axis=-1) + R2 = 1 - SSres / SStot + R2[~self.mask] = 0 + return R2 + + def mean_squared_error(self, data): + "Calculates the mean squared error of the model fit." + if self.model.scheme.TE is None: + data_ = data / self.S0[..., None] + else: + data_ = data / self.S0 + + y_hat = self.predict(S0=1.) + mse = np.mean((data_ - y_hat) ** 2, axis=-1) + mse[~self.mask] = 0 + return mse + +class FittedMultiCompartmentAMICOSphericalMeanModel: + """ + The FittedMultiCompartmentAMICOModel instance contains information about the + original MultiCompartmentModel, the estimated S0 values, the fitting mask + and the fitted model parameters. + + Parameters + ---------- + model : MultiCompartmentModel instance, + A dmipy MultiCompartmentModel. + S0 : array of size (Ndata,) or (N_data, N_DWIs), + Array containing the estimated S0 values of the data. If data is 4D, + then S0 is 3D if there is only one TE, and the same 4D size of the data + if there are multiple TEs. + mask : array of size (N_data,), + boolean mask of voxels that were fitted. + fitted_parameters_vector : array of size (N_data, Nparameters), + fitted model parameters array. + """ + + def __init__(self, model, S0, mask, fitted_parameters_vector, + fitted_multi_tissue_fractions_vector=None): + self.model = model + self.S0 = S0 + self.mask = mask + self.fitted_parameters_vector = fitted_parameters_vector + self.fitted_multi_tissue_fractions_vector = ( + fitted_multi_tissue_fractions_vector) + + @property + def fitted_distribution(self): + """Returns the distribution of each parameter.""" + # TODO: we have to find a way to return the distribution of each + # parameter obtained from the optimization + raise NotImplementedError + + @property + def fitted_parameters(self): + "Returns the fitted parameters as a dictionary." + # TODO: this could also rely on the fitted_distribution property + return NotImplementedError + + @property + def fitted_and_linked_parameters(self): + "Returns the fitted and linked parameters as a dictionary." + fitted_parameters = self.model.parameter_vector_to_parameters( + self.fitted_parameters_vector) + return self.model.add_linked_parameters_to_parameters( + fitted_parameters) + + @property + def fitted_multi_tissue_fractions(self): + "Returns the fitted multi tissue fractions as a dictionary." + return self._return_fitted_multi_tissue_fractions() + + @property + def fitted_multi_tissue_fractions_normalized(self): + "Returns the normalized fitted multi tissue fractions as a dictionary" + return self._return_fitted_multi_tissue_fractions(normalized=True) + + def _return_fitted_multi_tissue_fractions(self, normalized=False): + """ + Returns the multi-tissue estimated volume fractions. + + Parameters + ---------- + normalized: boolean, + whether or not to normalize returned multi-tissue volume fractions. + NOTE: This does not mean the unity constraint was enforced during + the estimation of the fractions - it is just a normalization after + the fact. + + Returns + ------- + mt_partial_volumes: dict, + contains the multi-tissue volume fractions by name. + NOTE: if the MC-model only consisted of 1 model, then the name will + be 'partial_volume_0', but will not have a counterpart in + self.fitted_parameters. + """ + mt_fract_vec = self.fitted_multi_tissue_fractions_vector.copy() + if normalized: + mt_fract_sum = np.sum(mt_fract_vec, axis=-1) + mt_fract_mask = mt_fract_sum > 0 + mt_fract_vec[mt_fract_mask] = ( + mt_fract_vec[mt_fract_mask] / + mt_fract_sum[mt_fract_mask][..., None]) + if self.model.N_models > 1: + fract_names = self.model.partial_volume_names + else: + fract_names = ['partial_volume_0'] + mt_partial_volumes = {} + for i, partial_volume_name in enumerate(fract_names): + mt_partial_volumes[partial_volume_name] = mt_fract_vec[..., i] + return mt_partial_volumes + + def predict(self, acquisition_scheme=None, S0=None, mask=None): + """ + simulates the dMRI signal of the fitted MultiCompartmentModel for the + estimated model parameters. If no acquisition_scheme is given, then + the same acquisition_scheme that was used for the fitting is used. If + no S0 is given then it is assumed to be the estimated one. If no mask + is given then all voxels are assumed to have been fitted. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy. + S0 : None or float, + Signal intensity without diffusion sensitization. If None, uses + estimated SO from fitting process. If float, uses that value. + mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), + mask of voxels to simulate data at. + + Returns + ------- + predicted_signal : array of size (Ndata, N_DWIS), + predicted DWIs for the given model parameters and acquisition + scheme. + """ + if acquisition_scheme is None: + acquisition_scheme = self.model.scheme + self.model._check_model_params_with_acquisition_params( + acquisition_scheme) + + dataset_shape = self.fitted_parameters_vector.shape[:-1] + if S0 is None: + S0 = self.S0 + elif isinstance(S0, float): + S0 = np.ones(dataset_shape) * S0 + if mask is None: + mask = self.mask + + N_samples = len(acquisition_scheme.shell_bvalues) + + predicted_signal = np.zeros(np.r_[dataset_shape, N_samples]) + mask_pos = np.where(mask) + for pos in zip(*mask_pos): + parameters = self.model.parameter_vector_to_parameters( + self.fitted_parameters_vector[pos]) + predicted_signal[pos] = self.model( + acquisition_scheme, **parameters) * S0[pos] + return predicted_signal + + def R2_coefficient_of_determination(self, data): + "Calculates the R-squared of the model fit." + Nshells = len(self.model.scheme.shell_bvalues) + data_ = np.zeros(np.r_[data.shape[:-1], Nshells]) + for pos in zip(*np.where(self.mask)): + data_[pos] = estimate_spherical_mean_multi_shell( + data[pos] / self.S0[pos], self.model.scheme) + + y_hat = self.predict(S0=1.) + y_bar = np.mean(data_, axis=-1) + SStot = np.sum((data_ - y_bar[..., None]) ** 2, axis=-1) + SSres = np.sum((data_ - y_hat) ** 2, axis=-1) + R2 = 1 - SSres / SStot + R2[~self.mask] = 0 + return R2 + + def mean_squared_error(self, data): + "Calculates the mean squared error of the model fit." + Nshells = len(self.model.scheme.shell_bvalues) + data_ = np.zeros(np.r_[data.shape[:-1], Nshells]) + for pos in zip(*np.where(self.mask)): + data_[pos] = estimate_spherical_mean_multi_shell( + data[pos] / self.S0[pos], self.model.scheme) + + y_hat = self.predict(S0=1.) + mse = np.mean((data_ - y_hat) ** 2, axis=-1) + mse[~self.mask] = 0 + return mse + + def return_parametric_fod_model( + self, distribution='watson', Ncompartments=1): + """ + Retuns parametric FOD model using the rotational harmonics of the + fitted spherical mean model as the convolution kernel. It can be called + with any implemented parametric distribution (Watson/Bingham) and for + any number of compartments. + + Internally, the input models to the spherical mean model are given to + a spherically distributed model where the parameter links are replayed + such that the distributed model has the same parameter constraints as + the spherical mean model. This distributed model now represents one + compartment of "bundle". This bundle representation is copied + Ncompartment times and given as input to a MultiCompartmentModel, where + now the non-linear are all linked such that each bundle has the same + convolution kernel. Finally, the FittedSphericalMeanModel parameters + are given as fixed parameters for the kernel (the kernel will not be + fitted while the FOD's distribution parameters are being optimized). + + The function returns a MultiCompartmentModel instance that can be + interacted with as usual to fit dMRI data. + + Parameters + ---------- + distribution: string, + Choice of parametric spherical distribution. + Can be 'watson', or 'bingham'. + Ncompartments: integer, + Number of bundles that will be fitted. Must be larger than zero. + + Returns + ------- + mc_bundles_model: Dmipy MultiCompartmentModel instance, + MultiCompartmentModel instance that can be used to estimate + parametric FODs using the fitted spherical mean model as a kernel. + """ + from .modeling_framework import MultiCompartmentModel + from ..distributions import distribute_models + + if not isinstance(Ncompartments, int) or Ncompartments < 1: + msg = 'Ncompartments must be integer larger or equal to one.' + raise ValueError(msg) + + if distribution == 'watson': + bundle = distribute_models.SD1WatsonDistributed( + self.model.models) + basename = 'SD1WatsonDistributed_' + elif distribution == 'bingham': + bundle = distribute_models.SD2BinghamDistributed( + self.model.models) + basename = 'SD2BinghamDistributed_' + else: + msg = '{} is not a valid distribution choice'.format( + distribute_models) + raise ValueError(msg) + + for link in self.model.parameter_links: + param_to_delete = self.model._inverted_parameter_map[link[0], + link[1]] + if link[2] is T1_tortuosity: + bundle.parameter_links.append( + [link[0], link[1], link[2], link[3][:-1]]) + elif link[2] is fractional_parameter: + new_parameter_name = param_to_delete + '_fraction' + bundle.parameter_ranges.update({new_parameter_name: [0., 1.]}) + bundle.parameter_scales.update({new_parameter_name: 1.}) + bundle.parameter_cardinality.update({new_parameter_name: 1}) + bundle.parameter_types.update({new_parameter_name: 'normal'}) + + bundle._parameter_map.update( + {new_parameter_name: (None, 'fraction')}) + bundle._inverted_parameter_map.update( + {(None, 'fraction'): new_parameter_name}) + + # add parmeter link to fractional parameter + param_larger_than = self.model._inverted_parameter_map[ + link[3][1][0], link[3][1][1]] + + model, name = bundle._parameter_map[param_to_delete] + bundle.parameter_links.append( + [model, name, fractional_parameter, [ + bundle._parameter_map[new_parameter_name], + bundle._parameter_map[param_larger_than]]]) + else: + bundle.parameter_links.append(link) + del bundle.parameter_ranges[param_to_delete] + del bundle.parameter_cardinality[param_to_delete] + del bundle.parameter_scales[param_to_delete] + del bundle.parameter_types[param_to_delete] + + bundles = [bundle.copy() for i in range(Ncompartments)] + mc_bundles_model = MultiCompartmentModel(bundles) + parameter_pairs = [] + for smt_par_name in self.model.parameter_names: + parameters = [] + parameters.append(smt_par_name) + for mc_par_name in mc_bundles_model.parameter_names: + if (mc_par_name.startswith(basename) and + mc_par_name.endswith(smt_par_name)): + parameters.append(mc_par_name) + if len(parameters) > 1: + parameter_pairs.append(parameters) + + for parameters in parameter_pairs: + for i in range(2, Ncompartments + 1): + mc_bundles_model.set_equal_parameter(parameters[1], + parameters[i]) + + for parameters in parameter_pairs: + smt_parameter_name = parameters[0] + mc_parameter_name = parameters[1] + mc_bundles_model.set_fixed_parameter( + mc_parameter_name, + self.fitted_parameters[smt_parameter_name]) + return mc_bundles_model + + def return_spherical_harmonics_fod_model(self, sh_order=8): + """ + Retuns spherical harmonics FOD model using the rotational harmonics of + the fitted spherical mean model as the convolution kernel. + + Internally, the input models to the spherical mean model are given to + a MultiCompartmentSphericalHarmonicsModel where the parameter links are + replayed such that the new model has the same parameter constraints as + the spherical mean model. The FittedSphericalMeanModel parameters + are given as fixed parameters for the kernel (the kernel will not be + fitted while the FOD's coefficients are being optimized). + + The function returns a MultiCompartmentSphericalHarmonicsModel instance + that can be interacted with as usual to fit dMRI data. + + Parameters + ---------- + sh_order: even, positive integer, + Spherical harmonics order of the FODs. + + Returns + ------- + mc_bundles_model: Dmipy MultiCompartmentModel instance, + MultiCompartmentModel instance that can be used to estimate + parametric FODs using the fitted spherical mean model as a kernel. + """ + from .modeling_framework import MultiCompartmentSphericalHarmonicsModel + + if sh_order < 0 or sh_order % 2 != 0: + msg = 'sh_order must be an even, positive integer.' + raise ValueError(msg) + + sh_model = MultiCompartmentSphericalHarmonicsModel( + self.model.models, sh_order=sh_order) + + for link in self.model.parameter_links: + param_to_delete = self.model._inverted_parameter_map[link[0], + link[1]] + if link[2] is T1_tortuosity: + sh_model.parameter_links.append( + [link[0], link[1], link[2], link[3][:-1]]) + elif link[2] is fractional_parameter: + new_parameter_name = param_to_delete + '_fraction' + sh_model.parameter_ranges.update( + {new_parameter_name: [0., 1.]}) + sh_model.parameter_scales.update({new_parameter_name: 1.}) + sh_model.parameter_cardinality.update({new_parameter_name: 1}) + sh_model.parameter_types.update({new_parameter_name: 'normal'}) + + sh_model._parameter_map.update( + {new_parameter_name: (None, 'fraction')}) + sh_model._inverted_parameter_map.update( + {(None, 'fraction'): new_parameter_name}) + + # add parmeter link to fractional parameter + param_larger_than = self.model._inverted_parameter_map[ + link[3][1][0], link[3][1][1]] + + model, name = sh_model._parameter_map[param_to_delete] + sh_model.parameter_links.append( + [model, name, fractional_parameter, [ + sh_model._parameter_map[new_parameter_name], + sh_model._parameter_map[param_larger_than]]]) + else: + sh_model.parameter_links.append(link) + del sh_model.parameter_ranges[param_to_delete] + del sh_model.parameter_cardinality[param_to_delete] + del sh_model.parameter_scales[param_to_delete] + del sh_model.parameter_types[param_to_delete] + del sh_model.parameter_optimization_flags[param_to_delete] + + for smt_par_name in self.model.parameter_names: + sh_model.set_fixed_parameter( + smt_par_name, self.fitted_parameters[smt_par_name]) + return sh_model diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 9f901c4d..19dbe170 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -18,7 +18,9 @@ from .fitted_modeling_framework import ( FittedMultiCompartmentModel, FittedMultiCompartmentSphericalMeanModel, - FittedMultiCompartmentSphericalHarmonicsModel) + FittedMultiCompartmentSphericalHarmonicsModel, + FittedMultiCompartmentAMICOModel, + FittedMultiCompartmentAMICOSphericalMeanModel) from ..optimizers.brute2fine import ( GlobalBruteOptimizer, Brute2FineOptimizer) from ..optimizers_fod.csd_tournier import CsdTournierOptimizer @@ -2191,6 +2193,735 @@ def __call__(self, acquisition_scheme, **kwargs): return E +class MultiCompartmentAMICOModel(MultiCompartmentModelProperties): + def __init__(self, models, S0_tissue_responses=None, parameter_links=None): + self.models = models + self.N_models = len(models) + if S0_tissue_responses is not None: + if len(S0_tissue_responses) != self.N_models: + msg = 'Number of S0_tissue responses {} must be same as '\ + 'number of input models {}.' + raise ValueError( + msg.format(len(S0_tissue_responses), self.N_models)) + self.S0_tissue_responses = S0_tissue_responses + self.parameter_links = parameter_links + if parameter_links is None: + self.parameter_links = [] + + self._prepare_parameters() + self._prepare_partial_volumes() + self._prepare_parameter_links() + self._prepare_model_properties() + self._check_for_double_model_class_instances() + self._prepare_parameters_to_optimize() + self._check_for_NMR_and_other_models() + self.x0_parameters = {} + self._forward_model_matrix = None + + if not have_numba: + msg = "We highly recommend installing numba for faster function " + msg += "execution and model fitting." + print(msg) + if not have_pathos: + msg = "We highly recommend installing pathos to take advantage of " + msg += "multicore processing." + print(msg) + + def _check_for_NMR_and_other_models(self): + model_types = [model._model_type for model in self.models] + if "NMRModel" in model_types: + if len(np.unique(model_types)) > 1: + msg = "Cannot combine 1D-NMR and other 3D model types together" + msg += " into a MultiCompartmentAMICOModel." + raise ValueError(msg) + + def _create_forward_model_matrix(self): + # TODO: implement the function that creates the forward model matrix + raise NotImplementedError + + @property + def forward_model(self): + if self._forward_model_matrix is None: + self._create_forward_model_matrix() + return self._forward_model_matrix + + def fit(self, acquisition_scheme, data, + mask=None, + solver=None, # TODO: define the solver + Ns=5, + maxiter=300, + N_sphere_samples=30, + use_parallel_processing=have_pathos, + number_of_processors=None): + """ The main data fitting function of a MultiCompartmentModel. + + This function can fit it to an N-dimensional dMRI data set, and returns + a FittedMultiCompartmentModel instance that contains the fitted + parameters and other useful functions to study the results. + + No initial guess needs to be given to fit a model, but a partial or + complete initial guess can be given if the user wants to have a + solution that is a local minimum close to that guess. The + parameter_initial_guess input can be created using + parameter_initial_guess_to_parameter_vector(). + + A mask can also be given to exclude voxels from fitting (e.g. voxels + that are outside the brain). If no mask is given then all voxels are + included. + + An optimization approach can be chosen as either 'brute2fine' or 'mix'. + - Choosing brute2fine will first use a brute-force optimization to find + an initial guess for parameters without one, and will then refine the + result using gradient-descent-based optimization. + + Note that given no initial guess will make brute2fine precompute an + global parameter grid that will be re-used for all voxels, which in + many cases is much faster than giving voxel-varying initial condition + that requires a grid to be estimated per voxel. + + - Choosing mix will use the recent MIX algorithm based on separation of + linear and non-linear parameters. MIX first uses a stochastic + algorithm to find the non-linear parameters (non-volume fractions), + then estimates the volume fractions while fixing the estimates of the + non-linear parameters, and then finally refines the solution using + a gradient-descent-based algorithm. + + The fitting process can be readily parallelized using the optional + "pathos" package. If it is installed then it will automatically use it, + but it can be turned off by setting use_parallel_processing=False. The + algorithm will automatically use all cores in the machine, unless + otherwise specified in number_of_processors. + + Data with multiple TE are normalized in separate segments using the + b0-values according that TE. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy. + data : N-dimensional array of size (N_x, N_y, ..., N_dwis), + The measured DWI signal attenuation array of either a single voxel + or an N-dimensional dataset. + mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), + Optional mask of voxels to be included in the optimization. + solver : string, + TODO: ? + Ns : integer, + for brute optimization, decised how many steps are sampled for + every parameter. TODO: ? + maxiter : integer, + for MIX optimization, how many iterations are allowed. TODO: ? + N_sphere_samples : integer, + for brute optimization, how many spherical orientations are sampled + for 'mu'. TODO: ? + use_parallel_processing : bool, + whether or not to use parallel processing using pathos. + number_of_processors : integer, + number of processors to use for parallel processing. Defaults to + the number of processors in the computer according to cpu_count(). + + Returns + ------- + FittedMultiCompartmentAMICOModel: class instance that contains fitted + parameters. Can be used to recover parameters themselves or other useful + functions. + """ + self._check_tissue_model_acquisition_scheme(acquisition_scheme) + self._check_model_params_with_acquisition_params(acquisition_scheme) + self._check_acquisition_scheme_has_b0s(acquisition_scheme) + self._check_if_volume_fractions_are_fixed() + + # estimate S0 + self.scheme = acquisition_scheme + data_ = np.atleast_2d(data) + if self.scheme.TE is None or len(np.unique(self.scheme.TE)) == 1: + S0 = np.mean(data_[..., self.scheme.b0_mask], axis=-1) + else: # if multiple TE are in the data + S0 = np.ones_like(data_) + for TE_ in self.scheme.shell_TE: + TE_mask = self.scheme.TE == TE_ + TE_b0_mask = np.all([self.scheme.b0_mask, TE_mask], axis=0) + S0[..., TE_mask] = np.mean( + data_[..., TE_b0_mask], axis=-1)[..., None] + + if mask is None: + mask = data_[..., 0] > 0 + else: + mask = np.all([mask, data_[..., 0] > 0], axis=0) + mask_pos = np.where(mask) + + N_parameters = len(self.bounds_for_optimization) + N_voxels = np.sum(mask) + + # make starting parameters and data the same size + x0_ = self.parameter_initial_guess_to_parameter_vector( + **self.x0_parameters) + x0_ = homogenize_x0_to_data( + data_, x0_) + x0_bool = np.all( + np.isnan(x0_), axis=tuple(np.arange(x0_.ndim - 1))) + x0_[..., ~x0_bool] /= self.scales_for_optimization[~x0_bool] + + if use_parallel_processing and not have_pathos: + msg = 'Cannot use parallel processing without pathos.' + raise ValueError(msg) + elif use_parallel_processing and have_pathos: + fitted_parameters_lin = [None] * N_voxels + if number_of_processors is None: + number_of_processors = cpu_count() + pool = pp.ProcessPool(number_of_processors) + print('Using parallel processing with {} workers.'.format( + number_of_processors)) + else: + fitted_parameters_lin = np.empty( + np.r_[N_voxels, N_parameters], dtype=float) + + # TODO: complete the setup of the optimization in model fitting + def fit_func(*args, **kwargs): + # This function will use the self.forward_model property + # TODO: implement the fit_func method for the model fitting + raise NotImplementedError + self.optimizer = fit_func + + start = time() + for idx, pos in enumerate(zip(*mask_pos)): + voxel_E = data_[pos] / S0[pos] + voxel_x0_vector = x0_[pos] + # TODO: preprocess the x0 as required by the solver + fit_args = (voxel_E, voxel_x0_vector) + + if use_parallel_processing: + fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) + else: + fitted_parameters_lin[idx] = fit_func(*fit_args) + if use_parallel_processing: + fitted_parameters_lin = np.array( + [p.get() for p in fitted_parameters_lin]) + pool.close() + pool.join() + pool.clear() + + fitting_time = time() - start + print('Fitting of {} voxels complete in {} seconds.'.format( + len(fitted_parameters_lin), fitting_time)) + print('Average of {} seconds per voxel.'.format( + fitting_time / N_voxels)) + + fitted_mt_fractions = None + if self.S0_tissue_responses: + # TODO: rescale the signal fractions to get the volume fractions + mt_fractions = None + msg = 'Multi-tissue fitting of {} voxels complete in {} seconds.' + print(msg.format(len(mt_fractions), fitting_time)) + fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) + fitted_mt_fractions[mask_pos] = mt_fractions + + fitted_parameters = np.zeros_like(x0_, dtype=float) + fitted_parameters[mask_pos] = ( + fitted_parameters_lin * self.scales_for_optimization) + + return FittedMultiCompartmentAMICOModel( + self, S0, mask, fitted_parameters, fitted_mt_fractions) + + def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): + """ + Function to simulate diffusion data for a given acquisition_scheme + and model parameters for the MultiCompartmentModel. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy + model_parameters_array : 1D array of size (N_parameters) or + N-dimensional array the same size as the data. + The model parameters of the MultiCompartmentModel model. + + Returns + ------- + E_simulated: 1D array of size (N_parameters) or N-dimensional + array the same size as x0. + The simulated signal of the microstructure model. + """ + self._check_model_params_with_acquisition_params(acquisition_scheme) + + Ndata = acquisition_scheme.number_of_measurements + if isinstance(parameters_array_or_dict, np.ndarray): + x0 = parameters_array_or_dict + elif isinstance(parameters_array_or_dict, dict): + x0 = self.parameters_to_parameter_vector( + **parameters_array_or_dict) + + x0_at_least_2d = np.atleast_2d(x0) + x0_2d = x0_at_least_2d.reshape(-1, x0_at_least_2d.shape[-1]) + E_2d = np.empty(np.r_[x0_2d.shape[:-1], Ndata]) + for i, x0_ in enumerate(x0_2d): + parameters = self.parameter_vector_to_parameters(x0_) + E_2d[i] = self(acquisition_scheme, **parameters) + E_simulated = E_2d.reshape( + np.r_[x0_at_least_2d.shape[:-1], Ndata]) + + if x0.ndim == 1: + return np.squeeze(E_simulated) + else: + return E_simulated + + def __call__(self, acquisition_scheme_or_vertices, + quantity="signal", **kwargs): + """ + The MultiCompartmentModel function call for to generate signal + attenuation for a given acquisition scheme and model parameters. + + First, the linked parameters are added to the optimized parameters. + + Then, every model in the MultiCompartmentModel is called with the right + parameters to recover the part of the signal attenuation of that model. + The resulting values are multiplied with the volume fractions and + finally the combined signal attenuation is returned. + + Aside from the signal, the function call can also return the Fiber + Orientation Distributions (FODs) when a dispersed model is used, and + can also return the stochastic cost function for the MIX algorithm. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy. + quantity : string + can be 'signal', 'FOD' or 'stochastic cost function' depending on + the need of the model. + kwargs: keyword arguments to the model parameter values, + Is internally given as **parameter_dictionary. + """ + if quantity == "signal" or quantity == "FOD": + values = 0 + elif quantity == "stochastic cost function": + values = np.empty(( + acquisition_scheme_or_vertices.number_of_measurements, + len(self.models) + )) + counter = 0 + + kwargs = self.add_linked_parameters_to_parameters( + kwargs + ) + if len(self.models) > 1: + partial_volumes = [ + kwargs[p] for p in self.partial_volume_names + ] + else: + partial_volumes = [1.] + + for model_name, model, partial_volume in zip( + self.model_names, self.models, partial_volumes + ): + parameters = {} + for parameter in model.parameter_ranges: + parameter_name = self._inverted_parameter_map[ + (model, parameter) + ] + parameters[parameter] = kwargs.get( + # , self.parameter_defaults.get(parameter_name) + parameter_name + ) + + if quantity == "signal": + values = ( + values + + partial_volume * model( + acquisition_scheme_or_vertices, **parameters) + ) + elif quantity == "FOD": + try: + values = ( + values + + partial_volume * model.fod( + acquisition_scheme_or_vertices, **parameters) + ) + except AttributeError: + continue + elif quantity == "stochastic cost function": + values[:, counter] = model(acquisition_scheme_or_vertices, + **parameters) + counter += 1 + return values + +class MultiCompartmentAMICOSphericalMeanModel(MultiCompartmentModelProperties): + r''' + The MultiCompartmentAMICOSphericalMeanModel class allows to combine any + number of CompartmentModels and DistributedModels into one combined model + that can be used to fit and simulate dMRI data. + + Parameters + ---------- + models : list of N CompartmentModel instances, + the models to combine into the MultiCompartmentModel. + parameter_links : list of iterables (model, parameter name, link function, + argument list), + deprecated, for testing only. + ''' + + def __init__(self, models, S0_tissue_responses=None, parameter_links=None): + self.models = models + self.N_models = len(models) + if S0_tissue_responses is not None: + if len(S0_tissue_responses) != self.N_models: + msg = 'Number of S0_tissue responses {} must be same as '\ + 'number of input models {}.' + raise ValueError( + msg.format(len(S0_tissue_responses), self.N_models)) + self.S0_tissue_responses = S0_tissue_responses + self.parameter_links = parameter_links + if parameter_links is None: + self.parameter_links = [] + + self._check_for_NMR_models() + self._prepare_parameters() + self._delete_orientation_parameters() + self._prepare_partial_volumes() + self._prepare_parameter_links() + self._prepare_model_properties() + self._check_for_double_model_class_instances() + self._prepare_parameters_to_optimize() + self.x0_parameters = {} + self._forward_model_matrix = None + + + if not have_numba: + msg = "We highly recommend installing numba for faster function " + msg += "execution and model fitting." + print(msg) + if not have_pathos: + msg = "We highly recommend installing pathos to take advantage of " + msg += "multicore processing." + print(msg) + + def _check_for_NMR_models(self): + for model in self.models: + if model._model_type == 'NMRModel': + msg = "Cannot estimate spherical mean of 1D-NMR models." + raise ValueError(msg) + + def _delete_orientation_parameters(self): + """ + Deletes orientation parameters from input models 'mu' since they're not + needed in spherical mean models. + """ + "Removes orientation parameters from input models." + for model in self.models: + for param_name, param_type in model.parameter_types.items(): + if param_type == 'orientation': + appended_param_name = self._inverted_parameter_map[ + model, param_name] + del self.parameter_ranges[appended_param_name] + del self.parameter_scales[appended_param_name] + del self.parameter_cardinality[appended_param_name] + del self.parameter_types[appended_param_name] + + def fit(self, acquisition_scheme, data, + mask=None, + solver=None, # TODO: define the solver + Ns=5, + maxiter=300, + N_sphere_samples=30, + use_parallel_processing=have_pathos, + number_of_processors=None): + """ The main data fitting function of a MultiCompartmentModel. + + This function can fit it to an N-dimensional dMRI data set, and returns + a FittedMultiCompartmentModel instance that contains the fitted + parameters and other useful functions to study the results. + + No initial guess needs to be given to fit a model, but a partial or + complete initial guess can be given if the user wants to have a + solution that is a local minimum close to that guess. The + parameter_initial_guess input can be created using + parameter_initial_guess_to_parameter_vector(). + + A mask can also be given to exclude voxels from fitting (e.g. voxels + that are outside the brain). If no mask is given then all voxels are + included. + + An optimization approach can be chosen as either 'brute2fine' or 'mix'. + - Choosing brute2fine will first use a brute-force optimization to find + an initial guess for parameters without one, and will then refine the + result using gradient-descent-based optimization. + + Note that given no initial guess will make brute2fine precompute an + global parameter grid that will be re-used for all voxels, which in + many cases is much faster than giving voxel-varying initial condition + that requires a grid to be estimated per voxel. + + - Choosing mix will use the recent MIX algorithm based on separation of + linear and non-linear parameters. MIX first uses a stochastic + algorithm to find the non-linear parameters (non-volume fractions), + then estimates the volume fractions while fixing the estimates of the + non-linear parameters, and then finally refines the solution using + a gradient-descent-based algorithm. + + The fitting process can be readily parallelized using the optional + "pathos" package. If it is installed then it will automatically use it, + but it can be turned off by setting use_parallel_processing=False. The + algorithm will automatically use all cores in the machine, unless + otherwise specified in number_of_processors. + + Data with multiple TE are normalized in separate segments using the + b0-values according that TE. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy. + data : N-dimensional array of size (N_x, N_y, ..., N_dwis), + The measured DWI signal attenuation array of either a single voxel + or an N-dimensional dataset. + mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), + Optional mask of voxels to be included in the optimization. + solver : string, + Selection of optimization algorithm. + - 'brute2fine' to use brute-force optimization. + - 'mix' to use Microstructure Imaging of Crossing (MIX) + optimization. + Ns : integer, + for brute optimization, decised how many steps are sampled for + every parameter. + maxiter : integer, + for MIX optimization, how many iterations are allowed. + N_sphere_samples : integer, + for brute optimization, how many spherical orientations are sampled + for 'mu'. + use_parallel_processing : bool, + whether or not to use parallel processing using pathos. + number_of_processors : integer, + number of processors to use for parallel processing. Defaults to + the number of processors in the computer according to cpu_count(). + + Returns + ------- + FittedCompartmentModel: class instance that contains fitted parameters, + Can be used to recover parameters themselves or other useful + functions. + """ + self._check_tissue_model_acquisition_scheme(acquisition_scheme) + self._check_model_params_with_acquisition_params(acquisition_scheme) + self._check_acquisition_scheme_has_b0s(acquisition_scheme) + self._check_if_volume_fractions_are_fixed() + + # estimate S0 + self.scheme = acquisition_scheme + data_ = np.atleast_2d(data) + if self.scheme.TE is None or len(np.unique(self.scheme.TE)) == 1: + S0 = np.mean(data_[..., self.scheme.b0_mask], axis=-1) + else: # if multiple TE are in the data + S0 = np.ones(np.r_[data_.shape[:-1], + len(acquisition_scheme.shell_TE)]) + for TE_ in self.scheme.shell_TE: + TE_mask = self.scheme.shell_TE == TE_ + TE_mask_shell = self.scheme.TE == TE_ + TE_b0_mask = np.all([self.scheme.b0_mask, TE_mask_shell], + axis=0) + S0[..., TE_mask] = np.mean( + data_[..., TE_b0_mask], axis=-1)[..., None] + + if mask is None: + mask = data_[..., 0] > 0 + else: + mask = np.all([mask, data_[..., 0] > 0], axis=0) + mask_pos = np.where(mask) + + N_parameters = len(self.bounds_for_optimization) + N_voxels = np.sum(mask) + + # make starting parameters and data the same size + x0_ = self.parameter_initial_guess_to_parameter_vector( + **self.x0_parameters) + x0_ = homogenize_x0_to_data( + data_, x0_) + x0_bool = np.all( + np.isnan(x0_), axis=tuple(np.arange(x0_.ndim - 1))) + x0_[..., ~x0_bool] /= self.scales_for_optimization[~x0_bool] + + if use_parallel_processing and not have_pathos: + msg = 'Cannot use parallel processing without pathos.' + raise ValueError(msg) + elif use_parallel_processing and have_pathos: + fitted_parameters_lin = [None] * N_voxels + if number_of_processors is None: + number_of_processors = cpu_count() + pool = pp.ProcessPool(number_of_processors) + print('Using parallel processing with {} workers.'.format( + number_of_processors)) + else: + fitted_parameters_lin = np.empty( + np.r_[N_voxels, N_parameters], dtype=float) + + # estimate the spherical mean of the data. + data_to_fit = np.zeros(np.r_[data_.shape[:-1], self.scheme.N_shells]) + for pos in zip(*mask_pos): + data_to_fit[pos] = estimate_spherical_mean_multi_shell( + data_[pos], self.scheme) + + # TODO: complete the setup of the optimization in model fitting + def fit_func(*args, **kwargs): + # This function will use the self.forward_model property + # TODO: implement the fit_func method for the model fitting + raise NotImplementedError + self.optimizer = fit_func + + start = time() + for idx, pos in enumerate(zip(*mask_pos)): + voxel_E = data_to_fit[pos] / S0[pos] + voxel_x0_vector = x0_[pos] + # TODO: preprocess the x0 as required by the solver + fit_args = (voxel_E, voxel_x0_vector) + + if use_parallel_processing: + fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) + else: + fitted_parameters_lin[idx] = fit_func(*fit_args) + if use_parallel_processing: + fitted_parameters_lin = np.array( + [p.get() for p in fitted_parameters_lin]) + pool.close() + pool.join() + pool.clear() + + fitting_time = time() - start + print('Fitting of {} voxels complete in {} seconds.'.format( + len(fitted_parameters_lin), fitting_time)) + print('Average of {} seconds per voxel.'.format( + fitting_time / N_voxels)) + + fitted_mt_fractions = None + if self.S0_tissue_responses: + # TODO: rescale the signal fractions to get the volume fractions + mt_fractions = None + msg = 'Multi-tissue fitting of {} voxels complete in {} seconds.' + print(msg.format(len(mt_fractions), fitting_time)) + fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) + fitted_mt_fractions[mask_pos] = mt_fractions + + fitted_parameters = np.zeros_like(x0_, dtype=float) + fitted_parameters[mask_pos] = ( + fitted_parameters_lin * self.scales_for_optimization) + + return FittedMultiCompartmentAMICOSphericalMeanModel( + self, S0, mask, fitted_parameters, fitted_mt_fractions) + + def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): + """ + Function to simulate diffusion data for a given acquisition_scheme + and model parameters for the MultiCompartmentModel. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy + model_parameters_array : 1D array of size (N_parameters) or + N-dimensional array the same size as the data. + The model parameters of the MultiCompartmentModel model. + + Returns + ------- + E_simulated: 1D array of size (N_parameters) or N-dimensional + array the same size as x0. + The simulated signal of the microstructure model. + """ + self._check_model_params_with_acquisition_params(acquisition_scheme) + + Ndata = acquisition_scheme.shell_indices.max() + 1 + if isinstance(parameters_array_or_dict, np.ndarray): + x0 = parameters_array_or_dict + elif isinstance(parameters_array_or_dict, dict): + x0 = self.parameters_to_parameter_vector( + **parameters_array_or_dict) + + x0_at_least_2d = np.atleast_2d(x0) + x0_2d = x0_at_least_2d.reshape(-1, x0_at_least_2d.shape[-1]) + E_2d = np.empty(np.r_[x0_2d.shape[:-1], Ndata]) + for i, x0_ in enumerate(x0_2d): + parameters = self.parameter_vector_to_parameters(x0_) + E_2d[i] = self(acquisition_scheme, **parameters) + E_simulated = E_2d.reshape( + np.r_[x0_at_least_2d.shape[:-1], Ndata]) + + if x0.ndim == 1: + return np.squeeze(E_simulated) + else: + return E_simulated + + def __call__(self, acquisition_scheme_or_vertices, + quantity="signal", **kwargs): + """ + The MultiCompartmentModel function call for to generate signal + attenuation for a given acquisition scheme and model parameters. + + First, the linked parameters are added to the optimized parameters. + + Then, every model in the MultiCompartmentModel is called with the right + parameters to recover the part of the signal attenuation of that model. + The resulting values are multiplied with the volume fractions and + finally the combined signal attenuation is returned. + + Aside from the signal, the function call can also return the Fiber + Orientation Distributions (FODs) when a dispersed model is used, and + can also return the stochastic cost function for the MIX algorithm. + + Parameters + ---------- + acquisition_scheme : DmipyAcquisitionScheme instance, + An acquisition scheme that has been instantiated using dMipy. + quantity : string + can be 'signal', 'FOD' or 'stochastic cost function' depending on + the need of the model. + kwargs: keyword arguments to the model parameter values, + Is internally given as **parameter_dictionary. + """ + if quantity == "signal": + values = 0 + elif quantity == "stochastic cost function": + values = np.empty(( + len(acquisition_scheme_or_vertices.shell_bvalues), + len(self.models) + )) + counter = 0 + + kwargs = self.add_linked_parameters_to_parameters( + kwargs + ) + if len(self.models) > 1: + partial_volumes = [ + kwargs[p] for p in self.partial_volume_names + ] + else: + partial_volumes = [1.] + + for model_name, model, partial_volume in zip( + self.model_names, self.models, partial_volumes + ): + parameters = {} + for parameter in model.parameter_ranges: + parameter_name = self._inverted_parameter_map[ + (model, parameter) + ] + parameters[parameter] = kwargs.get( + # , self.parameter_defaults.get(parameter_name) + parameter_name + ) + + if quantity == "signal": + values = ( + values + + partial_volume * model.spherical_mean( + acquisition_scheme_or_vertices, **parameters) + ) + elif quantity == "stochastic cost function": + values[:, counter] = model.spherical_mean( + acquisition_scheme_or_vertices, + **parameters) + counter += 1 + return values + + def homogenize_x0_to_data(data, x0): """ Function that checks if data and initial guess x0 are of the same size. diff --git a/dmipy/core/tests/test_amico_models.py b/dmipy/core/tests/test_amico_models.py new file mode 100644 index 00000000..1fb39e92 --- /dev/null +++ b/dmipy/core/tests/test_amico_models.py @@ -0,0 +1,99 @@ +from dmipy.distributions import distribute_models +from dmipy.signal_models import cylinder_models, gaussian_models +from dmipy.core import modeling_framework +from dmipy.data.saved_acquisition_schemes import ( + wu_minn_hcp_acquisition_scheme) +import numpy as np +from numpy.testing import assert_, assert_raises + +scheme = wu_minn_hcp_acquisition_scheme() + + +def one_compartment_amico_model(): + """ + Provides a simple Stick AMICO-like model. + + This function builds an AMICO model with a single compartment given by a + stick convolved with a watson distribution. + + All the parameters are set randomly. + + Returns: + tuple of length 2 with + * MultiCompartmentAMICOModel instance + * sample data simulated on the hcp acquisition scheme + """ + stick = cylinder_models.C1Stick() + watsonstick = distribute_models.SD1WatsonDistributed( + [stick]) + params = {} + for parameter, card, in watsonstick.parameter_cardinality.items(): + params[parameter] = (np.random.rand(card) * + watsonstick.parameter_scales[parameter]) + + model = modeling_framework.MultiCompartmentAMICOModel([watsonstick]) + data = np.atleast_2d(watsonstick(scheme, **params)) + return model, data + + +def two_compartment_amico_model(): + """ + Provides a simple Stick-and-Ball AMICO-like model. + + This function builds an AMICO model with two compartments given by a + stick convolved with a watson distribution and a ball. + + All the parameters are set randomly. + + Returns: + tuple of length 2 with + * MultiCompartmentAMICOModel instance + * sample data simulated on the hcp acquisition scheme + """ + stick = cylinder_models.C1Stick() + ball = gaussian_models.G1Ball() + watsonstick = distribute_models.SD1WatsonDistributed( + [stick]) + model = modeling_framework.MultiCompartmentAMICOModel([watsonstick, ball]) + + params = {} + for parameter, card, in model.parameter_cardinality.items(): + params[parameter] = (np.random.rand(card) * + model.parameter_scales[parameter]) + + for key, value in params.items(): + model.set_fixed_parameter(key, value) + + data = np.atleast_2d(model(scheme, **params)) + return model, data + + +def test_forward_model(): + # TODO: test that the observation matrix is defined as expected + pass + + +def test_nnls(): + # TODO: test the non-regularized version of the inverse problem + pass + + +def test_l2_regularization(): + # TODO: test the effect of L2 regularization + pass + + +def test_l1_regularization(): + # TODO: test the effect of L1 regularization + pass + + +def test_l1_and_l2_regularization(): + # TODO: test the effect of using both L1 and L2 regularization + pass + + +def test_fitted_amico_model_properties(): + # TODO: check that the FittedMultiCompartmentAMICOModel class has the + # expected properties + pass \ No newline at end of file diff --git a/dmipy/core/tests/test_amico_spherical_mean_models.py b/dmipy/core/tests/test_amico_spherical_mean_models.py new file mode 100644 index 00000000..a233d8e0 --- /dev/null +++ b/dmipy/core/tests/test_amico_spherical_mean_models.py @@ -0,0 +1,100 @@ +from dmipy.distributions import distribute_models +from dmipy.signal_models import cylinder_models, gaussian_models +from dmipy.core import modeling_framework +from dmipy.data.saved_acquisition_schemes import ( + wu_minn_hcp_acquisition_scheme) +import numpy as np +from numpy.testing import assert_, assert_raises + +scheme = wu_minn_hcp_acquisition_scheme() + + +def one_compartment_amico_model(): + """ + Provides a simple Stick AMICO-like model. + + This function builds an AMICO model with a single compartment given by a + stick convolved with a watson distribution. + + All the parameters are set randomly. + + Returns: + tuple of length 2 with + * MultiCompartmentAMICOModel instance + * sample data simulated on the hcp acquisition scheme + """ + stick = cylinder_models.C1Stick() + watsonstick = distribute_models.SD1WatsonDistributed( + [stick]) + params = {} + for parameter, card, in watsonstick.parameter_cardinality.items(): + params[parameter] = (np.random.rand(card) * + watsonstick.parameter_scales[parameter]) + + model = modeling_framework.MultiCompartmentAMICOModel([watsonstick]) + data = np.atleast_2d(watsonstick(scheme, **params)) + return model, data + + +def two_compartment_amico_model(): + """ + Provides a simple Stick-and-Ball AMICO-like model. + + This function builds an AMICO model with two compartments given by a + stick convolved with a watson distribution and a ball. + + All the parameters are set randomly. + + Returns: + tuple of length 2 with + * MultiCompartmentAMICOModel instance + * sample data simulated on the hcp acquisition scheme + """ + stick = cylinder_models.C1Stick() + ball = gaussian_models.G1Ball() + watsonstick = distribute_models.SD1WatsonDistributed( + [stick]) + model = modeling_framework.MultiCompartmentAMICOSphericalMeanModel( + [watsonstick, ball]) + + params = {} + for parameter, card, in model.parameter_cardinality.items(): + params[parameter] = (np.random.rand(card) * + model.parameter_scales[parameter]) + + for key, value in params.items(): + model.set_fixed_parameter(key, value) + + data = np.atleast_2d(model(scheme, **params)) + return model, data + + +def test_forward_model(): + # TODO: test that the observation matrix is defined as expected + pass + + +def test_nnls(): + # TODO: test the non-regularized version of the inverse problem + pass + + +def test_l2_regularization(): + # TODO: test the effect of L2 regularization + pass + + +def test_l1_regularization(): + # TODO: test the effect of L1 regularization + pass + + +def test_l1_and_l2_regularization(): + # TODO: test the effect of using both L1 and L2 regularization + pass + + +def test_fitted_amico_model_properties(): + # TODO: check that the FittedMultiCompartmentAMICOSphericalMeanModel class + # has the expected properties + pass \ No newline at end of file From 1a72f0d02933535e2f67120b41f11037096c82c8 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 30 Apr 2020 17:56:45 +0200 Subject: [PATCH 02/28] add dmipy.utils.build_sphere.get_sphere function This commit adds the mentioned function, which builds a set of N directions uniformly dispersed on a unitary hemisphere. --- dmipy/core/modeling_framework.py | 12 ++++++-- dmipy/utils/build_sphere.py | 37 +++++++++++++++++++++++++ dmipy/utils/tests/test_build_sphere.py | 38 ++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 2 deletions(-) create mode 100644 dmipy/utils/build_sphere.py create mode 100644 dmipy/utils/tests/test_build_sphere.py diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 19dbe170..aaa67f21 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -29,6 +29,7 @@ from ..optimizers.multi_tissue_convex_optimizer import ( MultiTissueConvexOptimizer) from dipy.utils.optpkg import optional_package +from dmipy.utils.build_sphere import get_hemisphere from uuid import uuid4 pathos, have_pathos, _ = optional_package("pathos") numba, have_numba, _ = optional_package("numba") @@ -2252,7 +2253,8 @@ def fit(self, acquisition_scheme, data, maxiter=300, N_sphere_samples=30, use_parallel_processing=have_pathos, - number_of_processors=None): + number_of_processors=None, + directions_for_lut=64): """ The main data fitting function of a MultiCompartmentModel. This function can fit it to an N-dimensional dMRI data set, and returns @@ -2319,6 +2321,11 @@ def fit(self, acquisition_scheme, data, number_of_processors : integer, number of processors to use for parallel processing. Defaults to the number of processors in the computer according to cpu_count(). + directions_for_lut : int or N-by-3 array + If integer, this corresponds to the number of points on the + hemisphere that will be used as directions for the look-up-table. + If array, it corresponds to the list of directions that will be + used for building the look up table. Default: 64. Returns ------- @@ -2376,6 +2383,8 @@ def fit(self, acquisition_scheme, data, fitted_parameters_lin = np.empty( np.r_[N_voxels, N_parameters], dtype=float) + # these are the directions to use for the look-up table + directions = get_hemisphere(directions_for_lut) # TODO: complete the setup of the optimization in model fitting def fit_func(*args, **kwargs): # This function will use the self.forward_model property @@ -2694,7 +2703,6 @@ def fit(self, acquisition_scheme, data, number_of_processors : integer, number of processors to use for parallel processing. Defaults to the number of processors in the computer according to cpu_count(). - Returns ------- FittedCompartmentModel: class instance that contains fitted parameters, diff --git a/dmipy/utils/build_sphere.py b/dmipy/utils/build_sphere.py new file mode 100644 index 00000000..c8cd3ed1 --- /dev/null +++ b/dmipy/utils/build_sphere.py @@ -0,0 +1,37 @@ +import numpy as np +from dipy.core.sphere import disperse_charges, HemiSphere + + +def get_hemisphere(directions=1000): + """ + Get a set of directions as a numpy array with 3 columns. + + It relies on the disperse_charges function of dipy. + + Args: + directions: int or N-by-3 array + If integer, this corresponds to the number of points on the + hemisphere that will be used as directions for the look-up-table. + If array, it corresponds to the list of directions that will be + used for building the look up table. + """ + if isinstance(directions, int) or isinstance(directions, float): + # this is a copypaste of the code from the dipy example + n_pts = int(directions) + theta = np.pi * np.random.rand(n_pts) + phi = 2 * np.pi * np.random.rand(n_pts) + hsph_initial = HemiSphere(theta=theta, phi=phi) + hsph_updated, _ = disperse_charges(hsph_initial, 5000) + return hsph_updated.vertices + elif isinstance(directions, np.ndarray) or isinstance(directions, list): + directions = np.squeeze(np.asarray(directions)) + if directions.ndim != 2: + raise ValueError('Directions must be passed as a 2d array.') + if 3 not in directions.shape: + raise ValueError('One of the directions must be 3.') + if directions.shape[0] == 3: + directions = directions.T + directions /= np.linalg.norm(directions, axis=1)[:, None] + return directions + else: + raise TypeError('Input argument must be an integer or a list/array.') \ No newline at end of file diff --git a/dmipy/utils/tests/test_build_sphere.py b/dmipy/utils/tests/test_build_sphere.py new file mode 100644 index 00000000..cb826a7e --- /dev/null +++ b/dmipy/utils/tests/test_build_sphere.py @@ -0,0 +1,38 @@ +from nose.tools import assert_raises, assert_equal, assert_true +from dmipy.utils.build_sphere import get_hemisphere +import numpy as np +from numpy.testing import assert_almost_equal + +def test_get_hemisphere(): + assert_raises(TypeError, get_hemisphere, 'string') + + array1d = np.random.rand(3) + assert_raises(ValueError, get_hemisphere, array1d) + + array_with_4_cols = np.random.rand(10, 4) + assert_raises(ValueError, get_hemisphere, array_with_4_cols) + + # setup dummy transposed set of directions and process it + proper_array_transposed = np.random.rand(3, 10) + proper_array_transposed /= np.linalg.norm(proper_array_transposed, axis=0) + processed = get_hemisphere(proper_array_transposed) + + # check dimensions + assert_equal(processed.shape[1], 3) + # check normalization + assert_almost_equal(np.linalg.norm(processed, axis=1), 1) + + # what if the input is a list? + dirlist = [d for d in processed] + directions = get_hemisphere(dirlist) + assert_true(isinstance(directions, np.ndarray)) + assert_equal(directions.shape[0], len(dirlist)) + assert_equal(directions.shape[1], 3) + + # what if the input is an int? + n = 64 + directions = get_hemisphere(n) + assert_true(isinstance(directions, np.ndarray)) + assert_equal(directions.shape[0], n) + assert_equal(directions.shape[1], 3) + From f12ce31ab6a1ea2ca6e9233e42f1dd6b43758d61 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Sun, 10 May 2020 19:33:35 +0200 Subject: [PATCH 03/28] remove AMICO-SM and add change definition of parameters in AMICO 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 --- dmipy/core/fitted_modeling_framework.py | 431 +----------------- dmipy/core/modeling_framework.py | 400 +--------------- dmipy/core/tests/test_amico_models.py | 4 + .../tests/test_amico_spherical_mean_models.py | 100 ---- 4 files changed, 48 insertions(+), 887 deletions(-) delete mode 100644 dmipy/core/tests/test_amico_spherical_mean_models.py diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index 434f4380..35a21933 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -937,8 +937,20 @@ def fitted_distribution(self): @property def fitted_parameters(self): "Returns the fitted parameters as a dictionary." - # TODO: this could also rely on the fitted_distribution property - return NotImplementedError + fitted = {} + for pname in self.model.parameter_names: + prange = self.model.parameter_ranges[pname] + idx = self.model.parameter_indices[pname] + px = self.fitted_parameters_vector[idx] + value = np.zeros_like(prange[0]) + for v, x in zip(prange, px): # TODO: can be optimized by using dot + value += v * x + value /= np.sum(px, axis=-1) + # By summing over the last axis, we assume the x have been saved + # in a compatible format, namely with one vector per voxel in the + # last direction. + fitted[pname] = value + return fitted @property def fitted_multi_tissue_fractions(self): @@ -1012,24 +1024,7 @@ def fod(self, vertices, visual_odi_lower_bound=0.): fods : array of size (Ndata, Nvertices), the FODs of the fitted model, scaled by volume fraction. """ - if not self.model.fod_available: - msg = ('FODs not available for current model.') - raise ValueError(msg) - dataset_shape = self.fitted_parameters_vector.shape[:-1] - N_samples = len(vertices) - fods = np.zeros(np.r_[dataset_shape, N_samples]) - mask_pos = np.where(self.mask) - for pos in zip(*mask_pos): - parameters = self.model.parameter_vector_to_parameters( - self.fitted_parameters_vector[pos]) - if visual_odi_lower_bound > 0: - param_keys = parameters.keys() - for key in param_keys: - if key[-3:] == 'odi': - parameters[key] = np.clip(parameters[key], - visual_odi_lower_bound, 1) - fods[pos] = self.model(vertices, quantity='FOD', **parameters) - return fods + raise NotImplementedError def fod_sh(self, sh_order=8, basis_type=None): """ @@ -1112,28 +1107,15 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None): scheme. """ if acquisition_scheme is None: - acquisition_scheme = self.model.scheme - self.model._check_model_params_with_acquisition_params( - acquisition_scheme) - - dataset_shape = self.fitted_parameters_vector.shape[:-1] - if S0 is None: - S0 = self.S0 - elif isinstance(S0, float): - S0 = np.ones(dataset_shape) * S0 - if mask is None: - mask = self.mask - - N_samples = len(acquisition_scheme.bvalues) - - predicted_signal = np.zeros(np.r_[dataset_shape, N_samples]) - mask_pos = np.where(mask) - for pos in zip(*mask_pos): - parameters = self.model.parameter_vector_to_parameters( - self.fitted_parameters_vector[pos]) - predicted_signal[pos] = self.model( - acquisition_scheme, **parameters) * S0[pos] - return predicted_signal + # TODO: for each voxel, multiply the forward model matrix by the + # corresponding x. + pass + else: + # TODO: recompute the forward model matrix on the given set of + # orientations. Also, the coefficients for each orientation must + # be resampled. + pass + raise NotImplementedError def R2_coefficient_of_determination(self, data): "Calculates the R-squared of the model fit." @@ -1158,368 +1140,3 @@ def mean_squared_error(self, data): mse = np.mean((data_ - y_hat) ** 2, axis=-1) mse[~self.mask] = 0 return mse - -class FittedMultiCompartmentAMICOSphericalMeanModel: - """ - The FittedMultiCompartmentAMICOModel instance contains information about the - original MultiCompartmentModel, the estimated S0 values, the fitting mask - and the fitted model parameters. - - Parameters - ---------- - model : MultiCompartmentModel instance, - A dmipy MultiCompartmentModel. - S0 : array of size (Ndata,) or (N_data, N_DWIs), - Array containing the estimated S0 values of the data. If data is 4D, - then S0 is 3D if there is only one TE, and the same 4D size of the data - if there are multiple TEs. - mask : array of size (N_data,), - boolean mask of voxels that were fitted. - fitted_parameters_vector : array of size (N_data, Nparameters), - fitted model parameters array. - """ - - def __init__(self, model, S0, mask, fitted_parameters_vector, - fitted_multi_tissue_fractions_vector=None): - self.model = model - self.S0 = S0 - self.mask = mask - self.fitted_parameters_vector = fitted_parameters_vector - self.fitted_multi_tissue_fractions_vector = ( - fitted_multi_tissue_fractions_vector) - - @property - def fitted_distribution(self): - """Returns the distribution of each parameter.""" - # TODO: we have to find a way to return the distribution of each - # parameter obtained from the optimization - raise NotImplementedError - - @property - def fitted_parameters(self): - "Returns the fitted parameters as a dictionary." - # TODO: this could also rely on the fitted_distribution property - return NotImplementedError - - @property - def fitted_and_linked_parameters(self): - "Returns the fitted and linked parameters as a dictionary." - fitted_parameters = self.model.parameter_vector_to_parameters( - self.fitted_parameters_vector) - return self.model.add_linked_parameters_to_parameters( - fitted_parameters) - - @property - def fitted_multi_tissue_fractions(self): - "Returns the fitted multi tissue fractions as a dictionary." - return self._return_fitted_multi_tissue_fractions() - - @property - def fitted_multi_tissue_fractions_normalized(self): - "Returns the normalized fitted multi tissue fractions as a dictionary" - return self._return_fitted_multi_tissue_fractions(normalized=True) - - def _return_fitted_multi_tissue_fractions(self, normalized=False): - """ - Returns the multi-tissue estimated volume fractions. - - Parameters - ---------- - normalized: boolean, - whether or not to normalize returned multi-tissue volume fractions. - NOTE: This does not mean the unity constraint was enforced during - the estimation of the fractions - it is just a normalization after - the fact. - - Returns - ------- - mt_partial_volumes: dict, - contains the multi-tissue volume fractions by name. - NOTE: if the MC-model only consisted of 1 model, then the name will - be 'partial_volume_0', but will not have a counterpart in - self.fitted_parameters. - """ - mt_fract_vec = self.fitted_multi_tissue_fractions_vector.copy() - if normalized: - mt_fract_sum = np.sum(mt_fract_vec, axis=-1) - mt_fract_mask = mt_fract_sum > 0 - mt_fract_vec[mt_fract_mask] = ( - mt_fract_vec[mt_fract_mask] / - mt_fract_sum[mt_fract_mask][..., None]) - if self.model.N_models > 1: - fract_names = self.model.partial_volume_names - else: - fract_names = ['partial_volume_0'] - mt_partial_volumes = {} - for i, partial_volume_name in enumerate(fract_names): - mt_partial_volumes[partial_volume_name] = mt_fract_vec[..., i] - return mt_partial_volumes - - def predict(self, acquisition_scheme=None, S0=None, mask=None): - """ - simulates the dMRI signal of the fitted MultiCompartmentModel for the - estimated model parameters. If no acquisition_scheme is given, then - the same acquisition_scheme that was used for the fitting is used. If - no S0 is given then it is assumed to be the estimated one. If no mask - is given then all voxels are assumed to have been fitted. - - Parameters - ---------- - acquisition_scheme : DmipyAcquisitionScheme instance, - An acquisition scheme that has been instantiated using dMipy. - S0 : None or float, - Signal intensity without diffusion sensitization. If None, uses - estimated SO from fitting process. If float, uses that value. - mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), - mask of voxels to simulate data at. - - Returns - ------- - predicted_signal : array of size (Ndata, N_DWIS), - predicted DWIs for the given model parameters and acquisition - scheme. - """ - if acquisition_scheme is None: - acquisition_scheme = self.model.scheme - self.model._check_model_params_with_acquisition_params( - acquisition_scheme) - - dataset_shape = self.fitted_parameters_vector.shape[:-1] - if S0 is None: - S0 = self.S0 - elif isinstance(S0, float): - S0 = np.ones(dataset_shape) * S0 - if mask is None: - mask = self.mask - - N_samples = len(acquisition_scheme.shell_bvalues) - - predicted_signal = np.zeros(np.r_[dataset_shape, N_samples]) - mask_pos = np.where(mask) - for pos in zip(*mask_pos): - parameters = self.model.parameter_vector_to_parameters( - self.fitted_parameters_vector[pos]) - predicted_signal[pos] = self.model( - acquisition_scheme, **parameters) * S0[pos] - return predicted_signal - - def R2_coefficient_of_determination(self, data): - "Calculates the R-squared of the model fit." - Nshells = len(self.model.scheme.shell_bvalues) - data_ = np.zeros(np.r_[data.shape[:-1], Nshells]) - for pos in zip(*np.where(self.mask)): - data_[pos] = estimate_spherical_mean_multi_shell( - data[pos] / self.S0[pos], self.model.scheme) - - y_hat = self.predict(S0=1.) - y_bar = np.mean(data_, axis=-1) - SStot = np.sum((data_ - y_bar[..., None]) ** 2, axis=-1) - SSres = np.sum((data_ - y_hat) ** 2, axis=-1) - R2 = 1 - SSres / SStot - R2[~self.mask] = 0 - return R2 - - def mean_squared_error(self, data): - "Calculates the mean squared error of the model fit." - Nshells = len(self.model.scheme.shell_bvalues) - data_ = np.zeros(np.r_[data.shape[:-1], Nshells]) - for pos in zip(*np.where(self.mask)): - data_[pos] = estimate_spherical_mean_multi_shell( - data[pos] / self.S0[pos], self.model.scheme) - - y_hat = self.predict(S0=1.) - mse = np.mean((data_ - y_hat) ** 2, axis=-1) - mse[~self.mask] = 0 - return mse - - def return_parametric_fod_model( - self, distribution='watson', Ncompartments=1): - """ - Retuns parametric FOD model using the rotational harmonics of the - fitted spherical mean model as the convolution kernel. It can be called - with any implemented parametric distribution (Watson/Bingham) and for - any number of compartments. - - Internally, the input models to the spherical mean model are given to - a spherically distributed model where the parameter links are replayed - such that the distributed model has the same parameter constraints as - the spherical mean model. This distributed model now represents one - compartment of "bundle". This bundle representation is copied - Ncompartment times and given as input to a MultiCompartmentModel, where - now the non-linear are all linked such that each bundle has the same - convolution kernel. Finally, the FittedSphericalMeanModel parameters - are given as fixed parameters for the kernel (the kernel will not be - fitted while the FOD's distribution parameters are being optimized). - - The function returns a MultiCompartmentModel instance that can be - interacted with as usual to fit dMRI data. - - Parameters - ---------- - distribution: string, - Choice of parametric spherical distribution. - Can be 'watson', or 'bingham'. - Ncompartments: integer, - Number of bundles that will be fitted. Must be larger than zero. - - Returns - ------- - mc_bundles_model: Dmipy MultiCompartmentModel instance, - MultiCompartmentModel instance that can be used to estimate - parametric FODs using the fitted spherical mean model as a kernel. - """ - from .modeling_framework import MultiCompartmentModel - from ..distributions import distribute_models - - if not isinstance(Ncompartments, int) or Ncompartments < 1: - msg = 'Ncompartments must be integer larger or equal to one.' - raise ValueError(msg) - - if distribution == 'watson': - bundle = distribute_models.SD1WatsonDistributed( - self.model.models) - basename = 'SD1WatsonDistributed_' - elif distribution == 'bingham': - bundle = distribute_models.SD2BinghamDistributed( - self.model.models) - basename = 'SD2BinghamDistributed_' - else: - msg = '{} is not a valid distribution choice'.format( - distribute_models) - raise ValueError(msg) - - for link in self.model.parameter_links: - param_to_delete = self.model._inverted_parameter_map[link[0], - link[1]] - if link[2] is T1_tortuosity: - bundle.parameter_links.append( - [link[0], link[1], link[2], link[3][:-1]]) - elif link[2] is fractional_parameter: - new_parameter_name = param_to_delete + '_fraction' - bundle.parameter_ranges.update({new_parameter_name: [0., 1.]}) - bundle.parameter_scales.update({new_parameter_name: 1.}) - bundle.parameter_cardinality.update({new_parameter_name: 1}) - bundle.parameter_types.update({new_parameter_name: 'normal'}) - - bundle._parameter_map.update( - {new_parameter_name: (None, 'fraction')}) - bundle._inverted_parameter_map.update( - {(None, 'fraction'): new_parameter_name}) - - # add parmeter link to fractional parameter - param_larger_than = self.model._inverted_parameter_map[ - link[3][1][0], link[3][1][1]] - - model, name = bundle._parameter_map[param_to_delete] - bundle.parameter_links.append( - [model, name, fractional_parameter, [ - bundle._parameter_map[new_parameter_name], - bundle._parameter_map[param_larger_than]]]) - else: - bundle.parameter_links.append(link) - del bundle.parameter_ranges[param_to_delete] - del bundle.parameter_cardinality[param_to_delete] - del bundle.parameter_scales[param_to_delete] - del bundle.parameter_types[param_to_delete] - - bundles = [bundle.copy() for i in range(Ncompartments)] - mc_bundles_model = MultiCompartmentModel(bundles) - parameter_pairs = [] - for smt_par_name in self.model.parameter_names: - parameters = [] - parameters.append(smt_par_name) - for mc_par_name in mc_bundles_model.parameter_names: - if (mc_par_name.startswith(basename) and - mc_par_name.endswith(smt_par_name)): - parameters.append(mc_par_name) - if len(parameters) > 1: - parameter_pairs.append(parameters) - - for parameters in parameter_pairs: - for i in range(2, Ncompartments + 1): - mc_bundles_model.set_equal_parameter(parameters[1], - parameters[i]) - - for parameters in parameter_pairs: - smt_parameter_name = parameters[0] - mc_parameter_name = parameters[1] - mc_bundles_model.set_fixed_parameter( - mc_parameter_name, - self.fitted_parameters[smt_parameter_name]) - return mc_bundles_model - - def return_spherical_harmonics_fod_model(self, sh_order=8): - """ - Retuns spherical harmonics FOD model using the rotational harmonics of - the fitted spherical mean model as the convolution kernel. - - Internally, the input models to the spherical mean model are given to - a MultiCompartmentSphericalHarmonicsModel where the parameter links are - replayed such that the new model has the same parameter constraints as - the spherical mean model. The FittedSphericalMeanModel parameters - are given as fixed parameters for the kernel (the kernel will not be - fitted while the FOD's coefficients are being optimized). - - The function returns a MultiCompartmentSphericalHarmonicsModel instance - that can be interacted with as usual to fit dMRI data. - - Parameters - ---------- - sh_order: even, positive integer, - Spherical harmonics order of the FODs. - - Returns - ------- - mc_bundles_model: Dmipy MultiCompartmentModel instance, - MultiCompartmentModel instance that can be used to estimate - parametric FODs using the fitted spherical mean model as a kernel. - """ - from .modeling_framework import MultiCompartmentSphericalHarmonicsModel - - if sh_order < 0 or sh_order % 2 != 0: - msg = 'sh_order must be an even, positive integer.' - raise ValueError(msg) - - sh_model = MultiCompartmentSphericalHarmonicsModel( - self.model.models, sh_order=sh_order) - - for link in self.model.parameter_links: - param_to_delete = self.model._inverted_parameter_map[link[0], - link[1]] - if link[2] is T1_tortuosity: - sh_model.parameter_links.append( - [link[0], link[1], link[2], link[3][:-1]]) - elif link[2] is fractional_parameter: - new_parameter_name = param_to_delete + '_fraction' - sh_model.parameter_ranges.update( - {new_parameter_name: [0., 1.]}) - sh_model.parameter_scales.update({new_parameter_name: 1.}) - sh_model.parameter_cardinality.update({new_parameter_name: 1}) - sh_model.parameter_types.update({new_parameter_name: 'normal'}) - - sh_model._parameter_map.update( - {new_parameter_name: (None, 'fraction')}) - sh_model._inverted_parameter_map.update( - {(None, 'fraction'): new_parameter_name}) - - # add parmeter link to fractional parameter - param_larger_than = self.model._inverted_parameter_map[ - link[3][1][0], link[3][1][1]] - - model, name = sh_model._parameter_map[param_to_delete] - sh_model.parameter_links.append( - [model, name, fractional_parameter, [ - sh_model._parameter_map[new_parameter_name], - sh_model._parameter_map[param_larger_than]]]) - else: - sh_model.parameter_links.append(link) - del sh_model.parameter_ranges[param_to_delete] - del sh_model.parameter_cardinality[param_to_delete] - del sh_model.parameter_scales[param_to_delete] - del sh_model.parameter_types[param_to_delete] - del sh_model.parameter_optimization_flags[param_to_delete] - - for smt_par_name in self.model.parameter_names: - sh_model.set_fixed_parameter( - smt_par_name, self.fitted_parameters[smt_par_name]) - return sh_model diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index aaa67f21..60859c6a 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -19,8 +19,7 @@ FittedMultiCompartmentModel, FittedMultiCompartmentSphericalMeanModel, FittedMultiCompartmentSphericalHarmonicsModel, - FittedMultiCompartmentAMICOModel, - FittedMultiCompartmentAMICOSphericalMeanModel) + FittedMultiCompartmentAMICOModel) from ..optimizers.brute2fine import ( GlobalBruteOptimizer, Brute2FineOptimizer) from ..optimizers_fod.csd_tournier import CsdTournierOptimizer @@ -2218,6 +2217,7 @@ def __init__(self, models, S0_tissue_responses=None, parameter_links=None): self._check_for_NMR_and_other_models() self.x0_parameters = {} self._forward_model_matrix = None + self._parameter_indices = {} if not have_numba: msg = "We highly recommend installing numba for faster function " @@ -2236,16 +2236,32 @@ def _check_for_NMR_and_other_models(self): msg += " into a MultiCompartmentAMICOModel." raise ValueError(msg) + def _check_if_model_orientations_are_fixed(self): + # TODO: return True if the orientations are fixed, False o/wise. + raise NotImplementedError + def _create_forward_model_matrix(self): - # TODO: implement the function that creates the forward model matrix + # TODO: implement the function that creates the forward model matrix. + # The function is expected to store the matrix in the + # self._forward_model_matrix field and at the same time it must define + # a dictionary whose keys are the parameter names and the values are + # the set of column indices in the matrix corresponding to the + # parameter. This dictionary must be stored in self._parameter_indices. raise NotImplementedError @property def forward_model(self): + """Return the forward model matrix associated to the AMICO model""" if self._forward_model_matrix is None: self._create_forward_model_matrix() return self._forward_model_matrix + @property + def parameter_indices(self): + """Return the dictionary containing the column indices associated to + each parameter in the forward model matrix""" + return self._parameter_indices + def fit(self, acquisition_scheme, data, mask=None, solver=None, # TODO: define the solver @@ -2496,8 +2512,7 @@ def __call__(self, acquisition_scheme_or_vertices, acquisition_scheme : DmipyAcquisitionScheme instance, An acquisition scheme that has been instantiated using dMipy. quantity : string - can be 'signal', 'FOD' or 'stochastic cost function' depending on - the need of the model. + can be 'signal' or 'FOD' depending on the need of the model. kwargs: keyword arguments to the model parameter values, Is internally given as **parameter_dictionary. """ @@ -2554,381 +2569,6 @@ def __call__(self, acquisition_scheme_or_vertices, counter += 1 return values -class MultiCompartmentAMICOSphericalMeanModel(MultiCompartmentModelProperties): - r''' - The MultiCompartmentAMICOSphericalMeanModel class allows to combine any - number of CompartmentModels and DistributedModels into one combined model - that can be used to fit and simulate dMRI data. - - Parameters - ---------- - models : list of N CompartmentModel instances, - the models to combine into the MultiCompartmentModel. - parameter_links : list of iterables (model, parameter name, link function, - argument list), - deprecated, for testing only. - ''' - - def __init__(self, models, S0_tissue_responses=None, parameter_links=None): - self.models = models - self.N_models = len(models) - if S0_tissue_responses is not None: - if len(S0_tissue_responses) != self.N_models: - msg = 'Number of S0_tissue responses {} must be same as '\ - 'number of input models {}.' - raise ValueError( - msg.format(len(S0_tissue_responses), self.N_models)) - self.S0_tissue_responses = S0_tissue_responses - self.parameter_links = parameter_links - if parameter_links is None: - self.parameter_links = [] - - self._check_for_NMR_models() - self._prepare_parameters() - self._delete_orientation_parameters() - self._prepare_partial_volumes() - self._prepare_parameter_links() - self._prepare_model_properties() - self._check_for_double_model_class_instances() - self._prepare_parameters_to_optimize() - self.x0_parameters = {} - self._forward_model_matrix = None - - - if not have_numba: - msg = "We highly recommend installing numba for faster function " - msg += "execution and model fitting." - print(msg) - if not have_pathos: - msg = "We highly recommend installing pathos to take advantage of " - msg += "multicore processing." - print(msg) - - def _check_for_NMR_models(self): - for model in self.models: - if model._model_type == 'NMRModel': - msg = "Cannot estimate spherical mean of 1D-NMR models." - raise ValueError(msg) - - def _delete_orientation_parameters(self): - """ - Deletes orientation parameters from input models 'mu' since they're not - needed in spherical mean models. - """ - "Removes orientation parameters from input models." - for model in self.models: - for param_name, param_type in model.parameter_types.items(): - if param_type == 'orientation': - appended_param_name = self._inverted_parameter_map[ - model, param_name] - del self.parameter_ranges[appended_param_name] - del self.parameter_scales[appended_param_name] - del self.parameter_cardinality[appended_param_name] - del self.parameter_types[appended_param_name] - - def fit(self, acquisition_scheme, data, - mask=None, - solver=None, # TODO: define the solver - Ns=5, - maxiter=300, - N_sphere_samples=30, - use_parallel_processing=have_pathos, - number_of_processors=None): - """ The main data fitting function of a MultiCompartmentModel. - - This function can fit it to an N-dimensional dMRI data set, and returns - a FittedMultiCompartmentModel instance that contains the fitted - parameters and other useful functions to study the results. - - No initial guess needs to be given to fit a model, but a partial or - complete initial guess can be given if the user wants to have a - solution that is a local minimum close to that guess. The - parameter_initial_guess input can be created using - parameter_initial_guess_to_parameter_vector(). - - A mask can also be given to exclude voxels from fitting (e.g. voxels - that are outside the brain). If no mask is given then all voxels are - included. - - An optimization approach can be chosen as either 'brute2fine' or 'mix'. - - Choosing brute2fine will first use a brute-force optimization to find - an initial guess for parameters without one, and will then refine the - result using gradient-descent-based optimization. - - Note that given no initial guess will make brute2fine precompute an - global parameter grid that will be re-used for all voxels, which in - many cases is much faster than giving voxel-varying initial condition - that requires a grid to be estimated per voxel. - - - Choosing mix will use the recent MIX algorithm based on separation of - linear and non-linear parameters. MIX first uses a stochastic - algorithm to find the non-linear parameters (non-volume fractions), - then estimates the volume fractions while fixing the estimates of the - non-linear parameters, and then finally refines the solution using - a gradient-descent-based algorithm. - - The fitting process can be readily parallelized using the optional - "pathos" package. If it is installed then it will automatically use it, - but it can be turned off by setting use_parallel_processing=False. The - algorithm will automatically use all cores in the machine, unless - otherwise specified in number_of_processors. - - Data with multiple TE are normalized in separate segments using the - b0-values according that TE. - - Parameters - ---------- - acquisition_scheme : DmipyAcquisitionScheme instance, - An acquisition scheme that has been instantiated using dMipy. - data : N-dimensional array of size (N_x, N_y, ..., N_dwis), - The measured DWI signal attenuation array of either a single voxel - or an N-dimensional dataset. - mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), - Optional mask of voxels to be included in the optimization. - solver : string, - Selection of optimization algorithm. - - 'brute2fine' to use brute-force optimization. - - 'mix' to use Microstructure Imaging of Crossing (MIX) - optimization. - Ns : integer, - for brute optimization, decised how many steps are sampled for - every parameter. - maxiter : integer, - for MIX optimization, how many iterations are allowed. - N_sphere_samples : integer, - for brute optimization, how many spherical orientations are sampled - for 'mu'. - use_parallel_processing : bool, - whether or not to use parallel processing using pathos. - number_of_processors : integer, - number of processors to use for parallel processing. Defaults to - the number of processors in the computer according to cpu_count(). - Returns - ------- - FittedCompartmentModel: class instance that contains fitted parameters, - Can be used to recover parameters themselves or other useful - functions. - """ - self._check_tissue_model_acquisition_scheme(acquisition_scheme) - self._check_model_params_with_acquisition_params(acquisition_scheme) - self._check_acquisition_scheme_has_b0s(acquisition_scheme) - self._check_if_volume_fractions_are_fixed() - - # estimate S0 - self.scheme = acquisition_scheme - data_ = np.atleast_2d(data) - if self.scheme.TE is None or len(np.unique(self.scheme.TE)) == 1: - S0 = np.mean(data_[..., self.scheme.b0_mask], axis=-1) - else: # if multiple TE are in the data - S0 = np.ones(np.r_[data_.shape[:-1], - len(acquisition_scheme.shell_TE)]) - for TE_ in self.scheme.shell_TE: - TE_mask = self.scheme.shell_TE == TE_ - TE_mask_shell = self.scheme.TE == TE_ - TE_b0_mask = np.all([self.scheme.b0_mask, TE_mask_shell], - axis=0) - S0[..., TE_mask] = np.mean( - data_[..., TE_b0_mask], axis=-1)[..., None] - - if mask is None: - mask = data_[..., 0] > 0 - else: - mask = np.all([mask, data_[..., 0] > 0], axis=0) - mask_pos = np.where(mask) - - N_parameters = len(self.bounds_for_optimization) - N_voxels = np.sum(mask) - - # make starting parameters and data the same size - x0_ = self.parameter_initial_guess_to_parameter_vector( - **self.x0_parameters) - x0_ = homogenize_x0_to_data( - data_, x0_) - x0_bool = np.all( - np.isnan(x0_), axis=tuple(np.arange(x0_.ndim - 1))) - x0_[..., ~x0_bool] /= self.scales_for_optimization[~x0_bool] - - if use_parallel_processing and not have_pathos: - msg = 'Cannot use parallel processing without pathos.' - raise ValueError(msg) - elif use_parallel_processing and have_pathos: - fitted_parameters_lin = [None] * N_voxels - if number_of_processors is None: - number_of_processors = cpu_count() - pool = pp.ProcessPool(number_of_processors) - print('Using parallel processing with {} workers.'.format( - number_of_processors)) - else: - fitted_parameters_lin = np.empty( - np.r_[N_voxels, N_parameters], dtype=float) - - # estimate the spherical mean of the data. - data_to_fit = np.zeros(np.r_[data_.shape[:-1], self.scheme.N_shells]) - for pos in zip(*mask_pos): - data_to_fit[pos] = estimate_spherical_mean_multi_shell( - data_[pos], self.scheme) - - # TODO: complete the setup of the optimization in model fitting - def fit_func(*args, **kwargs): - # This function will use the self.forward_model property - # TODO: implement the fit_func method for the model fitting - raise NotImplementedError - self.optimizer = fit_func - - start = time() - for idx, pos in enumerate(zip(*mask_pos)): - voxel_E = data_to_fit[pos] / S0[pos] - voxel_x0_vector = x0_[pos] - # TODO: preprocess the x0 as required by the solver - fit_args = (voxel_E, voxel_x0_vector) - - if use_parallel_processing: - fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) - else: - fitted_parameters_lin[idx] = fit_func(*fit_args) - if use_parallel_processing: - fitted_parameters_lin = np.array( - [p.get() for p in fitted_parameters_lin]) - pool.close() - pool.join() - pool.clear() - - fitting_time = time() - start - print('Fitting of {} voxels complete in {} seconds.'.format( - len(fitted_parameters_lin), fitting_time)) - print('Average of {} seconds per voxel.'.format( - fitting_time / N_voxels)) - - fitted_mt_fractions = None - if self.S0_tissue_responses: - # TODO: rescale the signal fractions to get the volume fractions - mt_fractions = None - msg = 'Multi-tissue fitting of {} voxels complete in {} seconds.' - print(msg.format(len(mt_fractions), fitting_time)) - fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) - fitted_mt_fractions[mask_pos] = mt_fractions - - fitted_parameters = np.zeros_like(x0_, dtype=float) - fitted_parameters[mask_pos] = ( - fitted_parameters_lin * self.scales_for_optimization) - - return FittedMultiCompartmentAMICOSphericalMeanModel( - self, S0, mask, fitted_parameters, fitted_mt_fractions) - - def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): - """ - Function to simulate diffusion data for a given acquisition_scheme - and model parameters for the MultiCompartmentModel. - - Parameters - ---------- - acquisition_scheme : DmipyAcquisitionScheme instance, - An acquisition scheme that has been instantiated using dMipy - model_parameters_array : 1D array of size (N_parameters) or - N-dimensional array the same size as the data. - The model parameters of the MultiCompartmentModel model. - - Returns - ------- - E_simulated: 1D array of size (N_parameters) or N-dimensional - array the same size as x0. - The simulated signal of the microstructure model. - """ - self._check_model_params_with_acquisition_params(acquisition_scheme) - - Ndata = acquisition_scheme.shell_indices.max() + 1 - if isinstance(parameters_array_or_dict, np.ndarray): - x0 = parameters_array_or_dict - elif isinstance(parameters_array_or_dict, dict): - x0 = self.parameters_to_parameter_vector( - **parameters_array_or_dict) - - x0_at_least_2d = np.atleast_2d(x0) - x0_2d = x0_at_least_2d.reshape(-1, x0_at_least_2d.shape[-1]) - E_2d = np.empty(np.r_[x0_2d.shape[:-1], Ndata]) - for i, x0_ in enumerate(x0_2d): - parameters = self.parameter_vector_to_parameters(x0_) - E_2d[i] = self(acquisition_scheme, **parameters) - E_simulated = E_2d.reshape( - np.r_[x0_at_least_2d.shape[:-1], Ndata]) - - if x0.ndim == 1: - return np.squeeze(E_simulated) - else: - return E_simulated - - def __call__(self, acquisition_scheme_or_vertices, - quantity="signal", **kwargs): - """ - The MultiCompartmentModel function call for to generate signal - attenuation for a given acquisition scheme and model parameters. - - First, the linked parameters are added to the optimized parameters. - - Then, every model in the MultiCompartmentModel is called with the right - parameters to recover the part of the signal attenuation of that model. - The resulting values are multiplied with the volume fractions and - finally the combined signal attenuation is returned. - - Aside from the signal, the function call can also return the Fiber - Orientation Distributions (FODs) when a dispersed model is used, and - can also return the stochastic cost function for the MIX algorithm. - - Parameters - ---------- - acquisition_scheme : DmipyAcquisitionScheme instance, - An acquisition scheme that has been instantiated using dMipy. - quantity : string - can be 'signal', 'FOD' or 'stochastic cost function' depending on - the need of the model. - kwargs: keyword arguments to the model parameter values, - Is internally given as **parameter_dictionary. - """ - if quantity == "signal": - values = 0 - elif quantity == "stochastic cost function": - values = np.empty(( - len(acquisition_scheme_or_vertices.shell_bvalues), - len(self.models) - )) - counter = 0 - - kwargs = self.add_linked_parameters_to_parameters( - kwargs - ) - if len(self.models) > 1: - partial_volumes = [ - kwargs[p] for p in self.partial_volume_names - ] - else: - partial_volumes = [1.] - - for model_name, model, partial_volume in zip( - self.model_names, self.models, partial_volumes - ): - parameters = {} - for parameter in model.parameter_ranges: - parameter_name = self._inverted_parameter_map[ - (model, parameter) - ] - parameters[parameter] = kwargs.get( - # , self.parameter_defaults.get(parameter_name) - parameter_name - ) - - if quantity == "signal": - values = ( - values + - partial_volume * model.spherical_mean( - acquisition_scheme_or_vertices, **parameters) - ) - elif quantity == "stochastic cost function": - values[:, counter] = model.spherical_mean( - acquisition_scheme_or_vertices, - **parameters) - counter += 1 - return values - def homogenize_x0_to_data(data, x0): """ diff --git a/dmipy/core/tests/test_amico_models.py b/dmipy/core/tests/test_amico_models.py index 1fb39e92..431e6147 100644 --- a/dmipy/core/tests/test_amico_models.py +++ b/dmipy/core/tests/test_amico_models.py @@ -96,4 +96,8 @@ def test_l1_and_l2_regularization(): def test_fitted_amico_model_properties(): # TODO: check that the FittedMultiCompartmentAMICOModel class has the # expected properties + pass + +def test_signal_fitting(): + # TODO: test the accuracy of the fitting on a synthetic signal pass \ No newline at end of file diff --git a/dmipy/core/tests/test_amico_spherical_mean_models.py b/dmipy/core/tests/test_amico_spherical_mean_models.py deleted file mode 100644 index a233d8e0..00000000 --- a/dmipy/core/tests/test_amico_spherical_mean_models.py +++ /dev/null @@ -1,100 +0,0 @@ -from dmipy.distributions import distribute_models -from dmipy.signal_models import cylinder_models, gaussian_models -from dmipy.core import modeling_framework -from dmipy.data.saved_acquisition_schemes import ( - wu_minn_hcp_acquisition_scheme) -import numpy as np -from numpy.testing import assert_, assert_raises - -scheme = wu_minn_hcp_acquisition_scheme() - - -def one_compartment_amico_model(): - """ - Provides a simple Stick AMICO-like model. - - This function builds an AMICO model with a single compartment given by a - stick convolved with a watson distribution. - - All the parameters are set randomly. - - Returns: - tuple of length 2 with - * MultiCompartmentAMICOModel instance - * sample data simulated on the hcp acquisition scheme - """ - stick = cylinder_models.C1Stick() - watsonstick = distribute_models.SD1WatsonDistributed( - [stick]) - params = {} - for parameter, card, in watsonstick.parameter_cardinality.items(): - params[parameter] = (np.random.rand(card) * - watsonstick.parameter_scales[parameter]) - - model = modeling_framework.MultiCompartmentAMICOModel([watsonstick]) - data = np.atleast_2d(watsonstick(scheme, **params)) - return model, data - - -def two_compartment_amico_model(): - """ - Provides a simple Stick-and-Ball AMICO-like model. - - This function builds an AMICO model with two compartments given by a - stick convolved with a watson distribution and a ball. - - All the parameters are set randomly. - - Returns: - tuple of length 2 with - * MultiCompartmentAMICOModel instance - * sample data simulated on the hcp acquisition scheme - """ - stick = cylinder_models.C1Stick() - ball = gaussian_models.G1Ball() - watsonstick = distribute_models.SD1WatsonDistributed( - [stick]) - model = modeling_framework.MultiCompartmentAMICOSphericalMeanModel( - [watsonstick, ball]) - - params = {} - for parameter, card, in model.parameter_cardinality.items(): - params[parameter] = (np.random.rand(card) * - model.parameter_scales[parameter]) - - for key, value in params.items(): - model.set_fixed_parameter(key, value) - - data = np.atleast_2d(model(scheme, **params)) - return model, data - - -def test_forward_model(): - # TODO: test that the observation matrix is defined as expected - pass - - -def test_nnls(): - # TODO: test the non-regularized version of the inverse problem - pass - - -def test_l2_regularization(): - # TODO: test the effect of L2 regularization - pass - - -def test_l1_regularization(): - # TODO: test the effect of L1 regularization - pass - - -def test_l1_and_l2_regularization(): - # TODO: test the effect of using both L1 and L2 regularization - pass - - -def test_fitted_amico_model_properties(): - # TODO: check that the FittedMultiCompartmentAMICOSphericalMeanModel class - # has the expected properties - pass \ No newline at end of file From 3fa0f81a2d1e1752ecb69a2fd4ff9b1cb44805bd Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Sun, 17 May 2020 19:10:53 +0200 Subject: [PATCH 04/28] add details to AMICO skeleton 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. --- dmipy/core/fitted_modeling_framework.py | 33 +---------- dmipy/core/modeling_framework.py | 74 ++++++++++--------------- dmipy/core/tests/test_amico_models.py | 36 ++++++++++++ dmipy/utils/build_sphere.py | 37 ------------- dmipy/utils/tests/test_build_sphere.py | 38 ------------- 5 files changed, 67 insertions(+), 151 deletions(-) delete mode 100644 dmipy/utils/build_sphere.py delete mode 100644 dmipy/utils/tests/test_build_sphere.py diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index 35a21933..d5fedb10 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -1046,41 +1046,14 @@ def fod_sh(self, sh_order=8, basis_type=None): spherical harmonics coefficients of the FODs, scaled by volume fraction. """ - if not self.model.fod_available: - msg = ('FODs not available for current model.') - raise ValueError(msg) - sphere = get_sphere(name='repulsion724') - vertices = sphere.vertices - _, inv_sh_matrix = sh_to_sf_matrix( - sphere, sh_order, basis_type=basis_type, return_inv=True) - fods_sf = self.fod(vertices) - - dataset_shape = self.fitted_parameters_vector.shape[:-1] - number_coef_used = int((sh_order + 2) * (sh_order + 1) // 2) - fods_sh = np.zeros(np.r_[dataset_shape, number_coef_used]) - mask_pos = np.where(self.mask) - for pos in zip(*mask_pos): - fods_sh[pos] = np.dot(inv_sh_matrix.T, fods_sf[pos]) - return fods_sh + raise NotImplementedError def peaks_spherical(self): "Returns the peak angles of the model." - mu_params = [] - for name, card in self.model.parameter_cardinality.items(): - if name[-2:] == 'mu' and card == 2: - mu_params.append(self.fitted_parameters[name]) - if len(mu_params) == 0: - msg = ('peaks not available for current model.') - raise ValueError(msg) - if len(mu_params) == 1: - return np.expand_dims(mu_params[0], axis=-2) - return np.concatenate([mu[..., None] for mu in mu_params], axis=-1) + raise NotImplementedError def peaks_cartesian(self): - "Returns the cartesian peak unit vectors of the model." - peaks_spherical = self.peaks_spherical() - peaks_cartesian = unitsphere2cart_Nd(peaks_spherical) - return peaks_cartesian + raise NotImplementedError def predict(self, acquisition_scheme=None, S0=None, mask=None): """ diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 60859c6a..7c1f4860 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -28,7 +28,6 @@ from ..optimizers.multi_tissue_convex_optimizer import ( MultiTissueConvexOptimizer) from dipy.utils.optpkg import optional_package -from dmipy.utils.build_sphere import get_hemisphere from uuid import uuid4 pathos, have_pathos, _ = optional_package("pathos") numba, have_numba, _ = optional_package("numba") @@ -2194,6 +2193,21 @@ def __call__(self, acquisition_scheme, **kwargs): class MultiCompartmentAMICOModel(MultiCompartmentModelProperties): + """ + The MultiCompartmentAmicoModel class allows to combine any number of + CompartmentModels and DistributedModels into one combined generalized + AMICO model that can be used to fit and simulate dMRI data. + + Parameters + ---------- + models : list of N CompartmentModel instances, + the models to combine into the MultiCompartmentModel. + S0_tissue_responses: list of N values, + the S0 response of the tissue modelled by each compartment. + parameter_links : list of iterables (model, parameter name, link function, + argument list), + deprecated, for testing only. + """ def __init__(self, models, S0_tissue_responses=None, parameter_links=None): self.models = models self.N_models = len(models) @@ -2237,7 +2251,7 @@ def _check_for_NMR_and_other_models(self): raise ValueError(msg) def _check_if_model_orientations_are_fixed(self): - # TODO: return True if the orientations are fixed, False o/wise. + # TODO: raise ValueError if orientations are not fixed raise NotImplementedError def _create_forward_model_matrix(self): @@ -2252,7 +2266,11 @@ def _create_forward_model_matrix(self): @property def forward_model(self): """Return the forward model matrix associated to the AMICO model""" - if self._forward_model_matrix is None: + # TODO: we have to find a way to reset the forward model matrix to None + # whenever a model parameter is changed. O/wise we will have a forward + # matrix stored in the attribute that does not correspond to the + # specified parameters. + if self._forward_model_matrix is None or self.parameter_indices is None: self._create_forward_model_matrix() return self._forward_model_matrix @@ -2260,21 +2278,19 @@ def forward_model(self): def parameter_indices(self): """Return the dictionary containing the column indices associated to each parameter in the forward model matrix""" + if self._forward_model_matrix is None or self.parameter_indices is None: + self._create_forward_model_matrix() return self._parameter_indices def fit(self, acquisition_scheme, data, mask=None, - solver=None, # TODO: define the solver - Ns=5, maxiter=300, - N_sphere_samples=30, use_parallel_processing=have_pathos, - number_of_processors=None, - directions_for_lut=64): + number_of_processors=None): """ The main data fitting function of a MultiCompartmentModel. This function can fit it to an N-dimensional dMRI data set, and returns - a FittedMultiCompartmentModel instance that contains the fitted + a FittedMultiCompartmentAMICOModel instance that contains the fitted parameters and other useful functions to study the results. No initial guess needs to be given to fit a model, but a partial or @@ -2287,23 +2303,6 @@ def fit(self, acquisition_scheme, data, that are outside the brain). If no mask is given then all voxels are included. - An optimization approach can be chosen as either 'brute2fine' or 'mix'. - - Choosing brute2fine will first use a brute-force optimization to find - an initial guess for parameters without one, and will then refine the - result using gradient-descent-based optimization. - - Note that given no initial guess will make brute2fine precompute an - global parameter grid that will be re-used for all voxels, which in - many cases is much faster than giving voxel-varying initial condition - that requires a grid to be estimated per voxel. - - - Choosing mix will use the recent MIX algorithm based on separation of - linear and non-linear parameters. MIX first uses a stochastic - algorithm to find the non-linear parameters (non-volume fractions), - then estimates the volume fractions while fixing the estimates of the - non-linear parameters, and then finally refines the solution using - a gradient-descent-based algorithm. - The fitting process can be readily parallelized using the optional "pathos" package. If it is installed then it will automatically use it, but it can be turned off by setting use_parallel_processing=False. The @@ -2322,26 +2321,14 @@ def fit(self, acquisition_scheme, data, or an N-dimensional dataset. mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), Optional mask of voxels to be included in the optimization. - solver : string, - TODO: ? - Ns : integer, - for brute optimization, decised how many steps are sampled for - every parameter. TODO: ? maxiter : integer, - for MIX optimization, how many iterations are allowed. TODO: ? - N_sphere_samples : integer, - for brute optimization, how many spherical orientations are sampled - for 'mu'. TODO: ? + How many iterations are allowed in the optimization process. + Defaults to 300. use_parallel_processing : bool, whether or not to use parallel processing using pathos. number_of_processors : integer, number of processors to use for parallel processing. Defaults to the number of processors in the computer according to cpu_count(). - directions_for_lut : int or N-by-3 array - If integer, this corresponds to the number of points on the - hemisphere that will be used as directions for the look-up-table. - If array, it corresponds to the list of directions that will be - used for building the look up table. Default: 64. Returns ------- @@ -2353,6 +2340,7 @@ def fit(self, acquisition_scheme, data, self._check_model_params_with_acquisition_params(acquisition_scheme) self._check_acquisition_scheme_has_b0s(acquisition_scheme) self._check_if_volume_fractions_are_fixed() + self._check_if_model_orientations_are_fixed() # estimate S0 self.scheme = acquisition_scheme @@ -2399,8 +2387,6 @@ def fit(self, acquisition_scheme, data, fitted_parameters_lin = np.empty( np.r_[N_voxels, N_parameters], dtype=float) - # these are the directions to use for the look-up table - directions = get_hemisphere(directions_for_lut) # TODO: complete the setup of the optimization in model fitting def fit_func(*args, **kwargs): # This function will use the self.forward_model property @@ -2563,10 +2549,6 @@ def __call__(self, acquisition_scheme_or_vertices, ) except AttributeError: continue - elif quantity == "stochastic cost function": - values[:, counter] = model(acquisition_scheme_or_vertices, - **parameters) - counter += 1 return values diff --git a/dmipy/core/tests/test_amico_models.py b/dmipy/core/tests/test_amico_models.py index 431e6147..6710c547 100644 --- a/dmipy/core/tests/test_amico_models.py +++ b/dmipy/core/tests/test_amico_models.py @@ -68,6 +68,42 @@ def two_compartment_amico_model(): return model, data +def two_sticks_one_ball_amico_model(): + """ + Provides a simple Stick-Stick-and-Ball AMICO-like model. + + This function builds an AMICO model with two compartments given by a + stick convolved with a watson distribution and a ball. + + All the parameters are set randomly. + + Returns: + tuple of length 2 with + * MultiCompartmentAMICOModel instance + * sample data simulated on the hcp acquisition scheme + """ + stick1 = cylinder_models.C1Stick() + stick2 = cylinder_models.C1Stick() + ball = gaussian_models.G1Ball() + watsonstick1 = distribute_models.SD1WatsonDistributed( + [stick1]) + watsonstick2 = distribute_models.SD1WatsonDistributed( + [stick2]) + compartments = [watsonstick1, watsonstick2, ball] + model = modeling_framework.MultiCompartmentAMICOModel(compartments) + + params = {} + for parameter, card, in model.parameter_cardinality.items(): + params[parameter] = (np.random.rand(card) * + model.parameter_scales[parameter]) + + for key, value in params.items(): + model.set_fixed_parameter(key, value) + + data = np.atleast_2d(model(scheme, **params)) + return model, data + + def test_forward_model(): # TODO: test that the observation matrix is defined as expected pass diff --git a/dmipy/utils/build_sphere.py b/dmipy/utils/build_sphere.py deleted file mode 100644 index c8cd3ed1..00000000 --- a/dmipy/utils/build_sphere.py +++ /dev/null @@ -1,37 +0,0 @@ -import numpy as np -from dipy.core.sphere import disperse_charges, HemiSphere - - -def get_hemisphere(directions=1000): - """ - Get a set of directions as a numpy array with 3 columns. - - It relies on the disperse_charges function of dipy. - - Args: - directions: int or N-by-3 array - If integer, this corresponds to the number of points on the - hemisphere that will be used as directions for the look-up-table. - If array, it corresponds to the list of directions that will be - used for building the look up table. - """ - if isinstance(directions, int) or isinstance(directions, float): - # this is a copypaste of the code from the dipy example - n_pts = int(directions) - theta = np.pi * np.random.rand(n_pts) - phi = 2 * np.pi * np.random.rand(n_pts) - hsph_initial = HemiSphere(theta=theta, phi=phi) - hsph_updated, _ = disperse_charges(hsph_initial, 5000) - return hsph_updated.vertices - elif isinstance(directions, np.ndarray) or isinstance(directions, list): - directions = np.squeeze(np.asarray(directions)) - if directions.ndim != 2: - raise ValueError('Directions must be passed as a 2d array.') - if 3 not in directions.shape: - raise ValueError('One of the directions must be 3.') - if directions.shape[0] == 3: - directions = directions.T - directions /= np.linalg.norm(directions, axis=1)[:, None] - return directions - else: - raise TypeError('Input argument must be an integer or a list/array.') \ No newline at end of file diff --git a/dmipy/utils/tests/test_build_sphere.py b/dmipy/utils/tests/test_build_sphere.py deleted file mode 100644 index cb826a7e..00000000 --- a/dmipy/utils/tests/test_build_sphere.py +++ /dev/null @@ -1,38 +0,0 @@ -from nose.tools import assert_raises, assert_equal, assert_true -from dmipy.utils.build_sphere import get_hemisphere -import numpy as np -from numpy.testing import assert_almost_equal - -def test_get_hemisphere(): - assert_raises(TypeError, get_hemisphere, 'string') - - array1d = np.random.rand(3) - assert_raises(ValueError, get_hemisphere, array1d) - - array_with_4_cols = np.random.rand(10, 4) - assert_raises(ValueError, get_hemisphere, array_with_4_cols) - - # setup dummy transposed set of directions and process it - proper_array_transposed = np.random.rand(3, 10) - proper_array_transposed /= np.linalg.norm(proper_array_transposed, axis=0) - processed = get_hemisphere(proper_array_transposed) - - # check dimensions - assert_equal(processed.shape[1], 3) - # check normalization - assert_almost_equal(np.linalg.norm(processed, axis=1), 1) - - # what if the input is a list? - dirlist = [d for d in processed] - directions = get_hemisphere(dirlist) - assert_true(isinstance(directions, np.ndarray)) - assert_equal(directions.shape[0], len(dirlist)) - assert_equal(directions.shape[1], 3) - - # what if the input is an int? - n = 64 - directions = get_hemisphere(n) - assert_true(isinstance(directions, np.ndarray)) - assert_equal(directions.shape[0], n) - assert_equal(directions.shape[1], 3) - From 6a101b255bb58295af70df66233b0d5999375e11 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Wed, 3 Jun 2020 14:04:47 +0200 Subject: [PATCH 05/28] implemented predict and check orientations for AMICO --- dmipy/core/fitted_modeling_framework.py | 49 +++++++++++---- dmipy/core/modeling_framework.py | 81 +++++-------------------- 2 files changed, 54 insertions(+), 76 deletions(-) diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index d5fedb10..d727d944 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -914,11 +914,17 @@ class FittedMultiCompartmentAMICOModel: if there are multiple TEs. mask : array of size (N_data,), boolean mask of voxels that were fitted. + forward_model_matrix : array of size N_DWIs-by-Nparameters, + forward model used in the fitting process. + parameter_indices : dictionary, + keys are parameter names and values are the column indices corresponding + to the parameter. fitted_parameters_vector : array of size (N_data, Nparameters), fitted model parameters array. """ def __init__(self, model, S0, mask, fitted_parameters_vector, + forward_model_matrix, parameter_indices, fitted_multi_tissue_fractions_vector=None): self.model = model self.S0 = S0 @@ -926,6 +932,18 @@ def __init__(self, model, S0, mask, fitted_parameters_vector, self.fitted_parameters_vector = fitted_parameters_vector self.fitted_multi_tissue_fractions_vector = ( fitted_multi_tissue_fractions_vector) + self._forward_model_matrix = forward_model_matrix + self._parameter_indices = parameter_indices + + @property + def forward_model_matrix(self): + """Return forward model matrix.""" + return self._forward_model_matrix + + @property + def parameter_indices(self): + """Return parameter indices.""" + return self._parameter_indices @property def fitted_distribution(self): @@ -1060,8 +1078,7 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None): simulates the dMRI signal of the fitted MultiCompartmentModel for the estimated model parameters. If no acquisition_scheme is given, then the same acquisition_scheme that was used for the fitting is used. If - no S0 is given then it is assumed to be the estimated one. If no mask - is given then all voxels are assumed to have been fitted. + no mask is given then all voxels are assumed to have been fitted. Parameters ---------- @@ -1080,15 +1097,25 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None): scheme. """ if acquisition_scheme is None: - # TODO: for each voxel, multiply the forward model matrix by the - # corresponding x. - pass - else: - # TODO: recompute the forward model matrix on the given set of - # orientations. Also, the coefficients for each orientation must - # be resampled. - pass - raise NotImplementedError + acquisition_scheme = self.model.scheme + self.model._check_model_params_with_acquisition_params( + acquisition_scheme) + + if S0 is None: + S0 = self.S0 + elif isinstance(S0, float): + S0 = np.ones(self.fitted_parameters_vector.shape[:-1]) * S0 + + if mask is None: + mask = self.mask + + E_shape = mask.shape + (acquisition_scheme.N_dwi,) + E = np.zeros(E_shape) + mask_pos = np.where(mask) + for pos in zip(*mask_pos): + parameter_array = self.fitted_parameters_vector[pos] + E[pos] = self.forward_model_matrix.dot(parameter_array) * S0[pos] + return E def R2_coefficient_of_determination(self, data): "Calculates the R-squared of the model fit." diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 7c1f4860..ff30558b 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2230,8 +2230,6 @@ def __init__(self, models, S0_tissue_responses=None, parameter_links=None): self._prepare_parameters_to_optimize() self._check_for_NMR_and_other_models() self.x0_parameters = {} - self._forward_model_matrix = None - self._parameter_indices = {} if not have_numba: msg = "We highly recommend installing numba for faster function " @@ -2251,36 +2249,9 @@ def _check_for_NMR_and_other_models(self): raise ValueError(msg) def _check_if_model_orientations_are_fixed(self): - # TODO: raise ValueError if orientations are not fixed - raise NotImplementedError - - def _create_forward_model_matrix(self): - # TODO: implement the function that creates the forward model matrix. - # The function is expected to store the matrix in the - # self._forward_model_matrix field and at the same time it must define - # a dictionary whose keys are the parameter names and the values are - # the set of column indices in the matrix corresponding to the - # parameter. This dictionary must be stored in self._parameter_indices. - raise NotImplementedError - - @property - def forward_model(self): - """Return the forward model matrix associated to the AMICO model""" - # TODO: we have to find a way to reset the forward model matrix to None - # whenever a model parameter is changed. O/wise we will have a forward - # matrix stored in the attribute that does not correspond to the - # specified parameters. - if self._forward_model_matrix is None or self.parameter_indices is None: - self._create_forward_model_matrix() - return self._forward_model_matrix - - @property - def parameter_indices(self): - """Return the dictionary containing the column indices associated to - each parameter in the forward model matrix""" - if self._forward_model_matrix is None or self.parameter_indices is None: - self._create_forward_model_matrix() - return self._parameter_indices + if 'orientation' in self.parameter_types.values: + msg = 'The orientation parameters must be fixed a priori.' + raise ValueError(msg) def fit(self, acquisition_scheme, data, mask=None, @@ -2431,6 +2402,8 @@ def fit_func(*args, **kwargs): fitted_parameters[mask_pos] = ( fitted_parameters_lin * self.scales_for_optimization) + # TODO: pass the forward model matrix and the parameter indices + # dictionary to the FittedMultiCompartmentAMICOModel return FittedMultiCompartmentAMICOModel( self, S0, mask, fitted_parameters, fitted_mt_fractions) @@ -2479,37 +2452,24 @@ def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): def __call__(self, acquisition_scheme_or_vertices, quantity="signal", **kwargs): """ - The MultiCompartmentModel function call for to generate signal + The MultiCompartmentAMICOModel function call for to generate signal attenuation for a given acquisition scheme and model parameters. First, the linked parameters are added to the optimized parameters. - Then, every model in the MultiCompartmentModel is called with the right - parameters to recover the part of the signal attenuation of that model. - The resulting values are multiplied with the volume fractions and + Then, every model in the MultiCompartmentAMICOModel is called with the + right parameters to recover the part of the signal attenuation of that + model. The resulting values are multiplied with the volume fractions and finally the combined signal attenuation is returned. - Aside from the signal, the function call can also return the Fiber - Orientation Distributions (FODs) when a dispersed model is used, and - can also return the stochastic cost function for the MIX algorithm. - Parameters ---------- acquisition_scheme : DmipyAcquisitionScheme instance, An acquisition scheme that has been instantiated using dMipy. - quantity : string - can be 'signal' or 'FOD' depending on the need of the model. kwargs: keyword arguments to the model parameter values, Is internally given as **parameter_dictionary. """ - if quantity == "signal" or quantity == "FOD": - values = 0 - elif quantity == "stochastic cost function": - values = np.empty(( - acquisition_scheme_or_vertices.number_of_measurements, - len(self.models) - )) - counter = 0 + values = 0 kwargs = self.add_linked_parameters_to_parameters( kwargs @@ -2534,21 +2494,12 @@ def __call__(self, acquisition_scheme_or_vertices, parameter_name ) - if quantity == "signal": - values = ( - values + - partial_volume * model( - acquisition_scheme_or_vertices, **parameters) - ) - elif quantity == "FOD": - try: - values = ( - values + - partial_volume * model.fod( - acquisition_scheme_or_vertices, **parameters) - ) - except AttributeError: - continue + values = ( + values + + partial_volume * model( + acquisition_scheme_or_vertices, **parameters) + ) + return values From 4ef1dbed474e9f785303ac090a4118d36b10c647 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 4 Jun 2020 16:05:57 +0200 Subject: [PATCH 06/28] add comments to prepare amico interface --- dmipy/core/modeling_framework.py | 7 ++++++- dmipy/optimizers/amico_cvxpy.py | 2 ++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index ff30558b..28060670 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2253,6 +2253,12 @@ def _check_if_model_orientations_are_fixed(self): msg = 'The orientation parameters must be fixed a priori.' raise ValueError(msg) + def forward_model_matrix(self, args, *kwargs): + # TODO: move the creation of the forward model matrix from the optimizer + # to here. At the same time, instantiate the parameter grid and + # indices. + raise NotImplementedError + def fit(self, acquisition_scheme, data, mask=None, maxiter=300, @@ -2310,7 +2316,6 @@ def fit(self, acquisition_scheme, data, self._check_tissue_model_acquisition_scheme(acquisition_scheme) self._check_model_params_with_acquisition_params(acquisition_scheme) self._check_acquisition_scheme_has_b0s(acquisition_scheme) - self._check_if_volume_fractions_are_fixed() self._check_if_model_orientations_are_fixed() # estimate S0 diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index dea2aed8..788133b8 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -80,6 +80,8 @@ def __init__(self, acquisition_scheme, model, x0_vector=None, x0_len += m_atoms # Creating parameter tessellation grids and corresponding indices + # TODO: move the matrix/grid/indices definition to the + # modeling_framework module. self.grids, self.idx = {}, {} for m_idx in range(self.model.N_models): model = self.model.models[m_idx] From 0278fa1728eb5e8569518dae748b41df678b5b49 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 4 Jun 2020 17:03:48 +0200 Subject: [PATCH 07/28] add specifications about interfacing amico solver --- dmipy/core/fitted_modeling_framework.py | 55 ++++++++++++----------- dmipy/core/modeling_framework.py | 60 ++++++++++++++++++++----- 2 files changed, 79 insertions(+), 36 deletions(-) diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index d727d944..9e294a7c 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -899,33 +899,35 @@ def mean_squared_error(self, data): return mse class FittedMultiCompartmentAMICOModel: - """ - The FittedMultiCompartmentModel instance contains information about the - original MultiCompartmentModel, the estimated S0 values, the fitting mask - and the fitted model parameters. - - Parameters - ---------- - model : MultiCompartmentAMICOModel instance, - A dmipy MultiCompartmentAMICOModel. - S0 : array of size (Ndata,) or (N_data, N_DWIs), - Array containing the estimated S0 values of the data. If data is 4D, - then S0 is 3D if there is only one TE, and the same 4D size of the data - if there are multiple TEs. - mask : array of size (N_data,), - boolean mask of voxels that were fitted. - forward_model_matrix : array of size N_DWIs-by-Nparameters, - forward model used in the fitting process. - parameter_indices : dictionary, - keys are parameter names and values are the column indices corresponding - to the parameter. - fitted_parameters_vector : array of size (N_data, Nparameters), - fitted model parameters array. - """ - def __init__(self, model, S0, mask, fitted_parameters_vector, - forward_model_matrix, parameter_indices, + forward_model_matrix, parameter_indices, optimizer, fitted_multi_tissue_fractions_vector=None): + """ + The FittedMultiCompartmentModel instance contains information about the + original MultiCompartmentModel, the estimated S0 values, the fitting mask + and the fitted model parameters. + + Parameters + ---------- + model : MultiCompartmentAMICOModel instance, + A dmipy MultiCompartmentAMICOModel. + S0 : array of size (Ndata,) or (N_data, N_DWIs), + Array containing the estimated S0 values of the data. If data is 4D, + then S0 is 3D if there is only one TE, and the same 4D size of the data + if there are multiple TEs. + mask : array of size (N_data,), + boolean mask of voxels that were fitted. + forward_model_matrix : array of size N_DWIs-by-Nparameters, + forward model used in the fitting process. + parameter_indices : dictionary, + keys are parameter names and values are the column indices + corresponding to the parameter. + optimizer : AmicoCvxpyOptimizer instance, + A dmipy AmicoCvxpyOptimizer instantiated object that will be used + for computing the parameter distributions. + fitted_parameters_vector : array of size (N_data, Nparameters), + fitted model parameters array. + """ self.model = model self.S0 = S0 self.mask = mask @@ -934,6 +936,7 @@ def __init__(self, model, S0, mask, fitted_parameters_vector, fitted_multi_tissue_fractions_vector) self._forward_model_matrix = forward_model_matrix self._parameter_indices = parameter_indices + self._optimizer = optimizer @property def forward_model_matrix(self): @@ -949,7 +952,7 @@ def parameter_indices(self): def fitted_distribution(self): """Returns the distribution of each parameter.""" # TODO: we have to find a way to return the distribution of each - # parameter obtained from the optimization + # parameter obtained from the optimization. Uses self._optimizer . raise NotImplementedError @property diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 28060670..8e14ea6b 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -22,6 +22,7 @@ FittedMultiCompartmentAMICOModel) from ..optimizers.brute2fine import ( GlobalBruteOptimizer, Brute2FineOptimizer) +from ..optimizers.amico_cvxpy import AmicoCvxpyOptimizer from ..optimizers_fod.csd_tournier import CsdTournierOptimizer from ..optimizers_fod.csd_cvxpy import CsdCvxpyOptimizer from ..optimizers.mix import MixOptimizer @@ -2230,6 +2231,9 @@ def __init__(self, models, S0_tissue_responses=None, parameter_links=None): self._prepare_parameters_to_optimize() self._check_for_NMR_and_other_models() self.x0_parameters = {} + self._amico_grid = None + self._amico_idx = None + self._freezed_parameters_vector = None if not have_numba: msg = "We highly recommend installing numba for faster function " @@ -2253,10 +2257,29 @@ def _check_if_model_orientations_are_fixed(self): msg = 'The orientation parameters must be fixed a priori.' raise ValueError(msg) - def forward_model_matrix(self, args, *kwargs): + @property + def amico_grid(self): + """Dictionary with parameter names as keys and values that are the + parameter grids used in the definition of the AMICO forward model + matrix.""" + return self._amico_grid + + @property + def amico_idx(self): + """Dictionary with parameter names as keys and values that are the + column indices corresponding to the parameter in the AMICO forward model + matrix.""" + return self._amico_idx + + def forward_model_matrix(self, *args, **kwargs): # TODO: move the creation of the forward model matrix from the optimizer # to here. At the same time, instantiate the parameter grid and # indices. + # - Create the forward model matrix that will be returned + # - Instantiate the self._amico_grid and self._amico_idx attributes + # - Save the "freezed" parameters vector as a hidden attribute and + # check if it has changed since the last call. + # Attribute: self._freezed_parameters_vector raise NotImplementedError def fit(self, acquisition_scheme, data, @@ -2363,11 +2386,17 @@ def fit(self, acquisition_scheme, data, fitted_parameters_lin = np.empty( np.r_[N_voxels, N_parameters], dtype=float) - # TODO: complete the setup of the optimization in model fitting - def fit_func(*args, **kwargs): - # This function will use the self.forward_model property - # TODO: implement the fit_func method for the model fitting + start = time() + opt = AmicoCvxpyOptimizer(self, self.scheme, x0_) # TODO: fix params + # setup self.amico_grids and self.amico_indices + _ = self.forward_model_matrix() + def fit_func(data): + # TODO: extract the parameter value for each param and return the + # parameter vector in the voxel. raise NotImplementedError + print('Setup AMICO optimizer in {} seconds'.format( + time() - start)) + self.optimizer = fit_func start = time() @@ -2375,7 +2404,7 @@ def fit_func(*args, **kwargs): voxel_E = data_[pos] / S0[pos] voxel_x0_vector = x0_[pos] # TODO: preprocess the x0 as required by the solver - fit_args = (voxel_E, voxel_x0_vector) + fit_args = (voxel_E, ) if use_parallel_processing: fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) @@ -2396,8 +2425,18 @@ def fit_func(*args, **kwargs): fitted_mt_fractions = None if self.S0_tissue_responses: - # TODO: rescale the signal fractions to get the volume fractions - mt_fractions = None + # secondary fitting including S0 responses + print('Starting secondary multi-tissue optimization.') + start = time() + mt_fractions = np.empty( + np.r_[N_voxels, self.N_models], dtype=float) + fit_func = MultiTissueConvexOptimizer( + acquisition_scheme, self, self.S0_tissue_responses) + for idx, pos in enumerate(zip(*mask_pos)): + voxel_S = data_[pos] + parameters = fitted_parameters_lin[idx] + mt_fractions[idx] = fit_func(voxel_S, parameters) + fitting_time = time() - start msg = 'Multi-tissue fitting of {} voxels complete in {} seconds.' print(msg.format(len(mt_fractions), fitting_time)) fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) @@ -2407,8 +2446,9 @@ def fit_func(*args, **kwargs): fitted_parameters[mask_pos] = ( fitted_parameters_lin * self.scales_for_optimization) - # TODO: pass the forward model matrix and the parameter indices - # dictionary to the FittedMultiCompartmentAMICOModel + # TODO: pass the forward model matrix, the AmicoCvxpyOptimizer object + # and the parameter indices dictionary to the + # FittedMultiCompartmentAMICOModel . return FittedMultiCompartmentAMICOModel( self, S0, mask, fitted_parameters, fitted_mt_fractions) From b398592098f78facc55f374c174f8088a5e70418 Mon Sep 17 00:00:00 2001 From: Sara04 Date: Mon, 8 Jun 2020 20:43:52 +0200 Subject: [PATCH 08/28] Move creation of forward model matrix --- dmipy/core/modeling_framework.py | 100 +++++++++++++++++++++++++++---- dmipy/optimizers/amico_cvxpy.py | 59 +----------------- 2 files changed, 91 insertions(+), 68 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 8e14ea6b..981540c0 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2205,11 +2205,15 @@ class MultiCompartmentAMICOModel(MultiCompartmentModelProperties): the models to combine into the MultiCompartmentModel. S0_tissue_responses: list of N values, the S0 response of the tissue modelled by each compartment. + Nt: int + Number of equidistant sampling points over each + parameter range used to create atoms of observation matrix M parameter_links : list of iterables (model, parameter name, link function, argument list), deprecated, for testing only. """ - def __init__(self, models, S0_tissue_responses=None, parameter_links=None): + def __init__(self, models, S0_tissue_responses=None, Nt=10, + parameter_links=None): self.models = models self.N_models = len(models) if S0_tissue_responses is not None: @@ -2219,6 +2223,7 @@ def __init__(self, models, S0_tissue_responses=None, parameter_links=None): raise ValueError( msg.format(len(S0_tissue_responses), self.N_models)) self.S0_tissue_responses = S0_tissue_responses + self.Nt = Nt self.parameter_links = parameter_links if parameter_links is None: self.parameter_links = [] @@ -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 + # somewhere in __init__ to avoid creating it every time we create + # forward model matrix + self.mc_model = MultiCompartmentModel(models=self.models) + self._prepare_parameters_to_optimize() self._check_for_NMR_and_other_models() self.x0_parameters = {} @@ -2271,16 +2282,83 @@ def amico_idx(self): matrix.""" return self._amico_idx - def forward_model_matrix(self, *args, **kwargs): - # TODO: move the creation of the forward model matrix from the optimizer - # to here. At the same time, instantiate the parameter grid and - # indices. - # - Create the forward model matrix that will be returned - # - Instantiate the self._amico_grid and self._amico_idx attributes - # - Save the "freezed" parameters vector as a hidden attribute and - # check if it has changed since the last call. - # Attribute: self._freezed_parameters_vector - raise NotImplementedError + def forward_model_matrix(self, acquisition_scheme, model_dirs, **kwargs): + """Creates forward model matrix, including parameter tessellation grid + and corresponding indices. + """ + """ + Arguments: + acquisition_scheme: instance containing acquisition protocol + model_dirs: list containing direction of all models in + multi-compartment model + Returns: + observation matrix M + """ + + dir_params = [p for p in self.mc_model.parameter_names + if p.endswith('mu')] + if len(dir_params) != len(model_dirs): + raise ValueError("Length of model_dirs should correspond " + "to the number of directional parameters!") + + if not self._freezed_parameters_vector: + self._amico_grid, self._amico_idx = {}, {} + + grid_params = \ + [p for p in self.mc_model.parameter_names + if not p.endswith('mu') and + not p.startswith('partial_volume')] + + # Compute length of the vector x0 + x0_len = 0 + for m_idx in range(self.N_models): + m_atoms = 1 + for p in self.models[m_idx].parameter_names: + if self.mc_model.model_names[m_idx] + p in grid_params: + m_atoms *= self.Nt + x0_len += m_atoms + + for m_idx in range(self.N_models): + model = self.models[m_idx] + model_name = self.mc_model.model_names[m_idx] + + param_sampling, grid_params_names = [], [] + m_atoms = 1 + for p in model.parameter_names: + if model_name + p not in grid_params: + continue + grid_params_names.append(model_name + p) + p_range = self.mc_model.parameter_ranges[model_name + p] + self._amico_grid[model_name + p] = \ + np.full(x0_len, np.mean(p_range)) + param_sampling.append(np.linspace(p_range[0], p_range[1], + self.Nt, endpoint=True)) + m_atoms *= self.Nt + + self._amico_idx[model_name] =\ + sum([len(self._amico_idx[k]) + for k in self._amico_idx]) + \ + np.arange(m_atoms) + params_mesh = np.meshgrid(*param_sampling) + for p_idx, p in enumerate(grid_params_names): + self.grids[p][self._amico_idx[model_name]] =\ + np.ravel(params_mesh[p_idx]) + + self._amico_grid['partial_volume_' + str(m_idx)] = \ + np.zeros(x0_len) + self._amico_grid['partial_volume_' + + str(m_idx)][self._amico_idx[model_name]] = 1. + self._freezed_parameters_vector = True + + for d_idx, dp in enumerate(dir_params): + self._amico_grids[dp] = model_dirs[d_idx] + + # 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 + return self.mc_model.simulate_signal(acquisition_scheme, + self._amico_grid).T def fit(self, acquisition_scheme, data, mask=None, diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index 788133b8..8d20d165 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -47,8 +47,8 @@ class AmicoCvxpyOptimizer: Learning Research 17.1 (2016): 2909-2913. """ - def __init__(self, acquisition_scheme, model, x0_vector=None, - lambda_1=None, lambda_2=None, Nt=10): + def __init__(self, acquisition_scheme, model, + lambda_1=None, lambda_2=None): self.model = model self.acquisition_scheme = acquisition_scheme @@ -56,64 +56,9 @@ def __init__(self, acquisition_scheme, model, x0_vector=None, len(lambda_2) != self.model.N_models: raise ValueError("Number of regularization weights should" "correspond to the number of compartments!") - self.lambda_1 = lambda_1 self.lambda_2 = lambda_2 - self.Nt = Nt - - # Make a list of parameters that are not fixed and that require - # tessellation of parameter ranges - dir_params = [p for p in self.model.parameter_names - if p.endswith('mu')] - grid_params =\ - [p for p in self.model.parameter_names - if not p.endswith('mu') and not p.startswith('partial_volume')] - - # Compute length of the vector x0 - x0_len = 0 - for m_idx in range(self.model.N_models): - m_atoms = 1 - for p in self.model.models[m_idx].parameter_names: - if self.model.model_names[m_idx] + p in grid_params: - m_atoms *= Nt - x0_len += m_atoms - - # Creating parameter tessellation grids and corresponding indices - # TODO: move the matrix/grid/indices definition to the - # modeling_framework module. - self.grids, self.idx = {}, {} - for m_idx in range(self.model.N_models): - model = self.model.models[m_idx] - model_name = self.model.model_names[m_idx] - - param_sampling, grid_params_names = [], [] - m_atoms = 1 - for p in model.parameter_names: - if model_name + p not in grid_params: - continue - grid_params_names.append(model_name + p) - p_range = self.model.parameter_ranges[model_name + p] - self.grids[model_name + p] = np.full(x0_len, np.mean(p_range)) - param_sampling.append(np.linspace(p_range[0], p_range[1], - self.Nt, endpoint=True)) - m_atoms *= self.Nt - - self.idx[model_name] =\ - sum([len(self.idx[k]) for k in self.idx]) + np.arange(m_atoms) - params_mesh = np.meshgrid(*param_sampling) - for p_idx, p in enumerate(grid_params_names): - self.grids[p][self.idx[model_name]] =\ - np.ravel(params_mesh[p_idx]) - - self.grids['partial_volume_' + str(m_idx)] = np.zeros(x0_len) - self.grids['partial_volume_' + - str(m_idx)][self.idx[model_name]] = 1. - - self.grids[dir_params[0]] = [0, 0] - self.M = self.model.simulate_signal(acquisition_scheme, self.grids) - self.M = self.M[:, ~acquisition_scheme.b0_mask].T - def __call__(self, data): """ The fitting function of AMICO optimizer. From 873fb4ef0b02acbc27f0cce3c97aa041bdd29960 Mon Sep 17 00:00:00 2001 From: Sara04 Date: Tue, 9 Jun 2020 16:11:50 +0200 Subject: [PATCH 09/28] Remove questions; fix typos --- dmipy/core/modeling_framework.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 981540c0..2061bf44 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2234,9 +2234,6 @@ def __init__(self, models, S0_tissue_responses=None, Nt=10, self._prepare_model_properties() self._check_for_double_model_class_instances() - # QUESTION: I think we should create multi-compartment model - # somewhere in __init__ to avoid creating it every time we create - # forward model matrix self.mc_model = MultiCompartmentModel(models=self.models) self._prepare_parameters_to_optimize() @@ -2341,7 +2338,7 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs, **kwargs): np.arange(m_atoms) params_mesh = np.meshgrid(*param_sampling) for p_idx, p in enumerate(grid_params_names): - self.grids[p][self._amico_idx[model_name]] =\ + self._amico_grid[p][self._amico_idx[model_name]] =\ np.ravel(params_mesh[p_idx]) self._amico_grid['partial_volume_' + str(m_idx)] = \ @@ -2351,12 +2348,8 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs, **kwargs): self._freezed_parameters_vector = True for d_idx, dp in enumerate(dir_params): - self._amico_grids[dp] = model_dirs[d_idx] + self._amico_grid[dp] = model_dirs[d_idx] - # 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 return self.mc_model.simulate_signal(acquisition_scheme, self._amico_grid).T From 02ad74e061049b4e7f4b4fec6408e2c472512791 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Tue, 18 Aug 2020 16:57:35 +0200 Subject: [PATCH 10/28] setup AmicoCvxpyOptimizer in MultiCompartmentAMICOModel 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. --- dmipy/core/modeling_framework.py | 102 ++++++++++++++++++------------- 1 file changed, 58 insertions(+), 44 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 25c43aed..750c3b8d 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -17,9 +17,9 @@ FittedMultiCompartmentSphericalMeanModel, FittedMultiCompartmentSphericalHarmonicsModel, FittedMultiCompartmentAMICOModel) +from ..optimizers.amico_cvxpy import AmicoCvxpyOptimizer from ..optimizers.brute2fine import ( GlobalBruteOptimizer, Brute2FineOptimizer) -from ..optimizers.amico_cvxpy import AmicoCvxpyOptimizer from ..optimizers.mix import MixOptimizer from ..optimizers.multi_tissue_convex_optimizer import ( MultiTissueConvexOptimizer) @@ -131,14 +131,14 @@ def parameter_vector_to_parameters(self, parameter_vector): if parameter_vector.ndim == 1: for parameter, card in self.parameter_cardinality.items(): parameters[parameter] = parameter_vector[ - current_pos: current_pos + card] + current_pos: current_pos + card] if card == 1: parameters[parameter] = parameters[parameter][0] current_pos += card else: for parameter, card in self.parameter_cardinality.items(): parameters[parameter] = parameter_vector[ - ..., current_pos: current_pos + card] + ..., current_pos: current_pos + card] if card == 1: parameters[parameter] = parameters[parameter][..., 0] current_pos += card @@ -1261,7 +1261,7 @@ def fit(self, acquisition_scheme, data, fitted_parameters = np.zeros_like(x0_, dtype=float) fitted_parameters[mask_pos] = ( - fitted_parameters_lin * self.scales_for_optimization) + fitted_parameters_lin * self.scales_for_optimization) return FittedMultiCompartmentModel(self, S0, mask, fitted_parameters, fitted_mt_fractions) @@ -1657,7 +1657,7 @@ def fit(self, acquisition_scheme, data, fitted_parameters = np.zeros_like(x0_, dtype=float) fitted_parameters[mask_pos] = ( - fitted_parameters_lin * self.scales_for_optimization) + fitted_parameters_lin * self.scales_for_optimization) return FittedMultiCompartmentSphericalMeanModel( self, S0, mask, fitted_parameters, fitted_mt_fractions) @@ -2199,7 +2199,7 @@ def __call__(self, acquisition_scheme, **kwargs): 'sh_coeff'] for i, name in enumerate(self.partial_volume_names): sh_coeff[self.optimizer.vf_indices[i]] = ( - kwargs[name] / (2 * np.sqrt(np.pi))) + kwargs[name] / (2 * np.sqrt(np.pi))) E = np.dot(A, sh_coeff) return E @@ -2223,13 +2223,14 @@ class MultiCompartmentAMICOModel(MultiCompartmentModelProperties): argument list), deprecated, for testing only. """ + def __init__(self, models, S0_tissue_responses=None, Nt=10, parameter_links=None): self.models = models self.N_models = len(models) if S0_tissue_responses is not None: if len(S0_tissue_responses) != self.N_models: - msg = 'Number of S0_tissue responses {} must be same as '\ + msg = 'Number of S0_tissue responses {} must be same as ' \ 'number of input models {}.' raise ValueError( msg.format(len(S0_tissue_responses), self.N_models)) @@ -2290,24 +2291,25 @@ def amico_idx(self): matrix.""" return self._amico_idx - def forward_model_matrix(self, acquisition_scheme, model_dirs, **kwargs): + def forward_model_matrix(self, acquisition_scheme, model_dirs=None, + **kwargs): """Creates forward model matrix, including parameter tessellation grid and corresponding indices. - """ - """ + Arguments: acquisition_scheme: instance containing acquisition protocol model_dirs: list containing direction of all models in - multi-compartment model + multi-compartment model. By default it uses the ones fixed + in the multi-compartment model. Default: None. Returns: observation matrix M """ - - dir_params = [p for p in self.mc_model.parameter_names - if p.endswith('mu')] - if len(dir_params) != len(model_dirs): - raise ValueError("Length of model_dirs should correspond " - "to the number of directional parameters!") + if model_dirs is not None: + dir_params = [p for p in self.mc_model.parameter_names + if p.endswith('mu')] + if len(dir_params) != len(model_dirs): + raise ValueError("Length of model_dirs should correspond " + "to the number of directional parameters!") if not self._freezed_parameters_vector: self._amico_grid, self._amico_idx = {}, {} @@ -2343,13 +2345,13 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs, **kwargs): self.Nt, endpoint=True)) m_atoms *= self.Nt - self._amico_idx[model_name] =\ + self._amico_idx[model_name] = \ sum([len(self._amico_idx[k]) for k in self._amico_idx]) + \ np.arange(m_atoms) params_mesh = np.meshgrid(*param_sampling) for p_idx, p in enumerate(grid_params_names): - self._amico_grid[p][self._amico_idx[model_name]] =\ + self._amico_grid[p][self._amico_idx[model_name]] = \ np.ravel(params_mesh[p_idx]) self._amico_grid['partial_volume_' + str(m_idx)] = \ @@ -2358,15 +2360,17 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs, **kwargs): str(m_idx)][self._amico_idx[model_name]] = 1. self._freezed_parameters_vector = True - for d_idx, dp in enumerate(dir_params): - self._amico_grid[dp] = model_dirs[d_idx] + if model_dirs is not None: + for d_idx, dp in enumerate(dir_params): + self._amico_grid[dp] = model_dirs[d_idx] return self.mc_model.simulate_signal(acquisition_scheme, - self._amico_grid).T + self.amico_grid).T def fit(self, acquisition_scheme, data, mask=None, - maxiter=300, + lambda_1=None, + lambda_2=None, use_parallel_processing=have_pathos, number_of_processors=None): """ The main data fitting function of a MultiCompartmentModel. @@ -2403,9 +2407,12 @@ def fit(self, acquisition_scheme, data, or an N-dimensional dataset. mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), Optional mask of voxels to be included in the optimization. - maxiter : integer, - How many iterations are allowed in the optimization process. - Defaults to 300. + lambda_1 : vector with one value per compartment. The + coefficients will be applied to the L1 regularization of each + compartment. Default: 0.0. + lambda_2 : vector with one value per compartment. The + coefficients will be applied to the L2 regularization of each + compartment. Default: 0.0. use_parallel_processing : bool, whether or not to use parallel processing using pathos. number_of_processors : integer, @@ -2469,13 +2476,23 @@ def fit(self, acquisition_scheme, data, np.r_[N_voxels, N_parameters], dtype=float) start = time() - opt = AmicoCvxpyOptimizer(self, self.scheme, x0_) # TODO: fix params - # setup self.amico_grids and self.amico_indices - _ = self.forward_model_matrix() + if lambda_1 is None: + l1 = np.zeros(self.N_models) + else: + l1 = lambda_1 + + if lambda_2 is None: + l2 = np.zeros(self.N_models) + else: + l2 = lambda_2 + + opt = AmicoCvxpyOptimizer(self, self.scheme, x0_, l1, l2) + + M = self.forward_model_matrix(self.scheme) + def fit_func(data): - # TODO: extract the parameter value for each param and return the - # parameter vector in the voxel. - raise NotImplementedError + return opt(data, M, self.amico_grid, self.amico_idx) + print('Setup AMICO optimizer in {} seconds'.format( time() - start)) @@ -2484,9 +2501,7 @@ def fit_func(data): start = time() for idx, pos in enumerate(zip(*mask_pos)): voxel_E = data_[pos] / S0[pos] - voxel_x0_vector = x0_[pos] - # TODO: preprocess the x0 as required by the solver - fit_args = (voxel_E, ) + fit_args = (voxel_E,) if use_parallel_processing: fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) @@ -2526,13 +2541,12 @@ def fit_func(data): fitted_parameters = np.zeros_like(x0_, dtype=float) fitted_parameters[mask_pos] = ( - fitted_parameters_lin * self.scales_for_optimization) + fitted_parameters_lin * self.scales_for_optimization) - # TODO: pass the forward model matrix, the AmicoCvxpyOptimizer object - # and the parameter indices dictionary to the - # FittedMultiCompartmentAMICOModel . return FittedMultiCompartmentAMICOModel( - self, S0, mask, fitted_parameters, fitted_mt_fractions) + self, S0, mask, fitted_parameters, M, self.amico_idx, opt, + fitted_mt_fractions + ) def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): """ @@ -2609,7 +2623,7 @@ def __call__(self, acquisition_scheme_or_vertices, partial_volumes = [1.] for model_name, model, partial_volume in zip( - self.model_names, self.models, partial_volumes + self.model_names, self.models, partial_volumes ): parameters = {} for parameter in model.parameter_ranges: @@ -2622,9 +2636,9 @@ def __call__(self, acquisition_scheme_or_vertices, ) values = ( - values + - partial_volume * model( - acquisition_scheme_or_vertices, **parameters) + values + + partial_volume * model( + acquisition_scheme_or_vertices, **parameters) ) return values From 85994d94b27d2febe10b10869b15a8c6bef6e69a Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Tue, 18 Aug 2020 16:59:02 +0200 Subject: [PATCH 11/28] fix regularization parameter bug in AMICO framework 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. --- dmipy/optimizers/amico_cvxpy.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index 670d7bae..ebdc140c 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -55,6 +55,12 @@ def __init__(self, model, acquisition_scheme, x0_vector=None, len(lambda_2) != self.model.N_models: raise ValueError("Number of regularization weights should" "correspond to the number of compartments!") + if lambda_1 is None: + lambda_1 = np.zeros(self.model.N_models) + + if lambda_1 is None: + lambda_1 = np.zeros(self.model.N_models) + self.lambda_1 = lambda_1 self.lambda_2 = lambda_2 From 55b76fe412e05923abf6181875158ff6b0785b46 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Wed, 19 Aug 2020 18:50:46 +0200 Subject: [PATCH 12/28] fix set_X_parameter in MultiCompartmentAMICOModel Any parameter that is set in an AMICO model, it must be set also in the underlying MultiCompartmentModel. This commit adds this operation. --- dmipy/core/modeling_framework.py | 47 +++++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 13 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 750c3b8d..9d5dd34e 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2246,7 +2246,8 @@ def __init__(self, models, S0_tissue_responses=None, Nt=10, self._prepare_model_properties() self._check_for_double_model_class_instances() - self.mc_model = MultiCompartmentModel(models=self.models) + self.mc_model = MultiCompartmentModel( + models=self.models, S0_tissue_responses=self.S0_tissue_responses) self._prepare_parameters_to_optimize() self._check_for_NMR_and_other_models() @@ -2273,7 +2274,7 @@ def _check_for_NMR_and_other_models(self): raise ValueError(msg) def _check_if_model_orientations_are_fixed(self): - if 'orientation' in self.parameter_types.values: + if 'orientation' in self.parameter_types.values(): msg = 'The orientation parameters must be fixed a priori.' raise ValueError(msg) @@ -2452,15 +2453,6 @@ def fit(self, acquisition_scheme, data, N_parameters = len(self.bounds_for_optimization) N_voxels = np.sum(mask) - # make starting parameters and data the same size - x0_ = self.parameter_initial_guess_to_parameter_vector( - **self.x0_parameters) - x0_ = homogenize_x0_to_data( - data_, x0_) - x0_bool = np.all( - np.isnan(x0_), axis=tuple(np.arange(x0_.ndim - 1))) - x0_[..., ~x0_bool] /= self.scales_for_optimization[~x0_bool] - if use_parallel_processing and not have_pathos: msg = 'Cannot use parallel processing without pathos.' raise ValueError(msg) @@ -2486,7 +2478,7 @@ def fit(self, acquisition_scheme, data, else: l2 = lambda_2 - opt = AmicoCvxpyOptimizer(self, self.scheme, x0_, l1, l2) + opt = AmicoCvxpyOptimizer(self, self.scheme, l1, l2) M = self.forward_model_matrix(self.scheme) @@ -2539,7 +2531,7 @@ def fit_func(data): fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) fitted_mt_fractions[mask_pos] = mt_fractions - fitted_parameters = np.zeros_like(x0_, dtype=float) + fitted_parameters = np.zeros_like(fitted_parameters_lin, dtype=float) fitted_parameters[mask_pos] = ( fitted_parameters_lin * self.scales_for_optimization) @@ -2548,6 +2540,35 @@ def fit_func(data): fitted_mt_fractions ) + def set_equal_parameter(self, parameter_name_in, parameter_name_out): + p = (parameter_name_in, parameter_name_out) + super(MultiCompartmentAMICOModel, self).set_equal_parameter(*p) + self.mc_model.set_equal_parameter(*p) + + def set_fixed_parameter(self, parameter_name, value): + p = (parameter_name, value) + super(MultiCompartmentAMICOModel, self).set_fixed_parameter(*p) + self.mc_model.set_fixed_parameter(*p) + + def set_fractional_parameter(self, parameter1_smaller_equal_than, + parameter2): + p = (parameter1_smaller_equal_than, parameter2) + super(MultiCompartmentAMICOModel, self).set_fractional_parameter(*p) + self.mc_model.set_fractional_parameter(*p) + + def set_tortuous_parameter(self, lambda_perp_parameter_name, + lambda_par_parameter_name, + volume_fraction_intra_parameter_name, + volume_fraction_extra_parameter_name, + S0_correction=False): + p = (lambda_perp_parameter_name, + lambda_par_parameter_name, + volume_fraction_intra_parameter_name, + volume_fraction_extra_parameter_name, + S0_correction) + super(MultiCompartmentAMICOModel, self).set_tortuous_parameter(*p) + self.mc_model.set_tortuous_parameter(*p) + def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): """ Function to simulate diffusion data for a given acquisition_scheme From 99582097f3ec726628044b1239c56c3e04aaab24 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Wed, 19 Aug 2020 18:56:04 +0200 Subject: [PATCH 13/28] fix parameter indices in FittedMultiCompartmentAMICOModel --- dmipy/core/fitted_modeling_framework.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index 3bab720d..12924d15 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -961,7 +961,7 @@ def fitted_parameters(self): fitted = {} for pname in self.model.parameter_names: prange = self.model.parameter_ranges[pname] - idx = self.model.parameter_indices[pname] + idx = self.model.amico_idx[pname] px = self.fitted_parameters_vector[idx] value = np.zeros_like(prange[0]) for v, x in zip(prange, px): # TODO: can be optimized by using dot From 626444f964f75540315c502e7f54abdf82a90a8f Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Wed, 19 Aug 2020 18:56:59 +0200 Subject: [PATCH 14/28] refactor distribution variable in AmicoCvxpyOptimizer --- dmipy/optimizers/amico_cvxpy.py | 53 ++++++++++++++++++--------------- 1 file changed, 29 insertions(+), 24 deletions(-) diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index ebdc140c..46baefee 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -24,9 +24,6 @@ class AmicoCvxpyOptimizer: An acquisition scheme that has been instantiated using Dmipy. model: Dmipy MultiCompartmentModel instance, Contains Dmipy multi-compartment model information. - x0_vector: Array of size (Nx), - Vector containing probability distributions of the parameters that are - being estimated lambda_1: Array of size (Nc) Vector containing L2 regularization weights for each compartment lambda_2: Array of size (Nc) @@ -45,11 +42,11 @@ class AmicoCvxpyOptimizer: modeling language for convex optimization." The Journal of Machine Learning Research 17.1 (2016): 2909-2913. """ - def __init__(self, model, acquisition_scheme, x0_vector=None, + def __init__(self, model, acquisition_scheme, lambda_1=None, lambda_2=None): self.model = model self.acquisition_scheme = acquisition_scheme - self.x0_vector = x0_vector + self._distribution = None if len(lambda_1) != self.model.N_models or \ len(lambda_2) != self.model.N_models: @@ -64,7 +61,11 @@ def __init__(self, model, acquisition_scheme, x0_vector=None, self.lambda_1 = lambda_1 self.lambda_2 = lambda_2 - def __call__(self, data, M, grid, idx, x0_th=1.e-4): + @property + def distribution(self): + return self._distribution + + def __call__(self, data, M, grid, idx, x_th=1.e-4): """ The fitting function of AMICO optimizer. Parameters @@ -79,7 +80,7 @@ def __call__(self, data, M, grid, idx, x0_th=1.e-4): idx: dict Dictionary containing indices that correspond to the parameters to be estimated for each model within multi-compartment model - x0_th: float + x_th: float Threshold for selecting important atoms after solving NNLS with L1 and L2 regularization terms @@ -98,27 +99,31 @@ def __call__(self, data, M, grid, idx, x0_th=1.e-4): # 2. Selecting important atoms by solving NNLS # regularized with L1 and L2 norms - x0 = cvxpy.Variable(len(self.x0_vector)) + x = cvxpy.Variable(shape=(M.shape[1], )) - cost = 0.5 * cvxpy.sum_squares(M * x0 - data) + cost = 0.5 * cvxpy.sum_squares(M @ x - data) for m_idx, model_name in enumerate(self.model.model_names): + # L1 regularization cost += self.lambda_1[m_idx] * \ - cvxpy.norm(x0[idx[model_name]], 1) + cvxpy.norm(x[idx[model_name]], 1) + # L2 regularization cost += 0.5 * self.lambda_2[m_idx] * \ - cvxpy.norm(x0[idx[model_name]], 2) ** 2 - problem = cvxpy.Problem(cvxpy.Minimize(cost), [x0 >= 0]) + cvxpy.norm(x[idx[model_name]], 2) ** 2 + + problem = cvxpy.Problem(cvxpy.Minimize(cost), [x >= 0]) problem.solve() - self.x0_vector = x0.value # 3. Computing distribution vector x0_vector by solving NNLS - x0_idx_i = self.x0_vector > x0_th - x0_i = cvxpy.Variable(sum(x0_idx_i)) - cost = cvxpy.sum_squares(M[:, x0_idx_i] * x0_i - data) - problem = cvxpy.Problem(cvxpy.Minimize(cost), [x0_i >= 0]) + dist = x.value + x_idx_i = dist > x_th + x_i = cvxpy.Variable(sum(x_idx_i)) + cost = cvxpy.sum_squares(M[:, x_idx_i] * x_i - data) + problem = cvxpy.Problem(cvxpy.Minimize(cost), [x_i >= 0]) problem.solve() - self.x0_vector[~x0_idx_i] = 0. - self.x0_vector[x0_idx_i] = x0_i.value - self.x0_vector /= (np.sum(self.x0_vector) + 1.e-8) + dist[~x_idx_i] = 0. + dist[x_idx_i] = x_i.value + dist /= (np.sum(dist) + 1.e-8) + self._distribution = dist # 4. Estimating parameters based using estimated distribution # vector and tessellation grids @@ -126,7 +131,7 @@ def __call__(self, data, M, grid, idx, x0_th=1.e-4): for m_idx, model_name in enumerate(self.model.model_names): m = self.model.models[m_idx] if 'partial_volume_' + str(m_idx) in grid: - p_estim = np.sum(self.x0_vector[idx[model_name]]) + p_estim = np.sum(self.distribution[idx[model_name]]) fitted_parameter_vector.append(p_estim) for p in m.parameter_names: if p.endswith('mu'): @@ -134,8 +139,8 @@ def __call__(self, data, M, grid, idx, x0_th=1.e-4): if model_name + p in grid: p_estim = \ np.sum(grid[model_name + p][idx[model_name]] * - self.x0_vector[idx[model_name]]) /\ - (np.sum(self.x0_vector[idx[model_name]]) + 1.e-8) + self.distribution[idx[model_name]]) / \ + (np.sum(self.distribution[idx[model_name]]) + 1.e-8) fitted_parameter_vector.append(p_estim) - return fitted_parameter_vector + return np.asarray(fitted_parameter_vector) From d1b71009fcf75a9d27777982422e895957226fe2 Mon Sep 17 00:00:00 2001 From: Sara04 Date: Thu, 20 Aug 2020 11:38:29 +0200 Subject: [PATCH 15/28] Fix fitted parameters conversion from array to dictionary --- dmipy/core/fitted_modeling_framework.py | 19 ++++--------------- dmipy/optimizers/amico_cvxpy.py | 7 ++++--- 2 files changed, 8 insertions(+), 18 deletions(-) diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index 12924d15..e0b5cdf7 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -898,6 +898,7 @@ def mean_squared_error(self, data): mse[~self.mask] = 0 return mse + class FittedMultiCompartmentAMICOModel: def __init__(self, model, S0, mask, fitted_parameters_vector, forward_model_matrix, parameter_indices, optimizer, @@ -957,21 +958,9 @@ def fitted_distribution(self): @property def fitted_parameters(self): - "Returns the fitted parameters as a dictionary." - fitted = {} - for pname in self.model.parameter_names: - prange = self.model.parameter_ranges[pname] - idx = self.model.amico_idx[pname] - px = self.fitted_parameters_vector[idx] - value = np.zeros_like(prange[0]) - for v, x in zip(prange, px): # TODO: can be optimized by using dot - value += v * x - value /= np.sum(px, axis=-1) - # By summing over the last axis, we assume the x have been saved - # in a compatible format, namely with one vector per voxel in the - # last direction. - fitted[pname] = value - return fitted + """Returns the fitted parameters as a dictionary.""" + return self.model.parameter_vector_to_parameters( + self.fitted_parameters_vector) @property def fitted_multi_tissue_fractions(self): diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index 46baefee..a7d46c93 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -130,9 +130,6 @@ def __call__(self, data, M, grid, idx, x_th=1.e-4): fitted_parameter_vector = [] for m_idx, model_name in enumerate(self.model.model_names): m = self.model.models[m_idx] - if 'partial_volume_' + str(m_idx) in grid: - p_estim = np.sum(self.distribution[idx[model_name]]) - fitted_parameter_vector.append(p_estim) for p in m.parameter_names: if p.endswith('mu'): continue @@ -143,4 +140,8 @@ def __call__(self, data, M, grid, idx, x_th=1.e-4): (np.sum(self.distribution[idx[model_name]]) + 1.e-8) fitted_parameter_vector.append(p_estim) + if 'partial_volume_' + str(m_idx) in grid: + p_estim = np.sum(self.distribution[idx[model_name]]) + fitted_parameter_vector.append(p_estim) + return np.asarray(fitted_parameter_vector) From 4e1a38d2e2f811266d15a2ed2d49377242cbaca9 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 20 Aug 2020 17:18:38 +0200 Subject: [PATCH 16/28] remove forward model matrix from FittedMultiCompartmentAMICOModel 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. --- dmipy/core/fitted_modeling_framework.py | 15 +++------------ examples/example_amico_test.ipynb | 6 +++--- 2 files changed, 6 insertions(+), 15 deletions(-) diff --git a/dmipy/core/fitted_modeling_framework.py b/dmipy/core/fitted_modeling_framework.py index e0b5cdf7..0b857e75 100644 --- a/dmipy/core/fitted_modeling_framework.py +++ b/dmipy/core/fitted_modeling_framework.py @@ -901,7 +901,7 @@ def mean_squared_error(self, data): class FittedMultiCompartmentAMICOModel: def __init__(self, model, S0, mask, fitted_parameters_vector, - forward_model_matrix, parameter_indices, optimizer, + parameter_indices, optimizer, fitted_multi_tissue_fractions_vector=None): """ The FittedMultiCompartmentModel instance contains information about the @@ -918,8 +918,6 @@ def __init__(self, model, S0, mask, fitted_parameters_vector, if there are multiple TEs. mask : array of size (N_data,), boolean mask of voxels that were fitted. - forward_model_matrix : array of size N_DWIs-by-Nparameters, - forward model used in the fitting process. parameter_indices : dictionary, keys are parameter names and values are the column indices corresponding to the parameter. @@ -935,15 +933,9 @@ def __init__(self, model, S0, mask, fitted_parameters_vector, self.fitted_parameters_vector = fitted_parameters_vector self.fitted_multi_tissue_fractions_vector = ( fitted_multi_tissue_fractions_vector) - self._forward_model_matrix = forward_model_matrix self._parameter_indices = parameter_indices self._optimizer = optimizer - @property - def forward_model_matrix(self): - """Return forward model matrix.""" - return self._forward_model_matrix - @property def parameter_indices(self): """Return parameter indices.""" @@ -1104,9 +1096,8 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None): E_shape = mask.shape + (acquisition_scheme.N_dwi,) E = np.zeros(E_shape) mask_pos = np.where(mask) - for pos in zip(*mask_pos): - parameter_array = self.fitted_parameters_vector[pos] - E[pos] = self.forward_model_matrix.dot(parameter_array) * S0[pos] + # TODO: get rid of the forward model matrix here. We should use the + # predict function of the underlying multi compartment model return E def R2_coefficient_of_determination(self, data): diff --git a/examples/example_amico_test.ipynb b/examples/example_amico_test.ipynb index 0c7a146c..658f0dd4 100644 --- a/examples/example_amico_test.ipynb +++ b/examples/example_amico_test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -12,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -160,7 +160,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.7.6" } }, "nbformat": 4, From 08d74e1a748390d4a8e689909a93b58dc4f8ebaf Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Fri, 21 Aug 2020 12:21:39 +0200 Subject: [PATCH 17/28] add directions to input of forward_model_matrix in AMICO models 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 --- dmipy/core/modeling_framework.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 9d5dd34e..e286e63c 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2274,9 +2274,11 @@ def _check_for_NMR_and_other_models(self): raise ValueError(msg) def _check_if_model_orientations_are_fixed(self): - if 'orientation' in self.parameter_types.values(): - msg = 'The orientation parameters must be fixed a priori.' - raise ValueError(msg) + msg = lambda p: 'Parameter {} must be fixed a priori.'.format(p) + for k, v in self.parameter_types.items(): + if v == 'orientation': + if self.parameter_optimization_flags[k]: + raise ValueError(msg) @property def amico_grid(self): @@ -2430,6 +2432,8 @@ def fit(self, acquisition_scheme, data, self._check_model_params_with_acquisition_params(acquisition_scheme) self._check_acquisition_scheme_has_b0s(acquisition_scheme) self._check_if_model_orientations_are_fixed() + # TODO: add check if some non-directional parameter is different for + # each voxel # estimate S0 self.scheme = acquisition_scheme @@ -2480,9 +2484,13 @@ def fit(self, acquisition_scheme, data, opt = AmicoCvxpyOptimizer(self, self.scheme, l1, l2) - M = self.forward_model_matrix(self.scheme) + dir_par_names = [] + for k, v in self.parameter_types.items(): + if v == 'orientation': + dir_par_names.append(k) - def fit_func(data): + def fit_func(data, directions): + M = self.forward_model_matrix(self.scheme, directions) return opt(data, M, self.amico_grid, self.amico_idx) print('Setup AMICO optimizer in {} seconds'.format( @@ -2493,7 +2501,10 @@ def fit_func(data): start = time() for idx, pos in enumerate(zip(*mask_pos)): voxel_E = data_[pos] / S0[pos] - fit_args = (voxel_E,) + directions = [] + for p in dir_par_names: + directions.append(self.x0_parameters[p][pos]) + fit_args = (voxel_E, directions) if use_parallel_processing: fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) @@ -2531,12 +2542,15 @@ def fit_func(data): fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) fitted_mt_fractions[mask_pos] = mt_fractions - fitted_parameters = np.zeros_like(fitted_parameters_lin, dtype=float) + shape = () + # TODO: the shape must be the one of the volume (data_) plus the + # number of parameters (including the directional ones) + fitted_parameters = np.zeros(shape=shape, dtype=float) fitted_parameters[mask_pos] = ( fitted_parameters_lin * self.scales_for_optimization) return FittedMultiCompartmentAMICOModel( - self, S0, mask, fitted_parameters, M, self.amico_idx, opt, + self, S0, mask, fitted_parameters, self.amico_idx, opt, fitted_mt_fractions ) From 735b2a172f746c0580c51407f36d3398f6373c7d Mon Sep 17 00:00:00 2001 From: Sara04 Date: Sat, 22 Aug 2020 00:17:56 +0200 Subject: [PATCH 18/28] Check if non-directional fixed parameters are unique for all voxels; add direction to the array of fitted parameters; update shape of fitted parameters --- dmipy/core/modeling_framework.py | 37 +++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index e286e63c..c1194ff1 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2309,7 +2309,8 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, """ if model_dirs is not None: dir_params = [p for p in self.mc_model.parameter_names - if p.endswith('mu')] + if self.mc_model.parameter_types[p] == + 'orientation'] if len(dir_params) != len(model_dirs): raise ValueError("Length of model_dirs should correspond " "to the number of directional parameters!") @@ -2319,7 +2320,7 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, grid_params = \ [p for p in self.mc_model.parameter_names - if not p.endswith('mu') and + if self.mc_model.parameter_types[p] == 'normal' and not p.startswith('partial_volume')] # Compute length of the vector x0 @@ -2432,8 +2433,6 @@ def fit(self, acquisition_scheme, data, self._check_model_params_with_acquisition_params(acquisition_scheme) self._check_acquisition_scheme_has_b0s(acquisition_scheme) self._check_if_model_orientations_are_fixed() - # TODO: add check if some non-directional parameter is different for - # each voxel # estimate S0 self.scheme = acquisition_scheme @@ -2491,7 +2490,8 @@ def fit(self, acquisition_scheme, data, def fit_func(data, directions): M = self.forward_model_matrix(self.scheme, directions) - return opt(data, M, self.amico_grid, self.amico_idx) + return np.concatenate((np.ravel(directions), + opt(data, M, self.amico_grid, self.amico_idx))) print('Setup AMICO optimizer in {} seconds'.format( time() - start)) @@ -2501,15 +2501,13 @@ def fit_func(data, directions): start = time() for idx, pos in enumerate(zip(*mask_pos)): voxel_E = data_[pos] / S0[pos] - directions = [] - for p in dir_par_names: - directions.append(self.x0_parameters[p][pos]) + directions = [self.x0_parameters[p][pos] for p in dir_par_names] fit_args = (voxel_E, directions) - if use_parallel_processing: fitted_parameters_lin[idx] = pool.apipe(fit_func, *fit_args) else: fitted_parameters_lin[idx] = fit_func(*fit_args) + if use_parallel_processing: fitted_parameters_lin = np.array( [p.get() for p in fitted_parameters_lin]) @@ -2542,10 +2540,10 @@ def fit_func(data, directions): fitted_mt_fractions = np.zeros(np.r_[mask.shape, self.N_models]) fitted_mt_fractions[mask_pos] = mt_fractions - shape = () - # TODO: the shape must be the one of the volume (data_) plus the - # number of parameters (including the directional ones) - fitted_parameters = np.zeros(shape=shape, dtype=float) + p_shape = tuple(list(data_.shape[:-1]) + + [fitted_parameters_lin.shape[-1]]) + + fitted_parameters = np.zeros(shape=p_shape, dtype=float) fitted_parameters[mask_pos] = ( fitted_parameters_lin * self.scales_for_optimization) @@ -2560,6 +2558,19 @@ def set_equal_parameter(self, parameter_name_in, parameter_name_out): self.mc_model.set_equal_parameter(*p) def set_fixed_parameter(self, parameter_name, value): + + if self.parameter_types[parameter_name] == 'normal': + if isinstance(value, list) or isinstance(value, np.ndarray): + if len(value) == 1: + value = value[0] + else: + raise ValueError("Parameter {} should be " + "unique for all voxels.". + format(parameter_name)) + if not isinstance(value, float): + raise ValueError("Parameter {} should be float". + format(parameter_name)) + p = (parameter_name, value) super(MultiCompartmentAMICOModel, self).set_fixed_parameter(*p) self.mc_model.set_fixed_parameter(*p) From 99168c03892e4eb3029467c259edb8e07cc30758 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Tue, 25 Aug 2020 16:35:14 +0200 Subject: [PATCH 19/28] normalize dictionary in AmicoCvxpyOptimizer --- dmipy/optimizers/amico_cvxpy.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index a7d46c93..35ae1d88 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -94,6 +94,10 @@ def __call__(self, data, M, grid, idx, x_th=1.e-4): # 1. Contracting matrix M and data to have one b=0 value M = np.vstack((np.mean(M[self.acquisition_scheme.b0_mask, :], axis=0), M[~self.acquisition_scheme.b0_mask, :])) + # normalize the columns of the matrix + norms = np.linalg.norm(M, 1, axis=0) + M /= norms + data = np.append(np.mean(data[self.acquisition_scheme.b0_mask]), data[~self.acquisition_scheme.b0_mask]) @@ -115,6 +119,7 @@ def __call__(self, data, M, grid, idx, x_th=1.e-4): # 3. Computing distribution vector x0_vector by solving NNLS dist = x.value + dist /= norms # rescale by the original norm of the columns x_idx_i = dist > x_th x_i = cvxpy.Variable(sum(x_idx_i)) cost = cvxpy.sum_squares(M[:, x_idx_i] * x_i - data) From fe3f35c4fe890e5abd8d178d5a89ea36a22cd911 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Tue, 25 Aug 2020 16:35:54 +0200 Subject: [PATCH 20/28] fix unexpected behaviour of set_fixed_parameter of AMICO models --- dmipy/core/modeling_framework.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index c1194ff1..6f6a0365 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2559,17 +2559,19 @@ def set_equal_parameter(self, parameter_name_in, parameter_name_out): def set_fixed_parameter(self, parameter_name, value): + message = "Parameter '{}' should be unique for all voxels.".format( + parameter_name) if self.parameter_types[parameter_name] == 'normal': - if isinstance(value, list) or isinstance(value, np.ndarray): - if len(value) == 1: - value = value[0] + if isinstance(value, list): + value = np.asarray(value) + + if isinstance(value, np.ndarray): + if value.size == 1: + value = value.ravel()[0] else: - raise ValueError("Parameter {} should be " - "unique for all voxels.". - format(parameter_name)) - if not isinstance(value, float): - raise ValueError("Parameter {} should be float". - format(parameter_name)) + raise ValueError(message) + + value = float(value) p = (parameter_name, value) super(MultiCompartmentAMICOModel, self).set_fixed_parameter(*p) From 6d6565814cee8539ad4975bf8dc736d4ac710c4f Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Tue, 25 Aug 2020 16:48:04 +0200 Subject: [PATCH 21/28] change dictionary normalization of AmicoCvxpyOptimizer from l1 to l2 --- dmipy/optimizers/amico_cvxpy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index 35ae1d88..f02ccc31 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -95,7 +95,7 @@ def __call__(self, data, M, grid, idx, x_th=1.e-4): M = np.vstack((np.mean(M[self.acquisition_scheme.b0_mask, :], axis=0), M[~self.acquisition_scheme.b0_mask, :])) # normalize the columns of the matrix - norms = np.linalg.norm(M, 1, axis=0) + norms = np.linalg.norm(M, 2, axis=0) M /= norms data = np.append(np.mean(data[self.acquisition_scheme.b0_mask]), From 51aa0abd54fe2757aa9f7cabc80ae29d88507da5 Mon Sep 17 00:00:00 2001 From: Sara04 Date: Thu, 27 Aug 2020 15:45:35 +0200 Subject: [PATCH 22/28] Check if partial volume parameters are set, check if there are parameters to be estimated --- dmipy/core/modeling_framework.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 6f6a0365..738fc958 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2358,16 +2358,21 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, self._amico_grid[p][self._amico_idx[model_name]] = \ np.ravel(params_mesh[p_idx]) - self._amico_grid['partial_volume_' + str(m_idx)] = \ - np.zeros(x0_len) - self._amico_grid['partial_volume_' + - str(m_idx)][self._amico_idx[model_name]] = 1. + if 'partial_volume_' + str(m_idx) in self.mc_model.parameter_names: + self._amico_grid['partial_volume_' + str(m_idx)] = \ + np.zeros(x0_len) + self._amico_grid['partial_volume_' + + str(m_idx)][self._amico_idx[model_name]] = 1. self._freezed_parameters_vector = True if model_dirs is not None: for d_idx, dp in enumerate(dir_params): self._amico_grid[dp] = model_dirs[d_idx] + if 'normal' not in self.mc_model.parameter_types.values(): + msg = 'No parameters to be estimated, all parameters are set.' + raise ValueError(msg) + return self.mc_model.simulate_signal(acquisition_scheme, self.amico_grid).T From 692bb1514fc28d16bff47098c7fb877a2b38261c Mon Sep 17 00:00:00 2001 From: Sara04 Date: Thu, 27 Aug 2020 20:12:21 +0200 Subject: [PATCH 23/28] Add scaling of parameters --- dmipy/core/modeling_framework.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 738fc958..d3686511 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2343,9 +2343,12 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, continue grid_params_names.append(model_name + p) p_range = self.mc_model.parameter_ranges[model_name + p] + p_scale = self.mc_model.parameter_scales[model_name + p] + self._amico_grid[model_name + p] = \ - np.full(x0_len, np.mean(p_range)) - param_sampling.append(np.linspace(p_range[0], p_range[1], + np.full(x0_len, np.mean(p_range) * p_scale) + param_sampling.append(np.linspace(p_range[0] * p_scale, + p_range[1] * p_scale, self.Nt, endpoint=True)) m_atoms *= self.Nt From 3cee2c8940afc233ead2beec04c94b995a9164dd Mon Sep 17 00:00:00 2001 From: Sara04 Date: Thu, 27 Aug 2020 21:23:29 +0200 Subject: [PATCH 24/28] Add possibility to estimate circular parameters (as in Bingham distribution); raise ValueError if number of atoms in forward model matrix is greater than max_atoms --- dmipy/core/modeling_framework.py | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index d3686511..19edd0d1 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2224,7 +2224,7 @@ class MultiCompartmentAMICOModel(MultiCompartmentModelProperties): deprecated, for testing only. """ - def __init__(self, models, S0_tissue_responses=None, Nt=10, + def __init__(self, models, S0_tissue_responses=None, Nt=10, max_atoms=2000, parameter_links=None): self.models = models self.N_models = len(models) @@ -2236,6 +2236,7 @@ def __init__(self, models, S0_tissue_responses=None, Nt=10, msg.format(len(S0_tissue_responses), self.N_models)) self.S0_tissue_responses = S0_tissue_responses self.Nt = Nt + self.max_atoms = max_atoms self.parameter_links = parameter_links if parameter_links is None: self.parameter_links = [] @@ -2320,7 +2321,7 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, grid_params = \ [p for p in self.mc_model.parameter_names - if self.mc_model.parameter_types[p] == 'normal' and + if self.mc_model.parameter_types[p] in ['normal', 'circular'] and not p.startswith('partial_volume')] # Compute length of the vector x0 @@ -2361,18 +2362,29 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, self._amico_grid[p][self._amico_idx[model_name]] = \ np.ravel(params_mesh[p_idx]) - if 'partial_volume_' + str(m_idx) in self.mc_model.parameter_names: - self._amico_grid['partial_volume_' + str(m_idx)] = \ - np.zeros(x0_len) - self._amico_grid['partial_volume_' + - str(m_idx)][self._amico_idx[model_name]] = 1. + if 'partial_volume_' + str(m_idx) in \ + self.mc_model.parameter_names: + p_volume_name = 'partial_volume_' + str(m_idx) + self._amico_grid[p_volume_name] = np.zeros(x0_len) + self._amico_grid[p_volume_name][ + self._amico_idx[model_name]] = 1. self._freezed_parameters_vector = True + if m_atoms > self.max_atoms: + raise ValueError("Large number of unknown parameters {} " + "resulted in large number of atoms in " + "forward model matrix.". + format(grid_params_names), + "Size of the forward model matrix is [{}, {}]". + format(acquisition_scheme.number_of_measurements, + m_atoms)) + if model_dirs is not None: for d_idx, dp in enumerate(dir_params): self._amico_grid[dp] = model_dirs[d_idx] - if 'normal' not in self.mc_model.parameter_types.values(): + if 'normal' not in self.mc_model.parameter_types.values() and \ + 'circular' not in self.mc_model.parameter_types.values(): msg = 'No parameters to be estimated, all parameters are set.' raise ValueError(msg) From 307958c53197eb89e964eec06bbfa2aeeb459bec Mon Sep 17 00:00:00 2001 From: Sara04 Date: Thu, 27 Aug 2020 21:51:40 +0200 Subject: [PATCH 25/28] Fix counting of total number of atoms; increase max number of atoms --- dmipy/core/modeling_framework.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 19edd0d1..444998a1 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2224,7 +2224,7 @@ class MultiCompartmentAMICOModel(MultiCompartmentModelProperties): deprecated, for testing only. """ - def __init__(self, models, S0_tissue_responses=None, Nt=10, max_atoms=2000, + def __init__(self, models, S0_tissue_responses=None, Nt=10, max_atoms=20000, parameter_links=None): self.models = models self.N_models = len(models) @@ -2370,14 +2370,15 @@ def forward_model_matrix(self, acquisition_scheme, model_dirs=None, self._amico_idx[model_name]] = 1. self._freezed_parameters_vector = True - if m_atoms > self.max_atoms: + n_atoms = sum([len(self._amico_idx[m]) for m in self._amico_idx]) + if n_atoms > self.max_atoms: raise ValueError("Large number of unknown parameters {} " "resulted in large number of atoms in " "forward model matrix.". format(grid_params_names), "Size of the forward model matrix is [{}, {}]". format(acquisition_scheme.number_of_measurements, - m_atoms)) + n_atoms)) if model_dirs is not None: for d_idx, dp in enumerate(dir_params): From 0e1f5bfa4f3145c252c5ef766c508aec19be4785 Mon Sep 17 00:00:00 2001 From: Sara04 Date: Thu, 27 Aug 2020 22:41:16 +0200 Subject: [PATCH 26/28] Add tests for forward model matrix for ball and stick model and NODDI Bingham model --- dmipy/core/tests/test_amico_models.py | 185 +++++++++++++++++++++++++- 1 file changed, 180 insertions(+), 5 deletions(-) diff --git a/dmipy/core/tests/test_amico_models.py b/dmipy/core/tests/test_amico_models.py index 6710c547..1ca75dbe 100644 --- a/dmipy/core/tests/test_amico_models.py +++ b/dmipy/core/tests/test_amico_models.py @@ -1,13 +1,17 @@ from dmipy.distributions import distribute_models from dmipy.signal_models import cylinder_models, gaussian_models +from dmipy.distributions.distribute_models import SD2BinghamDistributed from dmipy.core import modeling_framework from dmipy.data.saved_acquisition_schemes import ( wu_minn_hcp_acquisition_scheme) import numpy as np +import numpy.linalg from numpy.testing import assert_, assert_raises +from dmipy.core.acquisition_scheme import acquisition_scheme_from_bvalues +from dmipy.core.modeling_framework import MultiCompartmentAMICOModel, MultiCompartmentModel scheme = wu_minn_hcp_acquisition_scheme() - +from numpy.testing import assert_almost_equal def one_compartment_amico_model(): """ @@ -104,9 +108,179 @@ def two_sticks_one_ball_amico_model(): return model, data -def test_forward_model(): - # TODO: test that the observation matrix is defined as expected - pass +def test_forward_model_matrix(): + """Provides tests for creation of forward model matrices.""" + + # Define downsampled acquisition scheme + nM = 60 + scheme_ds = \ + acquisition_scheme_from_bvalues(scheme.bvalues[0:nM], + scheme.gradient_directions[0:nM, :], + scheme.delta[0:nM], + scheme.Delta[0:nM], + scheme.min_b_shell_distance, + scheme.b0_threshold) + # ________________________________________________________________________ + # 1. Ball and stick model + print("Testing forward model matrix for ball and stick model with " + "{} and {} unknown...\n\n".format('G1Ball_1_lambda_iso', + 'C1Stick_1_lambda_par')) + ball = gaussian_models.G1Ball() + stick = cylinder_models.C1Stick() + amico_mc = MultiCompartmentAMICOModel(models=[ball, stick]) + M = amico_mc.forward_model_matrix(scheme_ds, model_dirs=[[0, 0]]) + + mc = MultiCompartmentModel(models=[ball, stick]) + p_range = mc.parameter_ranges['G1Ball_1_lambda_iso'] + p_scale = mc.parameter_scales['G1Ball_1_lambda_iso'] + G1Ball_1_lambda_iso_grid = \ + np.linspace(p_range[0], p_range[1], amico_mc.Nt, endpoint=True) * \ + p_scale + G1Ball_1_lambda_iso_mean = np.mean(p_range) * p_scale + + p_range = mc.parameter_ranges['C1Stick_1_lambda_par'] + p_scale = mc.parameter_scales['C1Stick_1_lambda_par'] + C1Stick_1_lambda_par_grid = \ + np.linspace(p_range[0], p_range[1], amico_mc.Nt, endpoint=True) * \ + p_scale + C1Stick_1_lambda_par_mean = np.mean(p_range) * p_scale + + M_test = np.zeros_like(M) + arguments = {} + arguments['G1Ball_1_lambda_iso'] = G1Ball_1_lambda_iso_grid + arguments['C1Stick_1_lambda_par'] = C1Stick_1_lambda_par_mean + arguments['partial_volume_0'] = 1 + arguments['partial_volume_1'] = 0 + arguments['C1Stick_1_mu'] = [0, 0] + M_test[:, 0:amico_mc.Nt] = mc.simulate_signal(scheme_ds, arguments).T + + arguments['G1Ball_1_lambda_iso'] = G1Ball_1_lambda_iso_mean + arguments['C1Stick_1_lambda_par'] = C1Stick_1_lambda_par_grid + arguments['partial_volume_0'] = 0 + arguments['partial_volume_1'] = 1 + arguments['C1Stick_1_mu'] = [0, 0] + M_test[:, amico_mc.Nt:] = mc.simulate_signal(scheme_ds, arguments).T + + assert_almost_equal(M_test, M) + + print("Testing forward model matrix for ball and stick model with " + "{} unknown...\n\n".format('C1Stick_1_lambda_par')) + amico_mc = MultiCompartmentAMICOModel(models=[ball, stick]) + amico_mc.set_fixed_parameter('G1Ball_1_lambda_iso', 3.e-9) + M = amico_mc.forward_model_matrix(scheme_ds, model_dirs=[[0, 0]]) + + mc.set_fixed_parameter('G1Ball_1_lambda_iso', 3.e-9) + M_test = np.zeros_like(M) + + arguments = {} + arguments['C1Stick_1_lambda_par'] = C1Stick_1_lambda_par_mean + arguments['partial_volume_0'] = 1. + arguments['partial_volume_1'] = 0 + arguments['C1Stick_1_mu'] = [0, 0] + M_test[:, 0] = mc.simulate_signal(scheme_ds, arguments).T + + arguments = {} + arguments['C1Stick_1_lambda_par'] = C1Stick_1_lambda_par_grid + arguments['partial_volume_0'] = 0 + arguments['partial_volume_1'] = 1. + arguments['C1Stick_1_mu'] = [0, 0] + M_test[:, 1:] = mc.simulate_signal(scheme_ds, arguments).T + + assert_almost_equal(M_test, M) + # ________________________________________________________________________ + # 3. NODDI Bingham model + print("Testing forward model matrix for NODDI Bingham model with " + "{}, {}, {} and {} unknown...\n\n". + format('SD2BinghamDistributed_1_SD2Bingham_1_psi', + 'SD2BinghamDistributed_1_SD2Bingham_1_odi', + 'SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction', + 'SD2BinghamDistributed_1_partial_volume_0')) + ball = gaussian_models.G1Ball() + stick = cylinder_models.C1Stick() + zeppelin = gaussian_models.G2Zeppelin() + bingham_dispersed_bundle = SD2BinghamDistributed(models=[stick, zeppelin]) + bingham_dispersed_bundle.set_tortuous_parameter('G2Zeppelin_1_lambda_perp', + 'C1Stick_1_lambda_par', + 'partial_volume_0') + bingham_dispersed_bundle.set_equal_parameter('G2Zeppelin_1_lambda_par', + 'C1Stick_1_lambda_par') + bingham_dispersed_bundle.set_fixed_parameter('G2Zeppelin_1_lambda_par', + 1.7e-9) + amico_mc = MultiCompartmentAMICOModel(models=[ball, + bingham_dispersed_bundle]) + amico_mc.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9) + M = amico_mc.forward_model_matrix(scheme_ds, model_dirs=[[0, 0]]) + + mc = MultiCompartmentModel(models=[ball, bingham_dispersed_bundle]) + mc.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9) + + p_range = mc.parameter_ranges['SD2BinghamDistributed_1_SD2Bingham_1_psi'] + p_scale = mc.parameter_scales['SD2BinghamDistributed_1_SD2Bingham_1_psi'] + SD2BinghamDistributed_1_SD2Bingham_1_psi_grid = \ + np.linspace(p_range[0], p_range[1], amico_mc.Nt, endpoint=True) * \ + p_scale + SD2BinghamDistributed_1_SD2Bingham_1_psi_mean = np.mean(p_range) * p_scale + + p_range = mc.parameter_ranges['SD2BinghamDistributed_1_SD2Bingham_1_odi'] + p_scale = mc.parameter_scales['SD2BinghamDistributed_1_SD2Bingham_1_odi'] + SD2BinghamDistributed_1_SD2Bingham_1_odi_grid = \ + np.linspace(p_range[0], p_range[1], amico_mc.Nt, endpoint=True) * \ + p_scale + SD2BinghamDistributed_1_SD2Bingham_1_odi_mean = np.mean(p_range) * p_scale + + p_range = mc.parameter_ranges[ + 'SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction'] + p_scale = mc.parameter_scales[ + 'SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction'] + SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction_grid = \ + np.linspace(p_range[0], p_range[1], amico_mc.Nt, endpoint=True) * \ + p_scale + SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction_mean = \ + np.mean(p_range) * p_scale + + p_range = mc.parameter_ranges['SD2BinghamDistributed_1_partial_volume_0'] + p_scale = mc.parameter_scales['SD2BinghamDistributed_1_partial_volume_0'] + SD2BinghamDistributed_1_partial_volume_0_grid = \ + np.linspace(p_range[0], p_range[1], amico_mc.Nt, endpoint=True) * \ + p_scale + SD2BinghamDistributed_1_partial_volume_0_mean = np.mean(p_range) * p_scale + + M_test = np.zeros_like(M) + arguments = {} + arguments['SD2BinghamDistributed_1_SD2Bingham_1_psi'] = \ + SD2BinghamDistributed_1_SD2Bingham_1_psi_mean + arguments['SD2BinghamDistributed_1_SD2Bingham_1_odi'] = \ + SD2BinghamDistributed_1_SD2Bingham_1_odi_mean + arguments['SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction'] = \ + SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction_mean + arguments['SD2BinghamDistributed_1_partial_volume_0'] = \ + SD2BinghamDistributed_1_partial_volume_0_mean + arguments['partial_volume_0'] = 1. + arguments['partial_volume_1'] = 0 + arguments['SD2BinghamDistributed_1_SD2Bingham_1_mu'] = [0, 0] + M_test[:, 0] = mc.simulate_signal(scheme_ds, arguments).T + + params_mesh = np.meshgrid(*[ + SD2BinghamDistributed_1_SD2Bingham_1_psi_grid, + SD2BinghamDistributed_1_SD2Bingham_1_odi_grid, + SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction_grid, + SD2BinghamDistributed_1_partial_volume_0_grid]) + + arguments['partial_volume_0'] = 0. + arguments['partial_volume_1'] = 1. + + arguments['SD2BinghamDistributed_1_SD2Bingham_1_psi'] = \ + np.ravel(params_mesh[0]) + arguments['SD2BinghamDistributed_1_SD2Bingham_1_odi'] = \ + np.ravel(params_mesh[1]) + arguments['SD2BinghamDistributed_1_SD2Bingham_1_beta_fraction'] = \ + np.ravel(params_mesh[2]) + arguments['SD2BinghamDistributed_1_partial_volume_0'] = \ + np.ravel(params_mesh[3]) + + M_test[:, 1:] = mc.simulate_signal(scheme_ds, arguments).T + assert_almost_equal(M_test, M) + # ________________________________________________________________________ def test_nnls(): @@ -134,6 +308,7 @@ def test_fitted_amico_model_properties(): # expected properties pass + def test_signal_fitting(): # TODO: test the accuracy of the fitting on a synthetic signal - pass \ No newline at end of file + pass From 70d56cfec1e192067ba7b03e89c499baef9499d1 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Tue, 8 Sep 2020 14:37:38 +0200 Subject: [PATCH 27/28] fix initial condition error test amico --- dmipy/optimizers/tests/test_amico.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/dmipy/optimizers/tests/test_amico.py b/dmipy/optimizers/tests/test_amico.py index 113ade0c..baa22513 100644 --- a/dmipy/optimizers/tests/test_amico.py +++ b/dmipy/optimizers/tests/test_amico.py @@ -181,10 +181,8 @@ def test_amico(lambda_1=[0, 0.0001], lambda_2=[0, 0.0000001], Nt=12): if mc_model.model_names[m_idx] + p in grid_params: m_atoms *= Nt x0_len += m_atoms - x0_vector = np.zeros(x0_len) amico_opt = amico_cvxpy.AmicoCvxpyOptimizer(mc_model, scheme_hcp, - x0_vector=x0_vector, lambda_1=lambda_1, lambda_2=lambda_2) From e7a506bd9f7041d40ddf968dd925fc76b5cde778 Mon Sep 17 00:00:00 2001 From: FRIGO Matteo Date: Thu, 15 Oct 2020 17:30:02 +0200 Subject: [PATCH 28/28] FIX: forbid setting tortuosity from MultiCompartmentAMICOModel 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. --- dmipy/core/modeling_framework.py | 10 +++------- dmipy/optimizers/amico_cvxpy.py | 16 ++++++---------- 2 files changed, 9 insertions(+), 17 deletions(-) diff --git a/dmipy/core/modeling_framework.py b/dmipy/core/modeling_framework.py index 444998a1..b033c535 100644 --- a/dmipy/core/modeling_framework.py +++ b/dmipy/core/modeling_framework.py @@ -2609,13 +2609,9 @@ def set_tortuous_parameter(self, lambda_perp_parameter_name, volume_fraction_intra_parameter_name, volume_fraction_extra_parameter_name, S0_correction=False): - p = (lambda_perp_parameter_name, - lambda_par_parameter_name, - volume_fraction_intra_parameter_name, - volume_fraction_extra_parameter_name, - S0_correction) - super(MultiCompartmentAMICOModel, self).set_tortuous_parameter(*p) - self.mc_model.set_tortuous_parameter(*p) + raise NotImplementedError('The tortuousity constraint must be set ' + 'from a distributed model, not from a ' + 'multi compartment model.') def simulate_signal(self, acquisition_scheme, parameters_array_or_dict): """ diff --git a/dmipy/optimizers/amico_cvxpy.py b/dmipy/optimizers/amico_cvxpy.py index f02ccc31..e4a5280d 100644 --- a/dmipy/optimizers/amico_cvxpy.py +++ b/dmipy/optimizers/amico_cvxpy.py @@ -1,6 +1,5 @@ -import numpy as np import cvxpy - +import numpy as np __all__ = [ 'AmicoCvxpyOptimizer' @@ -42,6 +41,7 @@ class AmicoCvxpyOptimizer: modeling language for convex optimization." The Journal of Machine Learning Research 17.1 (2016): 2909-2913. """ + def __init__(self, model, acquisition_scheme, lambda_1=None, lambda_2=None): self.model = model @@ -93,33 +93,29 @@ def __call__(self, data, M, grid, idx, x_th=1.e-4): # 1. Contracting matrix M and data to have one b=0 value M = np.vstack((np.mean(M[self.acquisition_scheme.b0_mask, :], axis=0), - M[~self.acquisition_scheme.b0_mask, :])) - # normalize the columns of the matrix - norms = np.linalg.norm(M, 2, axis=0) - M /= norms + M[~self.acquisition_scheme.b0_mask, :])) data = np.append(np.mean(data[self.acquisition_scheme.b0_mask]), data[~self.acquisition_scheme.b0_mask]) # 2. Selecting important atoms by solving NNLS # regularized with L1 and L2 norms - x = cvxpy.Variable(shape=(M.shape[1], )) + x = cvxpy.Variable(shape=(M.shape[1],)) cost = 0.5 * cvxpy.sum_squares(M @ x - data) for m_idx, model_name in enumerate(self.model.model_names): # L1 regularization cost += self.lambda_1[m_idx] * \ - cvxpy.norm(x[idx[model_name]], 1) + cvxpy.norm(x[idx[model_name]], 1) # L2 regularization cost += 0.5 * self.lambda_2[m_idx] * \ - cvxpy.norm(x[idx[model_name]], 2) ** 2 + cvxpy.norm(x[idx[model_name]], 2) ** 2 problem = cvxpy.Problem(cvxpy.Minimize(cost), [x >= 0]) problem.solve() # 3. Computing distribution vector x0_vector by solving NNLS dist = x.value - dist /= norms # rescale by the original norm of the columns x_idx_i = dist > x_th x_i = cvxpy.Variable(sum(x_idx_i)) cost = cvxpy.sum_squares(M[:, x_idx_i] * x_i - data)