Update acquisition scheme and model handling for clinical data (delta…

…s unknown) (#34)

* clinical acquisition schemes from bvals/bvecs can now be generated without delta/Delta when they are unknown.
* acquisition scheme summary can be printed for clinical schemes and N/A is displayed for G/Delta/delta.
* all distributed and compartment models now have _required_acquisition_parameters
* added acquisition scheme checks at any point an acquisition scheme is given in (fitted) modeling frameworks.
* removed superfluous acquisition scheme function.
* updated acquisition scheme tutorial.
* added tests to check acquisition scheme catches for compartment models that need delta/Delta but they are not available.

.travis.yml

@@ -55,7 +55,7 @@ install:
 - python setup.py install --record installed_files.txt
 
 script:
-- flake8 --ignore N802,N806,E731,F401 `find . -name \*.py | grep -v setup.py | grep -v /doc/`
+- flake8 --ignore N802,N806,E731,F401,W504 `find . -name \*.py | grep -v setup.py | grep -v /doc/`
 - nosetests --with-coverage -v
   #- travis-sphinx build -s doc
   #- for a in examples/*ipynb; do 

dmipy/core/fitted_modeling_framework.py

@@ -174,6 +174,9 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
         """
         if acquisition_scheme is None:
             acquisition_scheme = self.model.scheme
+        self.model._check_model_params_with_acquisition_params(
+            acquisition_scheme)
+
         dataset_shape = self.fitted_parameters_vector.shape[:-1]
         if S0 is None:
             S0 = self.S0
@@ -284,6 +287,9 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
         """
         if acquisition_scheme is None:
             acquisition_scheme = self.model.scheme
+        self.model._check_model_params_with_acquisition_params(
+            acquisition_scheme)
+
         dataset_shape = self.fitted_parameters_vector.shape[:-1]
         if S0 is None:
             S0 = self.S0
@@ -737,6 +743,9 @@ def predict(self, acquisition_scheme=None, S0=None, mask=None):
         """
         if acquisition_scheme is None:
             acquisition_scheme = self.model.scheme
+        self.model._check_model_params_with_acquisition_params(
+            acquisition_scheme)
+
         dataset_shape = self.fitted_parameters_vector.shape[:-1]
         if S0 is None:
             if self.fit_S0_response:

dmipy/core/modeling_framework.py

@@ -708,6 +708,16 @@ def _add_optimization_parameter(
         self._inverted_parameter_map.update(
             {(None, 'fraction'): parameter_name})
 
+    def _check_model_params_with_acquisition_params(self, acquisition_scheme):
+        for model in self.models:
+            for parameter in model._required_acquisition_parameters:
+                if getattr(acquisition_scheme, parameter) is None:
+                    msg = "{} is not compatible with ".format(
+                        model.__class__.__name__)
+                    msg += "given acquisition scheme because it needs "
+                    msg += "{} as an acquisition parameter.".format(parameter)
+                    raise ValueError(msg)
+
     def visualize_model_setup(
             self, view=True, cleanup=True, with_parameters=False,
             im_format='png'):
@@ -956,6 +966,7 @@ def fit(self, acquisition_scheme, data,
             functions.
         """
         self._check_tissue_model_acquisition_scheme(acquisition_scheme)
+        self._check_model_params_with_acquisition_params(acquisition_scheme)
 
         # estimate S0
         self.scheme = acquisition_scheme
@@ -1068,6 +1079,8 @@ def simulate_signal(self, acquisition_scheme, parameters_array_or_dict):
             array the same size as x0.
             The simulated signal of the microstructure model.
         """
+        self._check_model_params_with_acquisition_params(acquisition_scheme)
+
         Ndata = acquisition_scheme.number_of_measurements
         if isinstance(parameters_array_or_dict, np.ndarray):
             x0 = parameters_array_or_dict
@@ -1318,6 +1331,7 @@ def fit(self, acquisition_scheme, data,
             functions.
         """
         self._check_tissue_model_acquisition_scheme(acquisition_scheme)
+        self._check_model_params_with_acquisition_params(acquisition_scheme)
 
         # estimate S0
         self.scheme = acquisition_scheme
@@ -1443,6 +1457,8 @@ def simulate_signal(self, acquisition_scheme, parameters_array_or_dict):
             array the same size as x0.
             The simulated signal of the microstructure model.
         """
+        self._check_model_params_with_acquisition_params(acquisition_scheme)
+
         Ndata = acquisition_scheme.shell_indices.max() + 1
         if isinstance(parameters_array_or_dict, np.ndarray):
             x0 = parameters_array_or_dict
@@ -1729,6 +1745,7 @@ def fit(self, acquisition_scheme, data, mask=None, solver='csd',
         """
         self._check_if_kernel_parameters_are_fixed()
         self._check_tissue_model_acquisition_scheme(acquisition_scheme)
+        self._check_model_params_with_acquisition_params(acquisition_scheme)
 
         self.voxel_varying_kernel = False
         if bool(self.x0_parameters):  # if the dictionary is not empty
@@ -1897,6 +1914,8 @@ def simulate_signal(self, acquisition_scheme, parameters_array_or_dict):
             array the same size as x0.
             The simulated signal of the microstructure model.
         """
+        self._check_model_params_with_acquisition_params(acquisition_scheme)
+
         Ndata = acquisition_scheme.number_of_measurements
         if isinstance(parameters_array_or_dict, np.ndarray):
             x0 = parameters_array_or_dict

dmipy/core/tests/test_acquisition_scheme.py

@@ -1,10 +1,17 @@
 import numpy as np
+from dmipy.data.saved_acquisition_schemes import (
+    duval_cat_spinal_cord_2d_acquisition_scheme,
+    wu_minn_hcp_acquisition_scheme)
 from dmipy.core.acquisition_scheme import (
     acquisition_scheme_from_bvalues,
     acquisition_scheme_from_qvalues,
     acquisition_scheme_from_gradient_strengths,
     calculate_shell_bvalues_and_indices,
-    gtab_dipy2mipy, gtab_mipy2dipy)
+    gtab_dipy2dmipy, gtab_dmipy2dipy)
+from dmipy.core.modeling_framework import (
+    MultiCompartmentModel)
+from dmipy.signal_models.cylinder_models import (
+    C4CylinderGaussianPhaseApproximation)
 from dipy.core.gradients import gradient_table
 from numpy.testing import (
     assert_raises, assert_equal, assert_array_equal)
@@ -143,7 +150,7 @@ def test_estimate_shell_indices():
     assert_array_equal(shell_indices, bvalues)
 
 
-def test_shell_indices_with_vayring_diffusion_times(Nsamples=10):
+def test_shell_indices_with_varying_diffusion_times(Nsamples=10):
     # tests whether measurements with the same bvalue but different diffusion
     # time are correctly classified in different shells
     bvalues = np.tile(1e9, Nsamples)
@@ -156,29 +163,29 @@ def test_shell_indices_with_vayring_diffusion_times(Nsamples=10):
     assert_equal(len(np.unique(scheme.shell_indices)), 2)
 
 
-def test_dipy2mipy_acquisition_converter(Nsamples=10):
+def test_dipy2dmipy_acquisition_converter(Nsamples=10):
     bvals = np.tile(1e3, Nsamples)
     bvecs = np.tile(np.r_[1., 0., 0.], (Nsamples, 1))
     big_delta = 0.03
     small_delta = 0.01
     gtab_dipy = gradient_table(
         bvals=bvals, bvecs=bvecs, small_delta=small_delta, big_delta=big_delta)
-    gtab_mipy = gtab_dipy2mipy(gtab_dipy)
+    gtab_mipy = gtab_dipy2dmipy(gtab_dipy)
     assert_array_equal(gtab_mipy.bvalues / 1e6, gtab_dipy.bvals)
     assert_array_equal(gtab_mipy.gradient_directions, gtab_dipy.bvecs)
     assert_equal(np.unique(gtab_mipy.Delta), gtab_dipy.big_delta)
     assert_equal(np.unique(gtab_mipy.delta), gtab_dipy.small_delta)
 
 
-def test_mipy2dipy_acquisition_converter(Nsamples=10):
+def test_dmipy2dipy_acquisition_converter(Nsamples=10):
     bvals = np.tile(1e9, Nsamples)
     bvecs = np.tile(np.r_[1., 0., 0.], (Nsamples, 1))
     big_delta = 0.03
     small_delta = 0.01
     gtab_mipy = acquisition_scheme_from_bvalues(
         bvalues=bvals, gradient_directions=bvecs,
         delta=small_delta, Delta=big_delta)
-    gtab_dipy = gtab_mipy2dipy(gtab_mipy)
+    gtab_dipy = gtab_dmipy2dipy(gtab_mipy)
     assert_array_equal(gtab_mipy.bvalues / 1e6, gtab_dipy.bvals)
     assert_array_equal(gtab_mipy.gradient_directions, gtab_dipy.bvecs)
     assert_equal(gtab_mipy.Delta, gtab_dipy.big_delta)
@@ -194,3 +201,32 @@ def test_acquisition_scheme_summary(Nsamples=10):
         bvalues=bvals, gradient_directions=bvecs,
         delta=small_delta, Delta=big_delta)
     gtab_mipy.print_acquisition_info
+
+
+def test_raise_dmipy2dmipy_multiple_delta_Delta():
+    scheme = duval_cat_spinal_cord_2d_acquisition_scheme()
+    assert_raises(ValueError, gtab_dmipy2dipy, scheme)
+
+
+def test_acquisition_scheme_pruning():
+    scheme = wu_minn_hcp_acquisition_scheme()
+    test_data = np.random.rand(len(scheme.bvalues))
+
+    scheme_pr, data_pr = scheme.return_pruned_acquisition_scheme(
+        [2], test_data)
+    assert_array_equal(
+        scheme_pr.bvalues,
+        scheme.bvalues[scheme.shell_indices == 2])
+    assert_array_equal(
+        data_pr,
+        test_data[scheme.shell_indices == 2])
+
+
+def test_acq_scheme_without_deltas_model_catch():
+    scheme = wu_minn_hcp_acquisition_scheme()
+    test_data = np.random.rand(len(scheme.bvalues))
+    scheme_clinical = acquisition_scheme_from_bvalues(
+        scheme.bvalues, scheme.gradient_directions)
+    mc_model = MultiCompartmentModel(
+        [C4CylinderGaussianPhaseApproximation()])
+    assert_raises(ValueError, mc_model.fit, scheme_clinical, test_data)

dmipy/distributions/distribute_models.py

@@ -356,6 +356,14 @@ def copy(self):
         """
         return copy.copy(self)
 
+    def _set_required_acquisition_parameters(self):
+        self._required_acquisition_parameters = []
+        for model in self.models:
+            self._required_acquisition_parameters += (
+                model._required_acquisition_parameters)
+        self._required_acquisition_parameters = np.unique(
+            self._required_acquisition_parameters).tolist()
+
     @property
     def parameter_names(self):
         "Retuns the DistributedModel parameter names."
@@ -584,6 +592,7 @@ class SD1WatsonDistributed(DistributedModel):
 
     def __init__(self, models, parameter_links=None):
         self.models = models
+        self._set_required_acquisition_parameters()
         self._check_for_double_model_class_instances()
         self._check_for_dispersable_models()
 
@@ -622,6 +631,7 @@ class SD2BinghamDistributed(DistributedModel):
 
     def __init__(self, models, parameter_links=None):
         self.models = models
+        self._set_required_acquisition_parameters()
         self._check_for_double_model_class_instances()
         self._check_for_dispersable_models()
 
@@ -661,6 +671,7 @@ class DD1GammaDistributed(DistributedModel):
     def __init__(self, models, parameter_links=None,
                  target_parameter='diameter'):
         self.models = models
+        self._set_required_acquisition_parameters()
         self.target_parameter = target_parameter
         self._check_for_double_model_class_instances()
         self._check_for_distributable_models()

dmipy/signal_models/capped_cylinder_models.py

@@ -37,7 +37,7 @@ class CC2CappedCylinderStejskalTannerApproximation(ModelProperties):
         spheres and between planes." Journal of Magnetic Resonance, Series A
         104.1 (1993): 17-25.
     """
-
+    _required_acquisition_parameters = ['gradient_directions', 'qvalues']
     _model_type = 'CompartmentModel'
 
     def __init__(
@@ -146,7 +146,8 @@ class CC3CappedCylinderCallaghanApproximation(ModelProperties):
         relaxation." Journal of magnetic resonance, Series A 113.1 (1995):
         53-59.
     """
-
+    _required_acquisition_parameters = [
+        'gradient_directions', 'qvalues', 'tau']
     _model_type = 'CompartmentModel'
 
     def __init__(

dmipy/signal_models/cylinder_models.py

@@ -45,6 +45,8 @@ class C1Stick(ModelProperties):
            Magnetic Resonance in Medicine (2003)
     """
 
+    _required_acquisition_parameters = ['bvalues', 'gradient_directions']
+
     _parameter_ranges = {
         'mu': ([0, np.pi], [-np.pi, np.pi]),
         'lambda_par': (.1, 3)
@@ -182,6 +184,9 @@ class C2CylinderStejskalTannerApproximation(ModelProperties):
             117.1 (1995): 94-97.
     """
 
+    _required_acquisition_parameters = [
+        'bvalues', 'gradient_directions', 'qvalues']
+
     _parameter_ranges = {
         'mu': ([0, np.pi], [-np.pi, np.pi]),
         'lambda_par': (.1, 3),
@@ -346,6 +351,9 @@ class C3CylinderCallaghanApproximation(ModelProperties):
             53-59.
     """
 
+    _required_acquisition_parameters = [
+        'bvalues', 'gradient_directions', 'qvalues', 'tau']
+
     _parameter_ranges = {
         'mu': ([0, np.pi], [-np.pi, np.pi]),
         'lambda_par': (.1, 3),
@@ -436,9 +444,7 @@ def __call__(self, acquisition_scheme, **kwargs):
         bvals = acquisition_scheme.bvalues
         n = acquisition_scheme.gradient_directions
         q = acquisition_scheme.qvalues
-        Delta = acquisition_scheme.Delta
-        delta = acquisition_scheme.delta
-        tau = Delta - delta / 3.0
+        tau = acquisition_scheme.tau
 
         diameter = kwargs.get('diameter', self.diameter)
         lambda_par = kwargs.get('lambda_par', self.lambda_par)
@@ -542,6 +548,10 @@ class C4CylinderGaussianPhaseApproximation(ModelProperties):
             Journal of Magnetic Resonance Series B (1994)
     """
 
+    _required_acquisition_parameters = [
+        'bvalues', 'gradient_directions',
+        'gradient_strengths', 'delta', 'Delta']
+
     _parameter_ranges = {
         'mu': ([0, np.pi], [-np.pi, np.pi]),
         'lambda_par': (.1, 3),

dmipy/signal_models/gaussian_models.py

@@ -38,6 +38,7 @@ class G1Ball(ModelProperties):
             diffusion-weighted MR imaging"
            Magnetic Resonance in Medicine (2003)
     """
+    _required_acquisition_parameters = ['bvalues']
 
     _parameter_ranges = {
         'lambda_iso': (.1, 3)
@@ -151,6 +152,7 @@ class G2Zeppelin(ModelProperties):
            "Compartment models of the diffusion MR signal in brain white
             matter: a taxonomy and comparison". NeuroImage (2012)
     """
+    _required_acquisition_parameters = ['bvalues', 'gradient_directions']
 
     _parameter_ranges = {
         'mu': ([0, np.pi], [-np.pi, np.pi]),
@@ -304,6 +306,8 @@ class G3TemporalZeppelin(ModelProperties):
         structure of neuronal tracts from time-dependent diffusion.
         NeuroImage 114, 18.
     """
+    _required_acquisition_parameters = [
+        'bvalues', 'gradient_directions', 'delta', 'Delta']
 
     _parameter_ranges = {
         'mu': ([0, np.pi], [-np.pi, np.pi]),

dmipy/signal_models/plane_models.py

@@ -24,6 +24,7 @@ class P2PlaneStejskalTannerApproximation(ModelProperties):
         spheres and between planes." Journal of Magnetic Resonance, Series A
         104.1 (1993): 17-25.
     """
+    _required_acquisition_parameters = ['qvalues']
     _parameter_ranges = {
         'diameter': (1e-2, 20)
     }
@@ -92,6 +93,7 @@ class P3PlaneCallaghanApproximation(ModelProperties):
     [1] Callaghan, "Pulsed-Gradient Spin-Echo NMR for Planar, Cylindrical,
         and Spherical Pores under Conditions of Wall Relaxation", JMR 1995
     """
+    _required_acquisition_parameters = ['qvalues', 'tau']
 
     _parameter_ranges = {
         'diameter': (1e-2, 20)

dmipy/signal_models/sphere_models.py

@@ -22,6 +22,7 @@ class S1Dot(ModelProperties):
            "Compartment models of the diffusion MR signal in brain white
             matter: a taxonomy and comparison". NeuroImage (2012)
     """
+    _required_acquisition_parameters = []
 
     _parameter_ranges = {
     }
@@ -119,6 +120,8 @@ class S2SphereStejskalTannerApproximation(ModelProperties):
         spheres and between planes." Journal of Magnetic Resonance, Series A
         104.1 (1993): 17-25.
     """
+    _required_acquisition_parameters = ['qvalues']
+
     _parameter_ranges = {
         'diameter': (1e-2, 20)
     }
@@ -241,6 +244,7 @@ class _S3SphereCallaghanApproximation(ModelProperties):
     [1] Callaghan, "Pulsed-Gradient Spin-Echo NMR for Planar, Cylindrical,
         and Spherical Pores under Conditions of Wall Relaxation", JMR 1995
     """
+    _required_acquisition_parameters = ['qvalues', 'tau']
 
     _parameter_ranges = {
         'diameter': (1e-2, 20)
@@ -333,6 +337,8 @@ class S4SphereGaussianPhaseApproximation(ModelProperties):
         spheres and between planes." Journal of Magnetic Resonance, Series A
         104.1 (1993): 17-25.
     """
+    _required_acquisition_parameters = ['gradient_strengths', 'delta', 'Delta']
+
     _parameter_ranges = {
         'diameter': (1e-2, 20)
     }