The stochastic Microstructure In Crossings (MIX) optimizer (Farooq et al. 2016) 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 this notebook we we first show an overall example of how the optimizer is used, and then go over the three steps as given by (Farooq et al. 2016).
# necessary imports
from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme
from dmipy.core.modeling_framework import MultiCompartmentModel
from dmipy.signal_models import gaussian_models, cylinder_models
import numpy as np
# Import some acquition scheme (Wu-Minn HCP scheme)
scheme = wu_minn_hcp_acquisition_scheme()
# Set up test Ball and Stick model
ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
mc_model = MultiCompartmentModel([stick, ball])
# Make some random model parameters with the appropriate scale
random_parameters = {}
for parameter, card in mc_model.parameter_cardinality.items():
random_parameters[parameter] = (
np.random.rand(card) * mc_model.parameter_scales[parameter])
# Making sure the volume fractions add up to 1.
random_parameters['partial_volume_1'] = 1. - random_parameters['partial_volume_0']
# Making a parameter vector out of the parameters so the mc_model understands it.
random_parameter_vector = mc_model.parameters_to_parameter_vector(**random_parameters)
# Generate data with the HCP scheme and the model parameters
random_signal = mc_model.simulate_signal(scheme, random_parameter_vector)
# Set up the MIX optimizer
from dmipy.optimizers.mix import MixOptimizer
mix = MixOptimizer(mc_model, scheme)
estimated_parameters = mix(random_signal)
# Print the true and estimated parameters in O(1) scale
print "True parameters: ", random_parameter_vector / mc_model.scales_for_optimization
print "Estimated parameters:", estimated_parameters
print "Note that the 'mu' parameters (second and third) might be pi away from each other (antipodal symmetry)."
True parameters: [ 0.46405481 0.30005209 0.03821492 0.64970014 0.48858711 0.51141289] Estimated parameters: [ 0.46407092 0.30005198 0.03821499 0.64970364 0.4886475 0.5113525 ] Note that the 'mu' parameters (second and third) might be pi away from each other (antipodal symmetry).
If you rerun this many times you'll see that most of the time the solver finds the ground truth solution. However, sometimes the random diffusivities will make the Ball and Stick look very similar in terms of signal, in which case it can still find the wrong solution.
In the code you should call MIX as follows:
fitted_mc_model = mc_model.fit(scheme, random_signal, solver='mix')
print "Estimated parameters:", fitted_mc_model.fitted_parameters_vector / mc_model.scales_for_optimization
Using parallel processing with 8 workers. Setup MIX optimizer in 3.09944152832e-05 seconds Fitting of 1 voxels complete in 4.15033698082 seconds. Average of 4.15033698082 seconds per voxel. Estimated parameters: [[ 0.46407092 2.84154068 -3.10337767 0.64970364 0.4886475 0.5113525 ]]
While it often manages to find the global solution for more complex models, it does require many more function evaluations than gradient descent methods.
In the first step a genetic algorithm (GA) is used to estimate the non-linear parameters of an MC model. GAs are inspired by the process of natural selection that belong to the larger class of evolutionary algorithms. Without going too much into the details, this means that GAs do not rely on using gradient-descent on some cost function, but instead stochastically sample pools in the parameter space looking for good candidates. For this we use scipy's differential_evolution (DE) algorithm.
The idea is to separate estimation of the linear from the non-linear parameters. This is done by only sampling the non-linear parameters of components of a multi-compartment model, and then estimating the linear volume fractions using least squares. Dmipy's MultiCompartmentModel instance can be called with the 'stochastic cost function' option to give us the signal attenuation for every model in the MultiCompartmentModel separately for the given parameters, ignoring the volume fraction values and just setting them to one for all models.
separate_model_signals = mc_model(scheme, quantity="stochastic cost function", **random_parameters)
# the first column is the signal for the Stick and the second for the Ball
# You can see that volume fractions are ignored and both models are normalized to 1.
separate_model_signals[:20]
array([[ 1. , 1. ],
[ 0.99384665, 0.62872909],
[ 0.77348182, 0.39530027],
[ 0.21384913, 0.24853678],
[ 0.9378883 , 0.62872909],
[ 0.89676035, 0.24853678],
[ 0.52003482, 0.39530027],
[ 0.5596586 , 0.62872909],
[ 0.68081258, 0.39530027],
[ 0.73836681, 0.24853678],
[ 0.93131368, 0.62872909],
[ 0.28058554, 0.39530027],
[ 0.53834826, 0.24853678],
[ 0.73879139, 0.62872909],
[ 0.90655213, 0.39530027],
[ 0.99504891, 0.24853678],
[ 1. , 1. ],
[ 0.79055492, 0.62872909],
[ 0.99993183, 0.39530027],
[ 0.21916649, 0.24853678]])
The volume fractions are then estimated using least squares
# doing this with the ground truth parmameters will give use the ground truth volume fractions
volume_fractions = np.dot(np.linalg.pinv(separate_model_signals), random_signal)
volume_fractions
array([ 0.48858711, 0.51141289])
# meaning we will find the right solution
estimated_signal = np.dot(separate_model_signals, volume_fractions)
print 'Sum Squared Difference:', np.sum((estimated_signal - random_signal) ** 2)
Sum Squared Difference: 5.06288463781e-30
Normally, the GE algorithms will sample the parameter space to look for the model parameters that will match the given signal. Since there are no other constraints on volume fractions in this step, this means that the least-squares solution for the volume fractions will probably not sum up to one (but will hopefully be pretty close). In the next step we use the found non-linear parameters to find a constrained solution for the volume fractions.
(Farooq et al. 2016) 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.
from scipy.optimize import fmin_cobyla
# defining the two constraints for the volume fractions
def cobyla_positivity_constraint(volume_fractions, *args):
"COBYLA positivity constraint on volume fractions"
return volume_fractions - 0.001
def cobyla_unity_constraint(volume_fractions, *args):
"COBYLA unity constraint on volume fractions"
return np.sum(volume_fractions) - 1
# generate slightly off initial guess
x0 = volume_fractions + (np.random.rand(2) - 0.5) / 5.
# do the COBYLA optimization
cobyla_fractions = fmin_cobyla(func=mix.cobyla_cost_function, x0=x0,
cons=[cobyla_positivity_constraint,cobyla_unity_constraint],
args=(separate_model_signals, random_signal))
print "initial fractions: ", x0
print "COBYLA fractions: ", cobyla_fractions
print "ground truth fractions:", volume_fractions
initial fractions: [ 0.53661881 0.50558657] COBYLA fractions: [ 0.48863674 0.51150103] ground truth fractions: [ 0.48858711 0.51141289]
The third and last step in 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. This is exactly the same as in the brute force algorithm, so we'll refer to that example for more info.