Skip to content

Commit

Permalink
multi-compartment mix optimization now allows fixed volume fractions …
Browse files Browse the repository at this point in the history 
…and error-catch for wrong solver name.
  • Loading branch information
@rutgerfick
rutgerfick committed Jun 14, 2018
1 parent dc43d12 commit 51e3e8a
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 24 deletions.
9 changes: 9 additions & 0 deletions dmipy/core/modeling_framework.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -885,6 +885,9 @@ def fit(self, acquisition_scheme, data,
fit_func = MixOptimizer(self, self.scheme, maxiter)
print('Setup MIX optimizer in {} seconds'.format(
time() - start))
else:
msg = "Unknown solver name {}".format(solver)
raise ValueError(msg)

start = time()
for idx, pos in enumerate(zip(*mask_pos)):
Expand Down Expand Up @@ -1255,6 +1258,9 @@ def fit(self, acquisition_scheme, data,
fit_func = MixOptimizer(self, self.scheme, maxiter)
print('Setup MIX optimizer in {} seconds'.format(
time() - start))
else:
msg = "Unknown solver name {}".format(solver)
raise ValueError(msg)

start = time()
for idx, pos in enumerate(zip(*mask_pos)):
Expand Down Expand Up @@ -1602,6 +1608,9 @@ def fit(self, acquisition_scheme, data, mask=None,
if solver == 'cvxpy':
fit_func = GeneralPurposeCSDOptimizer(
acquisition_scheme, self, x0_, self.sh_order, unity_constraint)
else:
msg = "Unknown solver name {}".format(solver)
raise ValueError(msg)
start = time()
for idx, pos in enumerate(zip(*mask_pos)):
voxel_E = data_[pos] / S0[pos]
Expand Down
68 changes: 44 additions & 24 deletions dmipy/optimizers/mix.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,23 +75,18 @@ def __call__(self, data, x0_vector=np.array([np.nan])):
The fitted MC model parameters using MIX.
"""
bounds = list(self.model.bounds_for_optimization)
bounds_de = bounds
bounds_minimize = bounds[:-1] # nested volume fractions
opt_flag_array = self.model.opt_params_for_optimization
if x0_vector is not None:
if not np.all(np.isnan(x0_vector)):
bounds_de = []
bounds_minimize = []
for i, x0_ in enumerate(x0_vector):
if np.isnan(x0_):
bounds_de.append(bounds[i])
if opt_flag_array[i] is True:
bounds_minimize.append(bounds[i])
bounds_minimize = bounds_minimize[:-1]

# if there is only one model then MIX only uses DE.
bounds = list(self.model.bounds_for_optimization)
if self.Nmodels == 1:
bounds_de = bounds
opt_flag_array = self.model.opt_params_for_optimization
if x0_vector is not None:
if not np.all(np.isnan(x0_vector)):
bounds_de = []
for i, x0_ in enumerate(x0_vector):
if np.isnan(x0_):
bounds_de.append(bounds[i])

# step 1: stochastic optimization on non-linear parameters.
res_de = differential_evolution(
self.stochastic_objective_function,
Expand All @@ -111,8 +106,23 @@ def __call__(self, data, x0_vector=np.array([np.nan])):
return fitted_parameters
# if there is more than 1 model then we do the 3 steps.
if self.Nmodels > 1:
bounds_de = bounds
bounds_minimize = bounds[:-1] # nested volume fractions
opt_flag_array = self.model.opt_params_for_optimization
if x0_vector is not None:
if not np.all(np.isnan(x0_vector)):
bounds_de = []
bounds_minimize = []
for i, x0_ in enumerate(x0_vector[:-self.Nmodels]):
if np.isnan(x0_):
bounds_de.append(bounds[i])
for i, x0_ in enumerate(x0_vector):
if opt_flag_array[i] is True:
bounds_minimize.append(bounds[i])
bounds_minimize = bounds_minimize[:-1]
else:
bounds_de = bounds_de[:-self.Nmodels]
# step 1: stochastic optimization on non-linear parameters.
bounds_de = bounds_de[:-self.Nmodels]
res_de = differential_evolution(self.stochastic_objective_function,
bounds=bounds_de,
maxiter=self.maxiter,
Expand All @@ -123,6 +133,12 @@ def __call__(self, data, x0_vector=np.array([np.nan])):
if np.all(np.isnan(x0_vector)):
parameter_vector = np.r_[optimized_parameter_vector,
np.ones(self.Nmodels)]
elif np.all(~np.isnan(x0_vector[-self.Nmodels:])):
x0_bool_array = ~np.isnan(x0_vector)
parameter_vector = np.ones(len(x0_bool_array))
parameter_vector[~x0_bool_array] = optimized_parameter_vector
parameter_vector[x0_bool_array] = x0_vector[x0_bool_array]
return parameter_vector
else:
x0_bool_array = ~np.isnan(x0_vector)
parameter_vector = np.ones(len(x0_bool_array))
Expand Down Expand Up @@ -209,7 +225,7 @@ def stochastic_objective_function(self, optimized_parameter_vector,
parameter_vector_no_vf[~x0_bool_n0_vf] = (
optimized_parameter_vector)
parameter_vector_no_vf[x0_bool_n0_vf] = x0_params[
x0_bool_array]
:-self.Nmodels][x0_bool_n0_vf]
parameter_vector[:-self.Nmodels] = parameter_vector_no_vf
parameter_vector = (
parameter_vector * self.model.scales_for_optimization)
Expand All @@ -218,13 +234,17 @@ def stochastic_objective_function(self, optimized_parameter_vector,
phi_x = self.model(acquisition_scheme,
quantity="stochastic cost function",
**parameters)
A = np.dot(phi_x.T, phi_x)
try:
phi_inv = np.dot(np.linalg.inv(A), phi_x.T)
vf = np.dot(phi_inv, data)
except np.linalg.linalg.LinAlgError:
# happens when models have the same signal attenuations.
vf = np.ones(self.Nmodels) / float(self.Nmodels)
if np.all(~np.isnan(x0_params[-self.Nmodels:])):
# if initial guess is given for volume fractions
vf = x0_params[-self.Nmodels:]
else:
A = np.dot(phi_x.T, phi_x)
try:
phi_inv = np.dot(np.linalg.inv(A), phi_x.T)
vf = np.dot(phi_inv, data)
except np.linalg.linalg.LinAlgError:
# happens when models have the same signal attenuations.
vf = np.ones(self.Nmodels) / float(self.Nmodels)
E_hat = np.dot(phi_x, vf)
objective = np.dot(data - E_hat, data - E_hat).squeeze()
return objective * 1e5
Expand Down

0 comments on commit 51e3e8a

Please sign in to comment.