Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[BF] IVIM fixes #1931

Merged
merged 14 commits into from Aug 1, 2019
Merged

[BF] IVIM fixes #1931

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
d676067
initial commit
ShreyasFadnavis Jul 18, 2019
9ddeded
changed lm to trf
ShreyasFadnavis Jul 18, 2019
085689b
del
ShreyasFadnavis Jul 18, 2019
ae472e1
Delete reconst_mcsd.py
ShreyasFadnavis Jul 18, 2019
8f41b28
rem conflicts
ShreyasFadnavis Jul 29, 2019
169658c
bounds
ShreyasFadnavis Jul 29, 2019
5e66b51
fixing pep8
ShreyasFadnavis Jul 29, 2019
69bde0a
fixes
ShreyasFadnavis Jul 29, 2019
ba978fa
Update dipy/reconst/ivim.py
ShreyasFadnavis Jul 30, 2019
e15274c
fixes suggested by Ariel
ShreyasFadnavis Jul 30, 2019
58a9d73
fixing vol frac
ShreyasFadnavis Jul 30, 2019
7682532
addressing Serge's comment
ShreyasFadnavis Jul 30, 2019
2cf1acf
addressing Ariel's comments
ShreyasFadnavis Jul 31, 2019
8442405
adding ref comment for bounds
ShreyasFadnavis Aug 1, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
51 changes: 33 additions & 18 deletions dipy/reconst/ivim.py
View file
Open in desktop
Expand Up @@ -10,6 +10,10 @@
cvxpy, have_cvxpy, _ = optional_package("cvxpy")


# global variable for bounding least_squares in both models
BOUNDS = ([0., 0., 0., 0.], [np.inf, .2, 1., 1.])


def ivim_prediction(params, gtab):
"""The Intravoxel incoherent motion (IVIM) model function.

Expand Down Expand Up @@ -120,37 +124,46 @@ def f_D_star_error(params, gtab, signal, S0, D):
return signal - f_D_star_prediction([f, D_star], gtab, S0, D)


def ivim_model_selector(gtab, fit_method='LM', **kwargs):
def ivim_model_selector(gtab, fit_method='trr', **kwargs):
"""
Selector function to switch between the 2-stage Levenberg-Marquardt based
NLLS fitting method (also containing the linear fit): `LM` and the Variable
Projections based fitting method: `VarPro`.
Selector function to switch between the 2-stage Trust-Region Reflective
based NLLS fitting method (also containing the linear fit): `trr` and the
Variable Projections based fitting method: `varpro`.

Parameters
----------
fit_method : string, optional
The value fit_method can either be 'LM' or 'VarPro'.
default : LM
The value fit_method can either be 'trr' or 'varpro'.
default : trr

"""
bounds_warning = 'Bounds for this fit have been set from experiments '
bounds_warning += 'and literature survey. To change the bounds, please '
bounds_warning += 'input your bounds in model definition...'

if fit_method.lower() == 'lm':
return IvimModelLM(gtab, **kwargs)
if fit_method.lower() == 'trr':
ivimmodel_trr = IvimModelTRR(gtab, **kwargs)
if 'bounds' not in kwargs:
warnings.warn(bounds_warning, UserWarning)
return ivimmodel_trr

elif fit_method.lower() == 'varpro':
return IvimModelVP(gtab, **kwargs)
ivimmodel_vp = IvimModelVP(gtab, **kwargs)
if 'bounds' not in kwargs:
warnings.warn(bounds_warning, UserWarning)
return ivimmodel_vp

else:
opt_msg = 'The fit_method option chosen was not correct. '
opt_msg += 'Using fit_method: LM instead...'
opt_msg += 'Using fit_method: TRR instead...'
warnings.warn(opt_msg, UserWarning)
return IvimModelLM(gtab, **kwargs)
return IvimModelTRR(gtab, **kwargs)


IvimModel = ivim_model_selector


class IvimModelLM(ReconstModel):
class IvimModelTRR(ReconstModel):
"""Ivim model
"""
def __init__(self, gtab, split_b_D=400.0, split_b_S0=200., bounds=None,
Expand Down Expand Up @@ -268,11 +281,11 @@ def __init__(self, gtab, split_b_D=400.0, split_b_S0=200., bounds=None,
'eps': eps, 'maxiter': maxiter}
self.x_scale = x_scale

self.bounds = bounds or ((0., 0., 0., 0.), (np.inf, .3, 1., 1.))
self.bounds = bounds or BOUNDS

@multi_voxel_fit
def fit(self, data):
""" Fit method of the IvimModelLM class.
""" Fit method of the IvimModelTRR class.

The fitting takes place in the following steps: Linear fitting for D
(bvals > `split_b_D` (default: 400)) and store S0_prime. Another linear
Expand Down Expand Up @@ -493,7 +506,7 @@ def _leastsq(self, data, x0):

class IvimModelVP(ReconstModel):

def __init__(self, gtab, maxiter=10, xtol=1e-8):
def __init__(self, gtab, bounds=None, maxiter=10, xtol=1e-8):
r""" Initialize an IvimModelVP class.

The IVIM model assumes that biological tissue includes a volume
Expand Down Expand Up @@ -542,6 +555,7 @@ def __init__(self, gtab, maxiter=10, xtol=1e-8):
self.yhat_perfusion = np.zeros(self.bvals.shape[0])
self.yhat_diffusion = np.zeros(self.bvals.shape[0])
self.exp_phi1 = np.zeros((self.bvals.shape[0], 2))
self.bounds = bounds or (BOUNDS[0][1:], BOUNDS[1][1:])

@multi_voxel_fit
def fit(self, data, bounds_de=None):
Expand Down Expand Up @@ -590,7 +604,7 @@ def fit(self, data, bounds_de=None):
x_f = self.x_and_f_to_x_f(x, f)

# Setting up the bounds for least_squares
bounds = ([0.01, 0.005, 10**-4], [0.3, 0.02, 0.003])
bounds = self.bounds

# Optimizer #3: Nonlinear-Least Squares
res = least_squares(self.nlls_cost, x_f, bounds=(bounds),
Expand Down Expand Up @@ -709,11 +723,12 @@ def cvx_fit(self, signal, phi):

# Create four scalar optimization variables.
f = cvxpy.Variable(2)
# Create four constraints.
# Constraints have been set similar to the MIX paper's
# Supplementary Note 2: Synthetic Data Experiments, experiment 2
constraints = [cvxpy.sum(f) == 1,
f[0] >= 0.011,
f[1] >= 0.011,
f[0] <= 0.29,
f[0] <= self.bounds[1][0],
arokem marked this conversation as resolved.
Show resolved Hide resolved
f[1] <= 0.89]

# Form objective.
Expand Down
87 changes: 44 additions & 43 deletions dipy/reconst/tests/test_ivim.py
View file
Open in desktop
Expand Up @@ -29,8 +29,8 @@


def setup_module():
global gtab, ivim_fit_single, ivim_model_LM, data_single, params_LM, \
data_multi, ivim_params_LM, D_star, D, f, S0, gtab_with_multiple_b0, \
global gtab, ivim_fit_single, ivim_model_trr, data_single, params_trr, \
data_multi, ivim_params_trr, D_star, D, f, S0, gtab_with_multiple_b0, \
noisy_single, mevals, gtab_no_b0, ivim_fit_multi, ivim_model_VP, \
f_VP, D_star_VP, D_VP, params_VP

Expand All @@ -44,7 +44,7 @@ def setup_module():

S0, f, D_star, D = 1000.0, 0.132, 0.00885, 0.000921
# params for a single voxel
params_LM = np.array([S0, f, D_star, D])
params_trr = np.array([S0, f, D_star, D])

mevals = np.array(([D_star, D_star, D_star], [D, D, D]))
# This gives an isotropic signal.
Expand All @@ -57,14 +57,14 @@ def setup_module():
data_multi[0, 0, 0] = data_multi[0, 1, 0] = data_multi[
1, 0, 0] = data_multi[1, 1, 0] = data_single

ivim_params_LM = np.zeros((2, 2, 1, 4))
ivim_params_LM[0, 0, 0] = ivim_params_LM[0, 1, 0] = params_LM
ivim_params_LM[1, 0, 0] = ivim_params_LM[1, 1, 0] = params_LM
ivim_params_trr = np.zeros((2, 2, 1, 4))
ivim_params_trr[0, 0, 0] = ivim_params_trr[0, 1, 0] = params_trr
ivim_params_trr[1, 0, 0] = ivim_params_trr[1, 1, 0] = params_trr

ivim_model_LM = IvimModel(gtab, fit_method='LM')
ivim_model_one_stage = IvimModel(gtab, fit_method='LM')
ivim_fit_single = ivim_model_LM.fit(data_single)
ivim_fit_multi = ivim_model_LM.fit(data_multi)
ivim_model_trr = IvimModel(gtab, fit_method='trr')
ivim_model_one_stage = IvimModel(gtab, fit_method='trr')
ivim_fit_single = ivim_model_trr.fit(data_single)
ivim_fit_multi = ivim_model_trr.fit(data_multi)

ivim_model_one_stage.fit(data_single)
ivim_model_one_stage.fit(data_multi)
Expand All @@ -77,8 +77,9 @@ def setup_module():
gtab_no_b0 = gradient_table(bvals_no_b0, bvecs.T, b0_threshold=0)

bvals_with_multiple_b0 = np.array([0., 0., 0., 0., 40., 60., 80., 100.,
120., 140., 160., 180., 200., 300., 400.,
500., 600., 700., 800., 900., 1000.])
120., 140., 160., 180., 200., 300.,
400., 500., 600., 700., 800., 900.,
1000.])

bvecs_with_multiple_b0 = generate_bvecs(N)
gtab_with_multiple_b0 = gradient_table(bvals_with_multiple_b0,
Expand Down Expand Up @@ -136,7 +137,7 @@ def test_single_voxel_fit():

assert_array_equal(est_signal.shape, data_single.shape)

assert_array_almost_equal(ivim_fit_single.model_params, params_LM)
assert_array_almost_equal(ivim_fit_single.model_params, params_trr)
assert_array_almost_equal(est_signal, data_single)

# Test predict function for single voxel
Expand All @@ -152,11 +153,11 @@ def test_multivoxel():
This is to ensure that the fitting routine takes care of signals packed as
1D, 2D or 3D arrays.
"""
ivim_fit_multi = ivim_model_LM.fit(data_multi)
ivim_fit_multi = ivim_model_trr.fit(data_multi)

est_signal = ivim_fit_multi.predict(gtab, S0=1.)
assert_array_equal(est_signal.shape, data_multi.shape)
assert_array_almost_equal(ivim_fit_multi.model_params, ivim_params_LM)
assert_array_almost_equal(ivim_fit_multi.model_params, ivim_params_trr)
assert_array_almost_equal(est_signal, data_multi)


Expand All @@ -168,13 +169,13 @@ def test_ivim_errors():
and is not supported by the older versions. Initializing an IvimModel
with bounds for older Scipy versions should raise an error.
"""
ivim_model_LM = IvimModel(gtab, bounds=([0., 0., 0., 0.],
[np.inf, 1., 1., 1.]),
fit_method='LM')
ivim_fit = ivim_model_LM.fit(data_multi)
ivim_model_trr = IvimModel(gtab, bounds=([0., 0., 0., 0.],
[np.inf, 1., 1., 1.]),
fit_method='trr')
ivim_fit = ivim_model_trr.fit(data_multi)
est_signal = ivim_fit.predict(gtab, S0=1.)
assert_array_equal(est_signal.shape, data_multi.shape)
assert_array_almost_equal(ivim_fit.model_params, ivim_params_LM)
assert_array_almost_equal(ivim_fit.model_params, ivim_params_trr)
assert_array_almost_equal(est_signal, data_multi)


Expand All @@ -185,12 +186,12 @@ def test_mask():
mask_correct = data_multi[..., 0] > 0.2
mask_not_correct = np.array([[False, True, False], [True, False]])

ivim_fit = ivim_model_LM.fit(data_multi, mask_correct)
ivim_fit = ivim_model_trr.fit(data_multi, mask_correct)
est_signal = ivim_fit.predict(gtab, S0=1.)
assert_array_equal(est_signal.shape, data_multi.shape)
assert_array_almost_equal(est_signal, data_multi)
assert_array_almost_equal(ivim_fit.model_params, ivim_params_LM)
assert_raises(ValueError, ivim_model_LM.fit, data_multi,
assert_array_almost_equal(ivim_fit.model_params, ivim_params_trr)
assert_raises(ValueError, ivim_model_trr.fit, data_multi,
mask=mask_not_correct)


Expand All @@ -208,7 +209,7 @@ def test_with_higher_S0():
# Single voxel data
data_single2 = signal2[0]

ivim_fit = ivim_model_LM.fit(data_single2)
ivim_fit = ivim_model_trr.fit(data_single2)

est_signal = ivim_fit.predict(gtab)
assert_array_equal(est_signal.shape, data_single2.shape)
Expand All @@ -228,7 +229,7 @@ def test_b0_threshold_greater_than0():
bvecs = generate_bvecs(N)
gtab = gradient_table(bvals_b0t, bvecs.T)
with assert_raises(ValueError) as vae:
_ = IvimModel(gtab, fit_method='LM')
_ = IvimModel(gtab, fit_method='trr')
b0_s = "The IVIM model requires a measurement at b==0. As of "
assert b0_s in vae.exception

Expand All @@ -253,7 +254,7 @@ def test_bounds_x0():
x0_test = np.array([1., 0.13, 0.001, 0.0001])
test_signal = ivim_prediction(x0_test, gtab)

ivim_fit = ivim_model_LM.fit(test_signal)
ivim_fit = ivim_model_trr.fit(test_signal)

est_signal = ivim_fit.predict(gtab)
assert_array_equal(est_signal.shape, test_signal.shape)
Expand All @@ -267,11 +268,11 @@ def test_predict():
"""
assert_array_almost_equal(ivim_fit_single.predict(gtab),
data_single)
assert_array_almost_equal(ivim_model_LM.predict
assert_array_almost_equal(ivim_model_trr.predict
(ivim_fit_single.model_params, gtab),
data_single)

ivim_fit_multi = ivim_model_LM.fit(data_multi)
ivim_fit_multi = ivim_model_trr.fit(data_multi)
assert_array_almost_equal(ivim_fit_multi.predict(gtab),
data_multi)

Expand All @@ -285,7 +286,7 @@ def test_fit_object():
assert_array_almost_equal(
ivim_fit_single.__getitem__(0).model_params, 1000.)

ivim_fit_multi = ivim_model_LM.fit(data_multi)
ivim_fit_multi = ivim_model_trr.fit(data_multi)
# Should raise a TypeError if the arguments are not passed as tuple
assert_raises(TypeError, ivim_fit_multi.__getitem__, -.1, 0)
# Should return IndexError if invalid indices are passed
Expand All @@ -305,7 +306,7 @@ def test_shape():
Test if `shape` in `IvimFit` class gives the correct output.
"""
assert_array_equal(ivim_fit_single.shape, ())
ivim_fit_multi = ivim_model_LM.fit(data_multi)
ivim_fit_multi = ivim_model_trr.fit(data_multi)
assert_array_equal(ivim_fit_multi.shape, (2, 2, 1))


Expand All @@ -318,7 +319,7 @@ def test_multiple_b0():
# Single voxel data
data_single = signal[0]

ivim_model_multiple_b0 = IvimModel(gtab_with_multiple_b0, fit_method='LM')
ivim_model_multiple_b0 = IvimModel(gtab_with_multiple_b0, fit_method='trr')

ivim_model_multiple_b0.fit(data_single)
# Test if all signals are positive
Expand All @@ -337,7 +338,7 @@ def test_noisy_fit():
around 135 and D and D_star values are equal. Hence doing a test based on
Scipy version.
"""
model_one_stage = IvimModel(gtab, fit_method='LM')
model_one_stage = IvimModel(gtab, fit_method='trr')
with warnings.catch_warnings(record=True) as w:
fit_one_stage = model_one_stage.fit(noisy_single)
assert_equal(len(w), 3)
Expand All @@ -358,7 +359,7 @@ def test_S0():
"""
assert_array_almost_equal(ivim_fit_single.S0_predicted, S0)
assert_array_almost_equal(ivim_fit_multi.S0_predicted,
ivim_params_LM[..., 0])
ivim_params_trr[..., 0])


def test_perfusion_fraction():
Expand All @@ -367,37 +368,37 @@ def test_perfusion_fraction():
"""
assert_array_almost_equal(ivim_fit_single.perfusion_fraction, f)
assert_array_almost_equal(
ivim_fit_multi.perfusion_fraction, ivim_params_LM[..., 1])
ivim_fit_multi.perfusion_fraction, ivim_params_trr[..., 1])


def test_D_star():
"""
Test if the `IvimFit` class returns the correct D_star
"""
assert_array_almost_equal(ivim_fit_single.D_star, D_star)
assert_array_almost_equal(ivim_fit_multi.D_star, ivim_params_LM[..., 2])
assert_array_almost_equal(ivim_fit_multi.D_star, ivim_params_trr[..., 2])


def test_D():
"""
Test if the `IvimFit` class returns the correct D
"""
assert_array_almost_equal(ivim_fit_single.D, D)
assert_array_almost_equal(ivim_fit_multi.D, ivim_params_LM[..., 3])
assert_array_almost_equal(ivim_fit_multi.D, ivim_params_trr[..., 3])


def test_estimate_linear_fit():
"""
Test the linear estimates considering a single exponential fit.
"""
data_single_exponential_D = single_exponential(S0, D, gtab.bvals)
assert_array_almost_equal(ivim_model_LM.estimate_linear_fit(
assert_array_almost_equal(ivim_model_trr.estimate_linear_fit(
data_single_exponential_D,
split_b=500.,
less_than=False),
(S0, D))
data_single_exponential_D_star = single_exponential(S0, D_star, gtab.bvals)
assert_array_almost_equal(ivim_model_LM.estimate_linear_fit(
assert_array_almost_equal(ivim_model_trr.estimate_linear_fit(
data_single_exponential_D_star,
split_b=100.,
less_than=True),
Expand All @@ -410,9 +411,9 @@ def test_estimate_f_D_star():
non-linear fit.
"""
params_f_D = f + 0.001, D + 0.0001
assert_array_almost_equal(ivim_model_LM.estimate_f_D_star(params_f_D,
data_single, S0,
D),
assert_array_almost_equal(ivim_model_trr.estimate_f_D_star(params_f_D,
data_single, S0,
D),
(f, D_star))


Expand Down Expand Up @@ -444,7 +445,7 @@ def test_leastsq_failing():
"""
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always", category=UserWarning)
fit_single = ivim_model_LM.fit(noisy_single)
fit_single = ivim_model_trr.fit(noisy_single)
assert_greater_equal(len(w), 3)
u_warn = [l_w for l_w in w if issubclass(l_w.category, UserWarning)]
assert_greater_equal(len(u_warn), 3)
Expand All @@ -467,7 +468,7 @@ def test_leastsq_error():
"""
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always", category=UserWarning)
fit = ivim_model_LM._leastsq(data_single, [-1, -1, -1, -1])
fit = ivim_model_trr._leastsq(data_single, [-1, -1, -1, -1])
assert_greater_equal(len(w), 1)
assert_(issubclass(w[-1].category, UserWarning))
assert_("" in str(w[-1].message))
Expand Down