core

Module: core.acquisition_scheme

dmipy.core.acquisition_scheme.get_sh_order_from_bval(bval)

Estimates minimum sh_order to represent data of given b-value.

class dmipy.core.acquisition_scheme.DmipyAcquisitionScheme(bvalues, gradient_directions, qvalues, gradient_strengths, delta, Delta, TE, min_b_shell_distance, b0_threshold)

Class that calculates and contains all information needed to simulate and fit data using microstructure models.

Attributes

Methods

print_acquisition_info

prints a small summary of the acquisition scheme. Is useful to check if the function correctly separated the shells and if the input parameters were given in the right scale.

to_schemefile(filename)

Exports acquisition scheme information in schemefile format, which can be used by the Camino Monte-Carlo simulator.

Parameters:

filename : string,

location at which to save the schemefile.

visualise_acquisition_G_Delta_rainbow(Delta_start=None, Delta_end=None, G_start=None, G_end=None, bval_isolines=array([ 0, 250, 1000, 2500, 5000, 7500, 10000, 14000]), alpha_shading=0.6)

This function visualizes a q-tau acquisition scheme as a function of gradient strength and pulse separation (big_delta). It represents every measurements at its G and big_delta position regardless of b-vector, with a background of b-value isolines for reference. It assumes there is only one unique pulse length (small_delta) in the acquisition scheme.

Parameters:

Delta_start : float,

optional minimum big_delta that is plotted in seconds

Delta_end : float,

optional maximum big_delta that is plotted in seconds

G_start : float,

optional minimum gradient strength that is plotted in T/m

G_end : float,

optional maximum gradient strength taht is plotted in T/m

bval_isolines : array,

optional array of bvalue isolines that are plotted in background given in s/mm^2

alpha_shading : float between [0-1]

optional shading of the bvalue colors in the background

class dmipy.core.acquisition_scheme.RotationalHarmonicsAcquisitionScheme(dmipy_acquisition_scheme, N_angular_samples=10)

AcquisitionScheme instance that contains the information necessary to calculate the rotational harmonics for a model for every acquisition shell. It is instantiated using a regular DmipyAcquisitionScheme and N_angular_samples determines how many samples are taken between mu=[0., 0.] and mu=[np.pi/2, 0.].

Parameters:

dmipy_acquisition_scheme: DmipyAcquisitionScheme instance :

An acquisition scheme that has been instantiated using dMipy.

N_angular_samples: int :

Integer representing the number of angular samples per shell.

class dmipy.core.acquisition_scheme.SphericalMeanAcquisitionScheme(bvalues, qvalues, gradient_strengths, Deltas, deltas)

Acquisition scheme for isotropic spherical mean models.

dmipy.core.acquisition_scheme.acquisition_scheme_from_bvalues(bvalues, gradient_directions, delta, Delta, TE=None, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Creates an acquisition scheme object from bvalues, gradient directions, pulse duration \(\delta\) and pulse separation time \(\Delta\).

Parameters:

bvalues: 1D numpy array of shape (Ndata) :

bvalues of the acquisition in s/m^2. e.g., a bvalue of 1000 s/mm^2 must be entered as 1000 * 1e6 s/m^2

gradient_directions: 2D numpy array of shape (Ndata, 3) :

gradient directions array of cartesian unit vectors.

delta: float or 1D numpy array of shape (Ndata) :

if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement.

Delta: float or 1D numpy array of shape (Ndata) :

if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement.

min_b_shell_distance : float

minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean.

b0_threshold : float

bvalue threshold for a measurement to be considered a b0 measurement.

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

dmipy.core.acquisition_scheme.acquisition_scheme_from_qvalues(qvalues, gradient_directions, delta, Delta, TE=None, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Creates an acquisition scheme object from qvalues, gradient directions, pulse duration \(\delta\) and pulse separation time \(\Delta\).

Parameters:

qvalues: 1D numpy array of shape (Ndata) :

diffusion sensitization of the acquisition in 1/m. e.g. a qvalue of 10 1/mm must be entered as 10 * 1e3 1/m

gradient_directions: 2D numpy array of shape (Ndata, 3) :

gradient directions array of cartesian unit vectors.

delta: float or 1D numpy array of shape (Ndata) :

if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement.

Delta: float or 1D numpy array of shape (Ndata) :

if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement.

min_b_shell_distance : float

minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean.

b0_threshold : float

bvalue threshold for a measurement to be considered a b0 measurement.

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

dmipy.core.acquisition_scheme.acquisition_scheme_from_gradient_strengths(gradient_strengths, gradient_directions, delta, Delta, TE=None, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Creates an acquisition scheme object from gradient strengths, gradient directions pulse duration \(\delta\) and pulse separation time \(\Delta\).

Parameters:

gradient_strengths: 1D numpy array of shape (Ndata) :

gradient strength of the acquisition in T/m. e.g., a gradient strength of 300 mT/m must be entered as 300 / 1e3 T/m

gradient_directions: 2D numpy array of shape (Ndata, 3) :

gradient directions array of cartesian unit vectors.

delta: float or 1D numpy array of shape (Ndata) :

if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement.

Delta: float or 1D numpy array of shape (Ndata) :

if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement.

min_b_shell_distance : float

minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean.

b0_threshold : float

bvalue threshold for a measurement to be considered a b0 measurement.

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

dmipy.core.acquisition_scheme.acquisition_scheme_from_schemefile(file_path, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Created an acquisition scheme object from a Camino scheme file, containing gradient directions, strengths, pulse duration \(\delta\) and pulse separation time \(\Delta\) and TE.

Parameters:

file_path: string :

absolute file path to schemefile location

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

dmipy.core.acquisition_scheme.unify_length_reference_delta_Delta(reference_array, delta, Delta, TE)

If either delta or Delta are given as float, makes them an array the same size as the reference array.

Parameters:

reference_array : array of size (Nsamples)

typically b-values, q-values or gradient strengths.

delta : float or array of size (Nsamples)

pulse duration in seconds.

Delta : float or array of size (Nsamples)

pulse separation in seconds.

TE : None, float or array of size (Nsamples)

Echo time of the acquisition in seconds.

Returns:

delta_ : array of size (Nsamples)

pulse duration copied to be same size as reference_array

Delta_ : array of size (Nsamples)

pulse separation copied to be same size as reference_array

TE_ : None or array of size (Nsamples)

Echo time copied to be same size as reference_array

dmipy.core.acquisition_scheme.calculate_shell_bvalues_and_indices(bvalues, max_distance=20000000.0)

Calculates which measurements belong to different acquisition shells. It uses scipy’s linkage clustering algorithm, which uses the max_distance input as a limit of including measurements in the same cluster.

For example, if bvalues were [1, 2, 3, 4, 5] and max_distance was 1, then all bvalues would belong to the same cluster. However, if bvalues were [1, 2, 4, 5] max max_distance was 1, then this would result in 2 clusters.

Parameters:

bvalues: 1D numpy array of shape (Ndata) :

bvalues of the acquisition in s/m^2.

max_distance: float :

maximum b-value distance for a measurement to be included in the same shell.

Returns:

shell_indices: 1D numpy array of shape (Ndata) :

array of integers, starting from 0, representing to which shell a measurement belongs. The number itself has no meaning other than just being different for different shells.

shell_bvalues: 1D numpy array of shape (Nshells) :

array of the mean bvalues for every acquisition shell.

dmipy.core.acquisition_scheme.check_acquisition_scheme(bqg_values, gradient_directions, delta, Delta, TE)

function to check the validity of the input parameters.

dmipy.core.acquisition_scheme.gtab_dipy2mipy(dipy_gradient_table)

Converts a dipy gradient_table to a dmipy acquisition_scheme.

dmipy.core.acquisition_scheme.gtab_mipy2dipy(dmipy_gradient_table)

Converts a dmipy acquisition scheme to a dipy gradient_table.

Module: core.constants

Module: core.fitted_modeling_framework

class dmipy.core.fitted_modeling_framework.FittedMultiCompartmentModel(model, S0, mask, fitted_parameters_vector)

The FittedMultiCompartmentModel 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.

Attributes

Methods

R2_coefficient_of_determination(data)

Calculates the R-squared of the model fit.

fitted_and_linked_parameters

Returns the fitted and linked parameters as a dictionary.

fitted_parameters

Returns the fitted parameters as a dictionary.

fod(vertices, visual_odi_lower_bound=0.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.

fod_sh(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.

mean_squared_error(data)

Calculates the mean squared error of the model fit.

peaks_cartesian()

Returns the cartesian peak unit vectors of the model.

peaks_spherical()

Returns the peak angles of the model.

predict(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.

class dmipy.core.fitted_modeling_framework.FittedMultiCompartmentSphericalMeanModel(model, S0, mask, fitted_parameters_vector)

The FittedMultiCompartmentModel 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.

Attributes

Methods

R2_coefficient_of_determination(data)

Calculates the R-squared of the model fit.

fitted_and_linked_parameters

Returns the fitted and linked parameters as a dictionary.

fitted_parameters

Returns the fitted parameters as a dictionary.

mean_squared_error(data)

Calculates the mean squared error of the model fit.

predict(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.

return_parametric_fod_optimizer(distribution='watson', Ncompartments=1)

Retuns parametric FOD optimizer 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 time 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 initial condition for the kernel, and the optimization flags for the kernel are turned off (the kernel will not be fitted while the FOD’s distribution parameters are being optimized).

The function returns a partially evaluated MultiCompartmentModel.fit() instance where the initial_parameters have been set as the spherical mean parameters. Any solver options can then be chosen as for a regular optimization.

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:

parametric_fod_optimizer: MultiCompartmentModel.fit() instance, :

Prepared fit instance of a regular MultiCompartmentModel that can be used to estimate parametric FODs using the fitted spherical mean model as a kernel.

Module: core.gradient_conversions

dmipy.core.gradient_conversions.q_from_b(b, delta, Delta)

Compute q-value from b-value. Units are standard units.

dmipy.core.gradient_conversions.b_from_q(q, delta, Delta)

Compute b-value from q-value. Units are standard units.

dmipy.core.gradient_conversions.q_from_g(g, delta, gyromagnetic_ratio=267513000.0)

Compute q-value from gradient strength. Units are standard units.

dmipy.core.gradient_conversions.g_from_q(q, delta, gyromagnetic_ratio=267513000.0)

Compute gradient strength from q-value. Units are standard units.

dmipy.core.gradient_conversions.b_from_g(g, delta, Delta, gyromagnetic_ratio=267513000.0)

Compute b-value from gradient strength. Units are standard units.

dmipy.core.gradient_conversions.g_from_b(b, delta, Delta, gyromagnetic_ratio=267513000.0)

Compute gradient strength from b-value. Units are standard units.

Module: core.modeling_framework

Document Module

class dmipy.core.modeling_framework.ModelProperties

Contains various properties for CompartmentModels.

Attributes

parameter_cardinality

Returns the cardinality of model parameters

parameter_names

Returns the names of model parameters.

parameter_ranges

Returns the optimization ranges of the model parameters. These ranges are given in O(1) scale so optimization algorithms don’t suffer from large scale differences in optimization parameters.

parameter_scales

Returns the optimization scales for the model parameters. The scales scale the parameter_ranges to their actual size inside optimization algorithms.

parameter_types

Returns the optimization scales for the model parameters. The scales scale the parameter_ranges to their actual size inside optimization algorithms.

class dmipy.core.modeling_framework.MultiCompartmentModelProperties

Class that contains various properties of MultiCompartmentModel instance.

Attributes

Methods

add_linked_parameters_to_parameters(parameters)

When making the MultiCompartmentModel function call, adds the linked parameter to the optimized parameters by evaluating the parameter link function.

bounds_for_optimization

Returns the linear parameter bounds for the model optimization.

opt_params_for_optimization

Returns the linear bools whether to optimize a model parameter.

parameter_initial_guess_to_parameter_vector(**parameters)

Function that returns a parameter_vector while allowing for partial input of model parameters, setting the ones that were not given to ‘None’. Such an array can be given to the fit() function to provide an initial parameter guess when fitting the data to the model.

Parameters:

parameters: keyword arguments of parameter names, :

parameter values of only the parameters you want to give as an initial condition for the optimizer.

Returns:

parameter_vector: array of size (Ndata_x, Ndata_y, …, Nparameters), :

array that contains the linearized model parameters for an ND-array of data voxels, with None’s for non-given parameters.

parameter_names

Returns the names of model parameters.

parameter_vector_to_parameters(parameter_vector)

Returns the model parameters in dictionary format according to their parameter_names. Takes parameter_vector as input, which is the same as the output of a FittedMultiCompartmentModel.fitted_parameter_vector.

Parameters:

parameter_vector: array of size (Ndata_x, Ndata_y, …, Nparameters), :

array that contains the linearized model parameters for an ND-array of data voxels.

Returns:

parameter: dictionary with parameter_names as parameter keys, :

contains the model parameters in dictionary format.

parameters_to_parameter_vector(**parameters)

Returns the model parameters in array format. The input is a parameters dictionary that has parameter_names as keys. This is also the output of a FittedMultiCompartmentModel.fitted_parameters.

It’s possible to give an array of values for one parameter and only a float for others. The function will automatically assume that that the float parameters are constant in the data set and broadcast them accordingly.

The output parameter_vector can be used in simulate_data() to generate data according to the given input parameters.

Parameters:

parameters: keyword arguments of parameter_names. :

Can be given as **parameter_dictionary that contains the model parameter values.

Returns:

parameter_vector: array of size (Ndata_x, Ndata_y, …, Nparameters), :

array that contains the linearized model parameters for an ND-array of data voxels.

scales_for_optimization

Returns the linear parameter scales for model optimization.

set_equal_parameter(parameter_name_in, parameter_name_out)

Allows the user to set two parameters equal to each other. This is used for example in the NODDI model to set the parallel diffusivities of the Stick and Zeppelin compartment to the same value.

The second input parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

parameter_name_in: string :

the first parameter name, see self.parameter_names.

parameter_name_out: string, :

the second parameter name, see self.parameter_names. This is the parameter that will be removed form the optimzed parameters.

set_fixed_parameter(parameter_name, value)

Allows the user to fix an optimization parameter to a static value. The fixed parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

parameter_name: string :

name of the to-be-fixed parameters, see self.parameter_names.

value: float or list of corresponding parameter_cardinality. :

the value to fix the parameter at in SI units.

set_fractional_parameter(parameter1_smaller_equal_than, parameter2)

Allows to impose a constraint to make one parameter smaller or equal to another parameter. This is done by replacing parameter1 with a new parameter that is defined as a fraction between 0 and 1 of parameter2. The new parameter will be the same as the old parameter name with “_fraction” appended to it.

Parameters:

parameter1_smaller_equal_than: string :

parameter name to be made a fraction of parameter2

parameter2: string :

the parameter that is larger or equal than parameter1

set_initial_guess_parameter(parameter_name, value)

Allows the user to fix an optimization parameter to a static value. The fixed parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

parameter_name: string :

name of the to-be-fixed parameters, see self.parameter_names.

value: float or list of corresponding parameter_cardinality. :

the value to fix the parameter at in SI units.

set_tortuous_parameter(lambda_perp_parameter_name, lambda_par_parameter_name, volume_fraction_intra_parameter_name, volume_fraction_extra_parameter_name)

Allows the user to set a tortuosity constraint on the perpendicular diffusivity of the extra-axonal compartment, which depends on the intra-axonal volume fraction and parallel diffusivity.

The perpendicular diffusivity parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

lambda_perp_parameter_name: string :

name of the perpendicular diffusivity parameter, see self.parameter_names.

lambda_par_parameter_name: string :

name of the parallel diffusivity parameter, see self.parameter_names.

volume_fraction_intra_parameter_name: string :

name of the intra-axonal volume fraction parameter, see self.parameter_names.

volume_fraction_extra_parameter_name: string :

name of the extra-axonal volume fraction parameter, see self.parameter_names.

class dmipy.core.modeling_framework.MultiCompartmentModel(models, parameter_links=None)

The MultiCompartmentModel 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.

Attributes

Methods

fit(acquisition_scheme, data, mask=None, solver='brute2fine', Ns=5, maxiter=300, N_sphere_samples=30, use_parallel_processing=False, 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.

simulate_signal(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.

class dmipy.core.modeling_framework.MultiCompartmentSphericalMeanModel(models, parameter_links=None)

The MultiCompartmentModel 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.

Attributes

Methods

fit(acquisition_scheme, data, mask=None, solver='brute2fine', Ns=5, maxiter=300, N_sphere_samples=30, use_parallel_processing=False, 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.

simulate_signal(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.

dmipy.core.modeling_framework.homogenize_x0_to_data(data, x0)

Function that checks if data and initial guess x0 are of the same size. If x0 is 1D, it will be tiled to be the same size as data.

class dmipy.core.modeling_framework.ReturnFixedValue(value)

Parameter fixing class for parameter links.

Methods

DmipyAcquisitionScheme

class dmipy.core.acquisition_scheme.DmipyAcquisitionScheme(bvalues, gradient_directions, qvalues, gradient_strengths, delta, Delta, TE, min_b_shell_distance, b0_threshold)

Class that calculates and contains all information needed to simulate and fit data using microstructure models.

Attributes

Methods

__init__(bvalues, gradient_directions, qvalues, gradient_strengths, delta, Delta, TE, min_b_shell_distance, b0_threshold)
print_acquisition_info

prints a small summary of the acquisition scheme. Is useful to check if the function correctly separated the shells and if the input parameters were given in the right scale.

to_schemefile(filename)

Exports acquisition scheme information in schemefile format, which can be used by the Camino Monte-Carlo simulator.

Parameters:

filename : string,

location at which to save the schemefile.

visualise_acquisition_G_Delta_rainbow(Delta_start=None, Delta_end=None, G_start=None, G_end=None, bval_isolines=array([ 0, 250, 1000, 2500, 5000, 7500, 10000, 14000]), alpha_shading=0.6)

This function visualizes a q-tau acquisition scheme as a function of gradient strength and pulse separation (big_delta). It represents every measurements at its G and big_delta position regardless of b-vector, with a background of b-value isolines for reference. It assumes there is only one unique pulse length (small_delta) in the acquisition scheme.

Parameters:

Delta_start : float,

optional minimum big_delta that is plotted in seconds

Delta_end : float,

optional maximum big_delta that is plotted in seconds

G_start : float,

optional minimum gradient strength that is plotted in T/m

G_end : float,

optional maximum gradient strength taht is plotted in T/m

bval_isolines : array,

optional array of bvalue isolines that are plotted in background given in s/mm^2

alpha_shading : float between [0-1]

optional shading of the bvalue colors in the background

GradientTable

class dmipy.core.acquisition_scheme.GradientTable(gradients, big_delta=None, small_delta=None, b0_threshold=0)

Bases: object

Diffusion gradient information

Parameters:

gradients : array_like (N, 3)

Diffusion gradients. The direction of each of these vectors corresponds to the b-vector, and the length corresponds to the b-value.

b0_threshold : float

Gradients with b-value less than or equal to b0_threshold are considered as b0s i.e. without diffusion weighting.

See also

gradient_table

Notes

The GradientTable object is immutable. Do NOT assign attributes. If you have your gradient table in a bval & bvec format, we recommend using the factory function gradient_table

Attributes

gradients (N,3) ndarray diffusion gradients
qvals: (N,) ndarray   The q-value for each gradient direction. Needs big and small delta.
b0_threshold float Gradients with b-value less than or equal to b0_threshold are considered to not have diffusion weighting.

Methods

__init__(gradients, big_delta=None, small_delta=None, b0_threshold=0)

Constructor for GradientTable class

b0s_mask()
bvals()
bvecs()
info
qvals()

RotationalHarmonicsAcquisitionScheme

class dmipy.core.acquisition_scheme.RotationalHarmonicsAcquisitionScheme(dmipy_acquisition_scheme, N_angular_samples=10)

AcquisitionScheme instance that contains the information necessary to calculate the rotational harmonics for a model for every acquisition shell. It is instantiated using a regular DmipyAcquisitionScheme and N_angular_samples determines how many samples are taken between mu=[0., 0.] and mu=[np.pi/2, 0.].

Parameters:

dmipy_acquisition_scheme: DmipyAcquisitionScheme instance :

An acquisition scheme that has been instantiated using dMipy.

N_angular_samples: int :

Integer representing the number of angular samples per shell.

__init__(dmipy_acquisition_scheme, N_angular_samples=10)

SphericalMeanAcquisitionScheme

class dmipy.core.acquisition_scheme.SphericalMeanAcquisitionScheme(bvalues, qvalues, gradient_strengths, Deltas, deltas)

Acquisition scheme for isotropic spherical mean models.

__init__(bvalues, qvalues, gradient_strengths, Deltas, deltas)

acquisition_scheme_from_bvalues

dmipy.core.acquisition_scheme.acquisition_scheme_from_bvalues(bvalues, gradient_directions, delta, Delta, TE=None, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Creates an acquisition scheme object from bvalues, gradient directions, pulse duration \(\delta\) and pulse separation time \(\Delta\).

Parameters:

bvalues: 1D numpy array of shape (Ndata) :

bvalues of the acquisition in s/m^2. e.g., a bvalue of 1000 s/mm^2 must be entered as 1000 * 1e6 s/m^2

gradient_directions: 2D numpy array of shape (Ndata, 3) :

gradient directions array of cartesian unit vectors.

delta: float or 1D numpy array of shape (Ndata) :

if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement.

Delta: float or 1D numpy array of shape (Ndata) :

if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement.

min_b_shell_distance : float

minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean.

b0_threshold : float

bvalue threshold for a measurement to be considered a b0 measurement.

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

acquisition_scheme_from_gradient_strengths

dmipy.core.acquisition_scheme.acquisition_scheme_from_gradient_strengths(gradient_strengths, gradient_directions, delta, Delta, TE=None, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Creates an acquisition scheme object from gradient strengths, gradient directions pulse duration \(\delta\) and pulse separation time \(\Delta\).

Parameters:

gradient_strengths: 1D numpy array of shape (Ndata) :

gradient strength of the acquisition in T/m. e.g., a gradient strength of 300 mT/m must be entered as 300 / 1e3 T/m

gradient_directions: 2D numpy array of shape (Ndata, 3) :

gradient directions array of cartesian unit vectors.

delta: float or 1D numpy array of shape (Ndata) :

if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement.

Delta: float or 1D numpy array of shape (Ndata) :

if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement.

min_b_shell_distance : float

minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean.

b0_threshold : float

bvalue threshold for a measurement to be considered a b0 measurement.

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

acquisition_scheme_from_qvalues

dmipy.core.acquisition_scheme.acquisition_scheme_from_qvalues(qvalues, gradient_directions, delta, Delta, TE=None, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Creates an acquisition scheme object from qvalues, gradient directions, pulse duration \(\delta\) and pulse separation time \(\Delta\).

Parameters:

qvalues: 1D numpy array of shape (Ndata) :

diffusion sensitization of the acquisition in 1/m. e.g. a qvalue of 10 1/mm must be entered as 10 * 1e3 1/m

gradient_directions: 2D numpy array of shape (Ndata, 3) :

gradient directions array of cartesian unit vectors.

delta: float or 1D numpy array of shape (Ndata) :

if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement.

Delta: float or 1D numpy array of shape (Ndata) :

if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement.

min_b_shell_distance : float

minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean.

b0_threshold : float

bvalue threshold for a measurement to be considered a b0 measurement.

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

acquisition_scheme_from_schemefile

dmipy.core.acquisition_scheme.acquisition_scheme_from_schemefile(file_path, min_b_shell_distance=50000000.0, b0_threshold=10000000.0)

Created an acquisition scheme object from a Camino scheme file, containing gradient directions, strengths, pulse duration \(\delta\) and pulse separation time \(\Delta\) and TE.

Parameters:

file_path: string :

absolute file path to schemefile location

Returns:

DmipyAcquisitionScheme: acquisition scheme object :

contains all information of the acquisition scheme to be used in any microstructure model.

b_from_g

dmipy.core.acquisition_scheme.b_from_g(g, delta, Delta, gyromagnetic_ratio=267513000.0)

Compute b-value from gradient strength. Units are standard units.

b_from_q

dmipy.core.acquisition_scheme.b_from_q(q, delta, Delta)

Compute b-value from q-value. Units are standard units.

calculate_shell_bvalues_and_indices

dmipy.core.acquisition_scheme.calculate_shell_bvalues_and_indices(bvalues, max_distance=20000000.0)

Calculates which measurements belong to different acquisition shells. It uses scipy’s linkage clustering algorithm, which uses the max_distance input as a limit of including measurements in the same cluster.

For example, if bvalues were [1, 2, 3, 4, 5] and max_distance was 1, then all bvalues would belong to the same cluster. However, if bvalues were [1, 2, 4, 5] max max_distance was 1, then this would result in 2 clusters.

Parameters:

bvalues: 1D numpy array of shape (Ndata) :

bvalues of the acquisition in s/m^2.

max_distance: float :

maximum b-value distance for a measurement to be included in the same shell.

Returns:

shell_indices: 1D numpy array of shape (Ndata) :

array of integers, starting from 0, representing to which shell a measurement belongs. The number itself has no meaning other than just being different for different shells.

shell_bvalues: 1D numpy array of shape (Nshells) :

array of the mean bvalues for every acquisition shell.

check_acquisition_scheme

dmipy.core.acquisition_scheme.check_acquisition_scheme(bqg_values, gradient_directions, delta, Delta, TE)

function to check the validity of the input parameters.

fcluster

dmipy.core.acquisition_scheme.fcluster(Z, t, criterion='inconsistent', depth=2, R=None, monocrit=None)

Form flat clusters from the hierarchical clustering defined by the given linkage matrix.

Parameters:

Z : ndarray

The hierarchical clustering encoded with the matrix returned by the linkage function.

t : float

The threshold to apply when forming flat clusters.

criterion : str, optional

The criterion to use in forming flat clusters. This can be any of the following values:

inconsistent : If a cluster node and all its

descendants have an inconsistent value less than or equal to t then all its leaf descendants belong to the same flat cluster. When no non-singleton cluster meets this criterion, every node is assigned to its own cluster. (Default)

distance : Forms flat clusters so that the original

observations in each flat cluster have no greater a cophenetic distance than t.

maxclust : Finds a minimum threshold r so that

the cophenetic distance between any two original observations in the same flat cluster is no more than r and no more than t flat clusters are formed.

monocrit : Forms a flat cluster from a cluster node c

with index i when monocrit[j] <= t.

For example, to threshold on the maximum mean distance as computed in the inconsistency matrix R with a threshold of 0.8 do:

MR = maxRstat(Z, R, 3)
cluster(Z, t=0.8, criterion='monocrit', monocrit=MR)
maxclust_monocrit : Forms a flat cluster from a

non-singleton cluster node c when monocrit[i] <= r for all cluster indices i below and including c. r is minimized such that no more than t flat clusters are formed. monocrit must be monotonic. For example, to minimize the threshold t on maximum inconsistency values so that no more than 3 flat clusters are formed, do:

MI = maxinconsts(Z, R)
cluster(Z, t=3, criterion='maxclust_monocrit', monocrit=MI)

depth : int, optional

The maximum depth to perform the inconsistency calculation. It has no meaning for the other criteria. Default is 2.

R : ndarray, optional

The inconsistency matrix to use for the ‘inconsistent’ criterion. This matrix is computed if not provided.

monocrit : ndarray, optional

An array of length n-1. monocrit[i] is the statistics upon which non-singleton i is thresholded. The monocrit vector must be monotonic, i.e. given a node c with index i, for all node indices j corresponding to nodes below c, monocrit[i] >= monocrit[j].

Returns:

fcluster : ndarray

An array of length n. T[i] is the flat cluster number to which original observation i belongs.

g_from_b

dmipy.core.acquisition_scheme.g_from_b(b, delta, Delta, gyromagnetic_ratio=267513000.0)

Compute gradient strength from b-value. Units are standard units.

g_from_q

dmipy.core.acquisition_scheme.g_from_q(q, delta, gyromagnetic_ratio=267513000.0)

Compute gradient strength from q-value. Units are standard units.

get_sh_order_from_bval

dmipy.core.acquisition_scheme.get_sh_order_from_bval(bval)

Estimates minimum sh_order to represent data of given b-value.

gradient_table

dmipy.core.acquisition_scheme.gradient_table(bvals, bvecs=None, big_delta=None, small_delta=None, b0_threshold=0, atol=0.01)

A general function for creating diffusion MR gradients.

It reads, loads and prepares scanner parameters like the b-values and b-vectors so that they can be useful during the reconstruction process.

Parameters:

bvals : can be any of the four options

  1. an array of shape (N,) or (1, N) or (N, 1) with the b-values.
  2. a path for the file which contains an array like the above (1).
  3. an array of shape (N, 4) or (4, N). Then this parameter is considered to be a b-table which contains both bvals and bvecs. In this case the next parameter is skipped.
  4. a path for the file which contains an array like the one at (3).

bvecs : can be any of two options

  1. an array of shape (N, 3) or (3, N) with the b-vectors.
  2. a path for the file which contains an array like the previous.

big_delta : float

acquisition timing duration (default None)

small_delta : float

acquisition timing duration (default None)

b0_threshold : float

All b-values with values less than or equal to bo_threshold are considered as b0s i.e. without diffusion weighting.

atol : float

All b-vectors need to be unit vectors up to a tolerance.

Returns:

gradients : GradientTable

A GradientTable with all the gradient information.

Notes

  1. Often b0s (b-values which correspond to images without diffusion weighting) have 0 values however in some cases the scanner cannot provide b0s of an exact 0 value and it gives a bit higher values e.g. 6 or 12. This is the purpose of the b0_threshold in the __init__.
  2. We assume that the minimum number of b-values is 7.
  3. B-vectors should be unit vectors.

Examples

>>> from dipy.core.gradients import gradient_table
>>> bvals=1500*np.ones(7)
>>> bvals[0]=0
>>> sq2=np.sqrt(2)/2
>>> bvecs=np.array([[0, 0, 0],
...                 [1, 0, 0],
...                 [0, 1, 0],
...                 [0, 0, 1],
...                 [sq2, sq2, 0],
...                 [sq2, 0, sq2],
...                 [0, sq2, sq2]])
>>> gt = gradient_table(bvals, bvecs)
>>> gt.bvecs.shape == bvecs.shape
True
>>> gt = gradient_table(bvals, bvecs.T)
>>> gt.bvecs.shape == bvecs.T.shape
False

gtab_dipy2mipy

dmipy.core.acquisition_scheme.gtab_dipy2mipy(dipy_gradient_table)

Converts a dipy gradient_table to a dmipy acquisition_scheme.

gtab_mipy2dipy

dmipy.core.acquisition_scheme.gtab_mipy2dipy(dmipy_gradient_table)

Converts a dmipy acquisition scheme to a dipy gradient_table.

linkage

dmipy.core.acquisition_scheme.linkage(y, method='single', metric='euclidean', optimal_ordering=False)

Perform hierarchical/agglomerative clustering.

The input y may be either a 1d compressed distance matrix or a 2d array of observation vectors.

If y is a 1d compressed distance matrix, then y must be a \({n \choose 2}\) sized vector where n is the number of original observations paired in the distance matrix. The behavior of this function is very similar to the MATLAB linkage function.

A \((n-1)\) by 4 matrix Z is returned. At the \(i\)-th iteration, clusters with indices Z[i, 0] and Z[i, 1] are combined to form cluster \(n + i\). A cluster with an index less than \(n\) corresponds to one of the \(n\) original observations. The distance between clusters Z[i, 0] and Z[i, 1] is given by Z[i, 2]. The fourth value Z[i, 3] represents the number of original observations in the newly formed cluster.

The following linkage methods are used to compute the distance \(d(s, t)\) between two clusters \(s\) and \(t\). The algorithm begins with a forest of clusters that have yet to be used in the hierarchy being formed. When two clusters \(s\) and \(t\) from this forest are combined into a single cluster \(u\), \(s\) and \(t\) are removed from the forest, and \(u\) is added to the forest. When only one cluster remains in the forest, the algorithm stops, and this cluster becomes the root.

A distance matrix is maintained at each iteration. The d[i,j] entry corresponds to the distance between cluster \(i\) and \(j\) in the original forest.

At each iteration, the algorithm must update the distance matrix to reflect the distance of the newly formed cluster u with the remaining clusters in the forest.

Suppose there are \(|u|\) original observations \(u[0], \ldots, u[|u|-1]\) in cluster \(u\) and \(|v|\) original objects \(v[0], \ldots, v[|v|-1]\) in cluster \(v\). Recall \(s\) and \(t\) are combined to form cluster \(u\). Let \(v\) be any remaining cluster in the forest that is not \(u\).

The following are methods for calculating the distance between the newly formed cluster \(u\) and each \(v\).

  • method=’single’ assigns

    \[d(u,v) = \min(dist(u[i],v[j]))\]

    for all points \(i\) in cluster \(u\) and \(j\) in cluster \(v\). This is also known as the Nearest Point Algorithm.

  • method=’complete’ assigns

    \[d(u, v) = \max(dist(u[i],v[j]))\]

    for all points \(i\) in cluster u and \(j\) in cluster \(v\). This is also known by the Farthest Point Algorithm or Voor Hees Algorithm.

  • method=’average’ assigns

    \[d(u,v) = \sum_{ij} \frac{d(u[i], v[j])} {(|u|*|v|)}\]

    for all points \(i\) and \(j\) where \(|u|\) and \(|v|\) are the cardinalities of clusters \(u\) and \(v\), respectively. This is also called the UPGMA algorithm.

  • method=’weighted’ assigns

    \[d(u,v) = (dist(s,v) + dist(t,v))/2\]

    where cluster u was formed with cluster s and t and v is a remaining cluster in the forest. (also called WPGMA)

  • method=’centroid’ assigns

    \[dist(s,t) = ||c_s-c_t||_2\]

    where \(c_s\) and \(c_t\) are the centroids of clusters \(s\) and \(t\), respectively. When two clusters \(s\) and \(t\) are combined into a new cluster \(u\), the new centroid is computed over all the original objects in clusters \(s\) and \(t\). The distance then becomes the Euclidean distance between the centroid of \(u\) and the centroid of a remaining cluster \(v\) in the forest. This is also known as the UPGMC algorithm.

  • method=’median’ assigns \(d(s,t)\) like the centroid method. When two clusters \(s\) and \(t\) are combined into a new cluster \(u\), the average of centroids s and t give the new centroid \(u\). This is also known as the WPGMC algorithm.

  • method=’ward’ uses the Ward variance minimization algorithm. The new entry \(d(u,v)\) is computed as follows,

    \[d(u,v) = \sqrt{\frac{|v|+|s|} {T}d(v,s)^2 + \frac{|v|+|t|} {T}d(v,t)^2 - \frac{|v|} {T}d(s,t)^2}\]

    where \(u\) is the newly joined cluster consisting of clusters \(s\) and \(t\), \(v\) is an unused cluster in the forest, \(T=|v|+|s|+|t|\), and \(|*|\) is the cardinality of its argument. This is also known as the incremental algorithm.

Warning: When the minimum distance pair in the forest is chosen, there may be two or more pairs with the same minimum distance. This implementation may choose a different minimum than the MATLAB version.

Parameters:

y : ndarray

A condensed distance matrix. A condensed distance matrix is a flat array containing the upper triangular of the distance matrix. This is the form that pdist returns. Alternatively, a collection of \(m\) observation vectors in \(n\) dimensions may be passed as an \(m\) by \(n\) array. All elements of the condensed distance matrix must be finite, i.e. no NaNs or infs.

method : str, optional

The linkage algorithm to use. See the Linkage Methods section below for full descriptions.

metric : str or function, optional

The distance metric to use in the case that y is a collection of observation vectors; ignored otherwise. See the pdist function for a list of valid distance metrics. A custom distance function can also be used.

optimal_ordering : bool, optional

If True, the linkage matrix will be reordered so that the distance between successive leaves is minimal. This results in a more intuitive tree structure when the data are visualized. defaults to False, because this algorithm can be slow, particularly on large datasets [R2]. See also the optimal_leaf_ordering function.

New in version 1.0.0.

Returns:

Z : ndarray

The hierarchical clustering encoded as a linkage matrix.

See also

scipy.spatial.distance.pdist
pairwise distance metrics

Notes

  1. For method ‘single’ an optimized algorithm based on minimum spanning tree is implemented. It has time complexity \(O(n^2)\). For methods ‘complete’, ‘average’, ‘weighted’ and ‘ward’ an algorithm called nearest-neighbors chain is implemented. It also has time complexity \(O(n^2)\). For other methods a naive algorithm is implemented with \(O(n^3)\) time complexity. All algorithms use \(O(n^2)\) memory. Refer to [R1] for details about the algorithms.
  2. Methods ‘centroid’, ‘median’ and ‘ward’ are correctly defined only if Euclidean pairwise metric is used. If y is passed as precomputed pairwise distances, then it is a user responsibility to assure that these distances are in fact Euclidean, otherwise the produced result will be incorrect.

References

[R1](1, 2) Daniel Mullner, “Modern hierarchical, agglomerative clustering algorithms”, :arXiv:`1109.2378v1`.
[R2](1, 2) Ziv Bar-Joseph, David K. Gifford, Tommi S. Jaakkola, “Fast optimal leaf ordering for hierarchical clustering”, 2001. Bioinformatics https://doi.org/10.1093/bioinformatics/17.suppl_1.S22

Examples

>>> from scipy.cluster.hierarchy import dendrogram, linkage
>>> from matplotlib import pyplot as plt
>>> X = [[i] for i in [2, 8, 0, 4, 1, 9, 9, 0]]
>>> Z = linkage(X, 'ward')
>>> fig = plt.figure(figsize=(25, 10))
>>> dn = dendrogram(Z)
>>> Z = linkage(X, 'single')
>>> fig = plt.figure(figsize=(25, 10))
>>> dn = dendrogram(Z)
>>> plt.show()

q_from_b

dmipy.core.acquisition_scheme.q_from_b(b, delta, Delta)

Compute q-value from b-value. Units are standard units.

q_from_g

dmipy.core.acquisition_scheme.q_from_g(g, delta, gyromagnetic_ratio=267513000.0)

Compute q-value from gradient strength. Units are standard units.

real_sym_rh_basis

dmipy.core.acquisition_scheme.real_sym_rh_basis(sh_order, theta, phi)

Samples a real symmetric rotational harmonic basis at point on the sphere

Samples the basis functions up to order sh_order at points on the sphere given by theta and phi. The basis functions are defined here the same way as in fibernavigator, where the real harmonic \(Y^m_n\) is defined to be:

\(Y^0_n\) if m = 0
Parameters:

sh_order : int

even int > 0, max spherical harmonic degree

theta : float [0, 2*pi]

The azimuthal (longitudinal) coordinate.

phi : float [0, pi]

The polar (colatitudinal) coordinate.

Returns:

real_rh_matrix : array of shape ()

The real harmonic \(Y^0_n\) sampled at theta and phi

real_sym_sh_mrtrix

dmipy.core.acquisition_scheme.real_sym_sh_mrtrix(sh_order, theta, phi)

Compute real spherical harmonics as in mrtrix, where the real harmonic \(Y^m_n\) is defined to be:

Real(:math:`Y^m_n`)       if m > 0
:math:`Y^0_n`             if m = 0
Imag(:math:`Y^|m|_n`)     if m < 0

This may take scalar or array arguments. The inputs will be broadcasted against each other.

Parameters:

sh_order : int

The maximum degree or the spherical harmonic basis.

theta : float [0, pi]

The polar (colatitudinal) coordinate.

phi : float [0, 2*pi]

The azimuthal (longitudinal) coordinate.

Returns:

y_mn : real float

The real harmonic \(Y^m_n\) sampled at theta and phi as implemented in mrtrix. Warning: the basis is Tournier et al 2004 and 2007 is slightly different.

m : array

The order of the harmonics.

n : array

The degree of the harmonics.

unify_length_reference_delta_Delta

dmipy.core.acquisition_scheme.unify_length_reference_delta_Delta(reference_array, delta, Delta, TE)

If either delta or Delta are given as float, makes them an array the same size as the reference array.

Parameters:

reference_array : array of size (Nsamples)

typically b-values, q-values or gradient strengths.

delta : float or array of size (Nsamples)

pulse duration in seconds.

Delta : float or array of size (Nsamples)

pulse separation in seconds.

TE : None, float or array of size (Nsamples)

Echo time of the acquisition in seconds.

Returns:

delta_ : array of size (Nsamples)

pulse duration copied to be same size as reference_array

Delta_ : array of size (Nsamples)

pulse separation copied to be same size as reference_array

TE_ : None or array of size (Nsamples)

Echo time copied to be same size as reference_array

warn

dmipy.core.acquisition_scheme.warn()

Issue a warning, or maybe ignore it or raise an exception.

FittedMultiCompartmentModel

class dmipy.core.fitted_modeling_framework.FittedMultiCompartmentModel(model, S0, mask, fitted_parameters_vector)

The FittedMultiCompartmentModel 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.

Attributes

Methods

__init__(model, S0, mask, fitted_parameters_vector)
R2_coefficient_of_determination(data)

Calculates the R-squared of the model fit.

fitted_and_linked_parameters

Returns the fitted and linked parameters as a dictionary.

fitted_parameters

Returns the fitted parameters as a dictionary.

fod(vertices, visual_odi_lower_bound=0.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.

fod_sh(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.

mean_squared_error(data)

Calculates the mean squared error of the model fit.

peaks_cartesian()

Returns the cartesian peak unit vectors of the model.

peaks_spherical()

Returns the peak angles of the model.

predict(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.

FittedMultiCompartmentSphericalMeanModel

class dmipy.core.fitted_modeling_framework.FittedMultiCompartmentSphericalMeanModel(model, S0, mask, fitted_parameters_vector)

The FittedMultiCompartmentModel 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.

Attributes

Methods

__init__(model, S0, mask, fitted_parameters_vector)
R2_coefficient_of_determination(data)

Calculates the R-squared of the model fit.

fitted_and_linked_parameters

Returns the fitted and linked parameters as a dictionary.

fitted_parameters

Returns the fitted parameters as a dictionary.

mean_squared_error(data)

Calculates the mean squared error of the model fit.

predict(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.

return_parametric_fod_optimizer(distribution='watson', Ncompartments=1)

Retuns parametric FOD optimizer 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 time 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 initial condition for the kernel, and the optimization flags for the kernel are turned off (the kernel will not be fitted while the FOD’s distribution parameters are being optimized).

The function returns a partially evaluated MultiCompartmentModel.fit() instance where the initial_parameters have been set as the spherical mean parameters. Any solver options can then be chosen as for a regular optimization.

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:

parametric_fod_optimizer: MultiCompartmentModel.fit() instance, :

Prepared fit instance of a regular MultiCompartmentModel that can be used to estimate parametric FODs using the fitted spherical mean model as a kernel.

partial

class dmipy.core.fitted_modeling_framework.partial

Bases: object

partial(func, *args, **keywords) - new function with partial application of the given arguments and keywords.

Methods

__init__()

x.__init__(…) initializes x; see help(type(x)) for signature

args

tuple of arguments to future partial calls

func

function object to use in future partial calls

keywords

dictionary of keyword arguments to future partial calls

T1_tortuosity

dmipy.core.fitted_modeling_framework.T1_tortuosity(lambda_par, vf_intra, vf_extra=None)

Tortuosity model for perpendicular extra-axonal diffusivity [1, 2, 3]. If vf_extra=None, then vf_intra must be a nested volume fraction, in the sense that E_bundle = vf_intra * E_intra + (1 - vf_intra) * E_extra, with vf_intra + (1 - vf_intra) = 1. If both vf_intra and vf_extra are given, then they have be be normalized fractions, in the sense that vf_intra + vf_extra <= 1.

Parameters:

lambda_par : float,

parallel diffusivity.

vf_intra : float,

intra-axonal volume fraction [0, 1].

vf_extra : float, (optional)

extra-axonal volume fraction [0, 1].

Returns:

lambda_perp : float,

Rotation matrix.

References :

——- :

.. [1] Bruggeman, Von DAG. “Berechnung verschiedener physikalischer :

Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen.” Annalen der physik 416.7 (1935): 636-664.

.. [2] Sen et al. “A self-similar model for sedimentary rocks with :

application to the dielectric constant of fused glass beads.” Geophysics 46.5 (1981): 781-795.

.. [3] Szafer et al. “Theoretical model for water diffusion in tissues.” :

Magnetic resonance in medicine 33.5 (1995): 697-712.

estimate_spherical_mean_multi_shell

dmipy.core.fitted_modeling_framework.estimate_spherical_mean_multi_shell(E_attenuation, acquisition_scheme, sh_order=6)

Estimates the spherical mean per shell of multi-shell acquisition. Uses spherical harmonics to do the estimation.

Parameters:

E_attenuation : array, shape (N),

signal attenuation.

bvecs : array, shape (N x 3),

x, y, z components of cartesian unit b-vectors.

b_shell_indices : array, shape (N)

array of integers indicating which measurement belongs to which shell. 0 should be used for b0 measurements, 1 for the lowest b-value shell, 2 for the second lowest etc.

Returns:

E_mean : array, shape (number of b-shells)

spherical means of the signal attenuation per shell of the. For example, if there are three shells in the acquisition then the array is of length 3.

fractional_parameter

dmipy.core.fitted_modeling_framework.fractional_parameter(fractional_parameter, parameter2)

Defines fractional parameter link where the original parameter ranges between 0-1 of parameter2.

Parameters:

fractional_parameter : string

new parameter that ranges between 0 and 1.

parameter2 : string

parameter that is larger or equal than the original parameter

get_sphere

dmipy.core.fitted_modeling_framework.get_sphere(name='symmetric362')

provide triangulated spheres

Parameters:

name : str

which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’

Returns:

sphere : a dipy.core.sphere.Sphere class instance

Examples

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') 
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

sh_to_sf_matrix

dmipy.core.fitted_modeling_framework.sh_to_sf_matrix(sphere, sh_order, basis_type=None, return_inv=True, smooth=0)

Matrix that transforms Spherical harmonics (SH) to spherical function (SF).

Parameters:

sphere : Sphere

The points on which to sample the spherical function.

sh_order : int, optional

Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order_2) / 2 SH coefficients (default 4).

basis_type : {None, ‘mrtrix’, ‘fibernav’}

None for the default dipy basis, mrtrix for the MRtrix basis, and fibernav for the FiberNavigator basis (default None).

return_inv : bool

If True then the inverse of the matrix is also returned

smooth : float, optional

Lambda-regularization in the SH fit (default 0.0).

Returns:

B : ndarray

Matrix that transforms spherical harmonics to spherical function sf = np.dot(sh, B).

invB : ndarray

Inverse of B.

unitsphere2cart_Nd

dmipy.core.fitted_modeling_framework.unitsphere2cart_Nd(mu)

Optimized function deicated to convert 1D unit sphere coordinates to cartesian coordinates.

Parameters:

mu : Nd array of size (…, 2)

unit sphere coordinates, as theta, phi = mu

Returns:

mu_cart, Nd array of size (…, 3) :

mu in cartesian coordinates, as x, y, z = mu_cart

b_from_g

dmipy.core.gradient_conversions.b_from_g(g, delta, Delta, gyromagnetic_ratio=267513000.0)

Compute b-value from gradient strength. Units are standard units.

b_from_q

dmipy.core.gradient_conversions.b_from_q(q, delta, Delta)

Compute b-value from q-value. Units are standard units.

g_from_b

dmipy.core.gradient_conversions.g_from_b(b, delta, Delta, gyromagnetic_ratio=267513000.0)

Compute gradient strength from b-value. Units are standard units.

g_from_q

dmipy.core.gradient_conversions.g_from_q(q, delta, gyromagnetic_ratio=267513000.0)

Compute gradient strength from q-value. Units are standard units.

q_from_b

dmipy.core.gradient_conversions.q_from_b(b, delta, Delta)

Compute q-value from b-value. Units are standard units.

q_from_g

dmipy.core.gradient_conversions.q_from_g(g, delta, gyromagnetic_ratio=267513000.0)

Compute q-value from gradient strength. Units are standard units.

Brute2FineOptimizer

class dmipy.core.modeling_framework.Brute2FineOptimizer(model, acquisition_scheme, Ns=5)

Brute force optimizer with refining. Essentially this function does both the brute force optimization like GlobalBruteOptimizer (without treating mu differently, which is currently suboptimal), but then follows it with a gradient-descent based refining step [1, 2] to find the local minimum.

All parameters are optimized within their parameter_bounds. Volume fraction are optimized by ‘nesting’ them, such that given a set of models \(m_1...m_N\), and the partial volume ratios math:v_1…v_{N-1}, the partial volume function is

\[v_1 m_1 + (1 - v_1) v_2 m_2 + ... + (1 - v_1)...(1-v_{N-1}) m_N\]
Parameters:

model: dmipy MultiCompartmentModel instance, :

Can be composed of any model combination.

acquisition_scheme: DmipyAcquisitionScheme instance, :

acquisition scheme of the to-be-fitted data.

Ns: integer, :

number of equally spaced sampling points for regular parameters.

References

[R3]Byrd, R H and P Lu and J. Nocedal. 1995. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM Journal on Scientific and Statistical Computing 16 (5): 1190-1208.
[R4]Zhu, C and R H Byrd and J Nocedal. 1997. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization. ACM Transactions on Mathematical Software 23 (4): 550-560.

Methods

__init__(model, acquisition_scheme, Ns=5)
objective_function(parameter_vector, data)

The objective function for brute-force and gradient-based optimizer.

FittedMultiCompartmentModel

class dmipy.core.modeling_framework.FittedMultiCompartmentModel(model, S0, mask, fitted_parameters_vector)

The FittedMultiCompartmentModel 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.

Attributes

Methods

__init__(model, S0, mask, fitted_parameters_vector)
R2_coefficient_of_determination(data)

Calculates the R-squared of the model fit.

fitted_and_linked_parameters

Returns the fitted and linked parameters as a dictionary.

fitted_parameters

Returns the fitted parameters as a dictionary.

fod(vertices, visual_odi_lower_bound=0.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.

fod_sh(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.

mean_squared_error(data)

Calculates the mean squared error of the model fit.

peaks_cartesian()

Returns the cartesian peak unit vectors of the model.

peaks_spherical()

Returns the peak angles of the model.

predict(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.

FittedMultiCompartmentSphericalMeanModel

class dmipy.core.modeling_framework.FittedMultiCompartmentSphericalMeanModel(model, S0, mask, fitted_parameters_vector)

The FittedMultiCompartmentModel 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.

Attributes

Methods

__init__(model, S0, mask, fitted_parameters_vector)
R2_coefficient_of_determination(data)

Calculates the R-squared of the model fit.

fitted_and_linked_parameters

Returns the fitted and linked parameters as a dictionary.

fitted_parameters

Returns the fitted parameters as a dictionary.

mean_squared_error(data)

Calculates the mean squared error of the model fit.

predict(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.

return_parametric_fod_optimizer(distribution='watson', Ncompartments=1)

Retuns parametric FOD optimizer 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 time 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 initial condition for the kernel, and the optimization flags for the kernel are turned off (the kernel will not be fitted while the FOD’s distribution parameters are being optimized).

The function returns a partially evaluated MultiCompartmentModel.fit() instance where the initial_parameters have been set as the spherical mean parameters. Any solver options can then be chosen as for a regular optimization.

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:

parametric_fod_optimizer: MultiCompartmentModel.fit() instance, :

Prepared fit instance of a regular MultiCompartmentModel that can be used to estimate parametric FODs using the fitted spherical mean model as a kernel.

GlobalBruteOptimizer

class dmipy.core.modeling_framework.GlobalBruteOptimizer(model, acquisition_scheme, x0_vector=None, Ns=5, N_sphere_samples=30)

Brute-Force optimizer. Given a model and an acquisition scheme, first computes a global grid of parameters and corresponding signal attenuations. All parameters except the spherical orientation parameter ‘mu’ are sampled between their corresponding parameter_ranges in ‘Ns’ equal steps. For ‘mu’ a spherical grid of ‘N_sphere_samples” points is used, which were generated using the work of Caruyer et al. [1].

When calling the function with measured data, the closest parameters are return based on the sum-squared error between the signal grid and the data.

Parameters:

model: dmipy MultiCompartmentModel instance, :

Can be composed of any model combination.

acquisition_scheme: DmipyAcquisitionScheme instance, :

acquisition scheme of the to-be-fitted data.

x0_vector: array of size (Nparameters,) :

optional initial guess parameters. As long as the initial guess does not vary voxel-by-voxel the parameter grid will be estimated, only including the initial guess value for the parameters that were included in x0_vector.

Ns: integer, :

number of equally spaced sampling points for regular parameters.

N_sphere_sampled: integer, :

number of sampled sphere points to sample orientation ‘mu’.

References

[R5]Caruyer, Emmanuel, et al. “Design of multishell sampling schemes with uniform coverage in diffusion MRI.” Magnetic resonance in medicine 69.6 (2013): 1534-1540.

Methods

__init__(model, acquisition_scheme, x0_vector=None, Ns=5, N_sphere_samples=30)
precompute_signal_grid(model, x0_vector, Ns, N_sphere_samples)

Function that estimates the parameter grid and corresponding signal attenuation.

NOTE: In the current implementation initial guesses for volume fractions are still ignored

MixOptimizer

class dmipy.core.modeling_framework.MixOptimizer(model, acquisition_scheme, maxiter=150)

The stochastic Microstructure In Crossings (MIX) optimizer [1] uses a three-step process to fit the parameters of a multi-compartment (MC) model to data. The key innovation is that they separate linear from non-linear parameters in the fitting process, meaning the linear volume fractions and non-linear other ones(e.g. diffusivities) are optimized at different stages in the process.

In the first step [1] describes using a genetic algorithm to estimate the non-linear parameters of an MC model. For this we use scipy’s differential_evolution algorithm.

In the second step [1] describes using CVX to estimate the linear volume fractions of an MC model. For this we use scipy’s COBYLA algorithm since it allows us to impose the parameter constraints we need for volume fractions; namely that they are positive and sum up to one.

The third and last step in [1] is a refining step to find a local minimum given the solutions of step one and two. For this we use scipy’s gradient-based L-BFGS-B algorithm with nested volume fractions.

The final result takes a model’s parameter_ranges into account, only yielding parameters within their allowed optimization domain.

Parameters:

model : MultiCompartmentModel instance,

A multicompartment model that has been instantiated using dMipy.

acquisition_scheme : DmipyAcquisitionScheme instance,

An acquisition scheme that has been instantiated using dMipy.

maxiter : integer

The maximum allowed iterations for the differential evolution algorithm

References

[R6]Farooq, Hamza, et al. “Microstructure Imaging of Crossing (MIX) White Matter Fibers from diffusion MRI.” Nature Scientific reports 6 (2016).

Methods

__init__(model, acquisition_scheme, maxiter=150)
cobyla_cost_function(fractions, phi, data)

Objective function of linear parameter estimation using COBYLA.

objective_function(optimized_parameter_vector, data, acquisition_scheme, x0_params)

Objective function of final refining step using L-BFGS-B

stochastic_objective_function(optimized_parameter_vector, data, acquisition_scheme, x0_params)

Objective function for stochastic non-linear parameter estimation using differential_evolution

ModelProperties

class dmipy.core.modeling_framework.ModelProperties

Contains various properties for CompartmentModels.

Attributes

parameter_cardinality

Returns the cardinality of model parameters

parameter_names

Returns the names of model parameters.

parameter_ranges

Returns the optimization ranges of the model parameters. These ranges are given in O(1) scale so optimization algorithms don’t suffer from large scale differences in optimization parameters.

parameter_scales

Returns the optimization scales for the model parameters. The scales scale the parameter_ranges to their actual size inside optimization algorithms.

parameter_types

Returns the optimization scales for the model parameters. The scales scale the parameter_ranges to their actual size inside optimization algorithms.

MultiCompartmentModel

class dmipy.core.modeling_framework.MultiCompartmentModel(models, parameter_links=None)

Bases: dmipy.core.modeling_framework.MultiCompartmentModelProperties

The MultiCompartmentModel 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.

Attributes

Methods

__init__(models, parameter_links=None)
fit(acquisition_scheme, data, mask=None, solver='brute2fine', Ns=5, maxiter=300, N_sphere_samples=30, use_parallel_processing=False, 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.

simulate_signal(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.

MultiCompartmentModelProperties

class dmipy.core.modeling_framework.MultiCompartmentModelProperties

Class that contains various properties of MultiCompartmentModel instance.

Attributes

Methods

add_linked_parameters_to_parameters(parameters)

When making the MultiCompartmentModel function call, adds the linked parameter to the optimized parameters by evaluating the parameter link function.

bounds_for_optimization

Returns the linear parameter bounds for the model optimization.

opt_params_for_optimization

Returns the linear bools whether to optimize a model parameter.

parameter_initial_guess_to_parameter_vector(**parameters)

Function that returns a parameter_vector while allowing for partial input of model parameters, setting the ones that were not given to ‘None’. Such an array can be given to the fit() function to provide an initial parameter guess when fitting the data to the model.

Parameters:

parameters: keyword arguments of parameter names, :

parameter values of only the parameters you want to give as an initial condition for the optimizer.

Returns:

parameter_vector: array of size (Ndata_x, Ndata_y, …, Nparameters), :

array that contains the linearized model parameters for an ND-array of data voxels, with None’s for non-given parameters.

parameter_names

Returns the names of model parameters.

parameter_vector_to_parameters(parameter_vector)

Returns the model parameters in dictionary format according to their parameter_names. Takes parameter_vector as input, which is the same as the output of a FittedMultiCompartmentModel.fitted_parameter_vector.

Parameters:

parameter_vector: array of size (Ndata_x, Ndata_y, …, Nparameters), :

array that contains the linearized model parameters for an ND-array of data voxels.

Returns:

parameter: dictionary with parameter_names as parameter keys, :

contains the model parameters in dictionary format.

parameters_to_parameter_vector(**parameters)

Returns the model parameters in array format. The input is a parameters dictionary that has parameter_names as keys. This is also the output of a FittedMultiCompartmentModel.fitted_parameters.

It’s possible to give an array of values for one parameter and only a float for others. The function will automatically assume that that the float parameters are constant in the data set and broadcast them accordingly.

The output parameter_vector can be used in simulate_data() to generate data according to the given input parameters.

Parameters:

parameters: keyword arguments of parameter_names. :

Can be given as **parameter_dictionary that contains the model parameter values.

Returns:

parameter_vector: array of size (Ndata_x, Ndata_y, …, Nparameters), :

array that contains the linearized model parameters for an ND-array of data voxels.

scales_for_optimization

Returns the linear parameter scales for model optimization.

set_equal_parameter(parameter_name_in, parameter_name_out)

Allows the user to set two parameters equal to each other. This is used for example in the NODDI model to set the parallel diffusivities of the Stick and Zeppelin compartment to the same value.

The second input parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

parameter_name_in: string :

the first parameter name, see self.parameter_names.

parameter_name_out: string, :

the second parameter name, see self.parameter_names. This is the parameter that will be removed form the optimzed parameters.

set_fixed_parameter(parameter_name, value)

Allows the user to fix an optimization parameter to a static value. The fixed parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

parameter_name: string :

name of the to-be-fixed parameters, see self.parameter_names.

value: float or list of corresponding parameter_cardinality. :

the value to fix the parameter at in SI units.

set_fractional_parameter(parameter1_smaller_equal_than, parameter2)

Allows to impose a constraint to make one parameter smaller or equal to another parameter. This is done by replacing parameter1 with a new parameter that is defined as a fraction between 0 and 1 of parameter2. The new parameter will be the same as the old parameter name with “_fraction” appended to it.

Parameters:

parameter1_smaller_equal_than: string :

parameter name to be made a fraction of parameter2

parameter2: string :

the parameter that is larger or equal than parameter1

set_initial_guess_parameter(parameter_name, value)

Allows the user to fix an optimization parameter to a static value. The fixed parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

parameter_name: string :

name of the to-be-fixed parameters, see self.parameter_names.

value: float or list of corresponding parameter_cardinality. :

the value to fix the parameter at in SI units.

set_tortuous_parameter(lambda_perp_parameter_name, lambda_par_parameter_name, volume_fraction_intra_parameter_name, volume_fraction_extra_parameter_name)

Allows the user to set a tortuosity constraint on the perpendicular diffusivity of the extra-axonal compartment, which depends on the intra-axonal volume fraction and parallel diffusivity.

The perpendicular diffusivity parameter will be removed from the optimized parameters and added as a linked parameter.

Parameters:

lambda_perp_parameter_name: string :

name of the perpendicular diffusivity parameter, see self.parameter_names.

lambda_par_parameter_name: string :

name of the parallel diffusivity parameter, see self.parameter_names.

volume_fraction_intra_parameter_name: string :

name of the intra-axonal volume fraction parameter, see self.parameter_names.

volume_fraction_extra_parameter_name: string :

name of the extra-axonal volume fraction parameter, see self.parameter_names.

MultiCompartmentSphericalMeanModel

class dmipy.core.modeling_framework.MultiCompartmentSphericalMeanModel(models, parameter_links=None)

Bases: dmipy.core.modeling_framework.MultiCompartmentModelProperties

The MultiCompartmentModel 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.

Attributes

Methods

__init__(models, parameter_links=None)
fit(acquisition_scheme, data, mask=None, solver='brute2fine', Ns=5, maxiter=300, N_sphere_samples=30, use_parallel_processing=False, 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.

simulate_signal(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.

OrderedDict

class dmipy.core.modeling_framework.OrderedDict(*args, **kwds)

Bases: dict

Dictionary that remembers insertion order

Methods

__init__(*args, **kwds)

Initialize an ordered dictionary. The signature is the same as regular dictionaries, but keyword arguments are not recommended because their insertion order is arbitrary.

clear() → None. Remove all items from od.
copy() → a shallow copy of od
classmethod fromkeys(S[, v]) → New ordered dictionary with keys from S.

If not specified, the value defaults to None.

items() → list of (key, value) pairs in od
iteritems()

od.iteritems -> an iterator over the (key, value) pairs in od

iterkeys() → an iterator over the keys in od
itervalues()

od.itervalues -> an iterator over the values in od

keys() → list of keys in od
pop(k[, d]) → v, remove specified key and return the corresponding

value. If key is not found, d is returned if given, otherwise KeyError is raised.

popitem() → (k, v), return and remove a (key, value) pair.

Pairs are returned in LIFO order if last is true or FIFO order if false.

setdefault(k[, d]) → od.get(k,d), also set od[k]=d if k not in od
update([E, ]**F) → None. Update D from mapping/iterable E and F.

If E present and has a .keys() method, does: for k in E: D[k] = E[k] If E present and lacks .keys() method, does: for (k, v) in E: D[k] = v In either case, this is followed by: for k, v in F.items(): D[k] = v

values() → list of values in od
viewitems() → a set-like object providing a view on od's items
viewkeys() → a set-like object providing a view on od's keys
viewvalues() → an object providing a view on od's values

ReturnFixedValue

class dmipy.core.modeling_framework.ReturnFixedValue(value)

Parameter fixing class for parameter links.

Methods

__init__(value)

T1_tortuosity

dmipy.core.modeling_framework.T1_tortuosity(lambda_par, vf_intra, vf_extra=None)

Tortuosity model for perpendicular extra-axonal diffusivity [1, 2, 3]. If vf_extra=None, then vf_intra must be a nested volume fraction, in the sense that E_bundle = vf_intra * E_intra + (1 - vf_intra) * E_extra, with vf_intra + (1 - vf_intra) = 1. If both vf_intra and vf_extra are given, then they have be be normalized fractions, in the sense that vf_intra + vf_extra <= 1.

Parameters:

lambda_par : float,

parallel diffusivity.

vf_intra : float,

intra-axonal volume fraction [0, 1].

vf_extra : float, (optional)

extra-axonal volume fraction [0, 1].

Returns:

lambda_perp : float,

Rotation matrix.

References :

——- :

.. [1] Bruggeman, Von DAG. “Berechnung verschiedener physikalischer :

Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen.” Annalen der physik 416.7 (1935): 636-664.

.. [2] Sen et al. “A self-similar model for sedimentary rocks with :

application to the dielectric constant of fused glass beads.” Geophysics 46.5 (1981): 781-795.

.. [3] Szafer et al. “Theoretical model for water diffusion in tissues.” :

Magnetic resonance in medicine 33.5 (1995): 697-712.

estimate_spherical_mean_multi_shell

dmipy.core.modeling_framework.estimate_spherical_mean_multi_shell(E_attenuation, acquisition_scheme, sh_order=6)

Estimates the spherical mean per shell of multi-shell acquisition. Uses spherical harmonics to do the estimation.

Parameters:

E_attenuation : array, shape (N),

signal attenuation.

bvecs : array, shape (N x 3),

x, y, z components of cartesian unit b-vectors.

b_shell_indices : array, shape (N)

array of integers indicating which measurement belongs to which shell. 0 should be used for b0 measurements, 1 for the lowest b-value shell, 2 for the second lowest etc.

Returns:

E_mean : array, shape (number of b-shells)

spherical means of the signal attenuation per shell of the. For example, if there are three shells in the acquisition then the array is of length 3.

fractional_parameter

dmipy.core.modeling_framework.fractional_parameter(fractional_parameter, parameter2)

Defines fractional parameter link where the original parameter ranges between 0-1 of parameter2.

Parameters:

fractional_parameter : string

new parameter that ranges between 0 and 1.

parameter2 : string

parameter that is larger or equal than the original parameter

homogenize_x0_to_data

dmipy.core.modeling_framework.homogenize_x0_to_data(data, x0)

Function that checks if data and initial guess x0 are of the same size. If x0 is 1D, it will be tiled to be the same size as data.

optional_package

dmipy.core.modeling_framework.optional_package(name, trip_msg=None)

Return package-like thing and module setup for package name

Parameters:

name : str

package name

trip_msg : None or str

message to give when someone tries to use the return package, but we could not import it, and have returned a TripWire object instead. Default message if None.

Returns:

pkg_like : module or TripWire instance

If we can import the package, return it. Otherwise return an object raising an error when accessed

have_pkg : bool

True if import for package was successful, false otherwise

module_setup : function

callable usually set as setup_module in calling namespace, to allow skipping tests.

parameter_equality

dmipy.core.modeling_framework.parameter_equality(param)

Function to force two model parameters to be equal in the optimization

time

dmipy.core.modeling_framework.time() → floating point number

Return the current time in seconds since the Epoch. Fractions of a second may be present if the system clock provides them.