MCMDI issues when b0s are not zeros (b=0)

[villalonreina]

Hi @matteofrigo @rutgerfick We have been trying to run MCMDI models on UKBB data and we were getting the following error:

File "/Users/jvillalo/opt/anaconda3/lib/python3.8/site-packages/dmipy/signal_models/cylinder_models.py", line 121, in spherical_mean
    E_mean[~acquisition_scheme.shell_b0_mask] = E_mean_
ValueError: NumPy boolean array indexing assignment cannot assign 3 input values to the 2 output values where the mask is true

After analyzing the error, we found out that the b0s in UKBB are not full b0s, but instead they are b=5.
Hence, there is an issue in dmipy/signal_models/cylinder_models.py because of the way it indexes bvals greater than zero (bval > 0). See below:

        bvals = acquisition_scheme.shell_bvalues
        bvals_ = bvals[~acquisition_scheme.shell_b0_mask]

        lambda_par = kwargs.get('lambda_par', self.lambda_par)

        E_mean = np.ones_like(bvals)
        bval_indices_above0 = bvals > 0
        bvals_ = bvals[bval_indices_above0]
        E_mean_ = ((np.sqrt(np.pi) * erf(np.sqrt(bvals_ * lambda_par))) /
                   (2 * np.sqrt(bvals_ * lambda_par)))
        E_mean[~acquisition_scheme.shell_b0_mask] = E_mean_
        return E_mean

[rutgerfick]

Hey @villalonreina,

When you create the dmipy gradient table you can indicate what is the b-value threshold for counting a measurement as a b0. You can avoid this error in the future this way :-)

https://github.com/AthenaEPI/dmipy/blob/dcb0916cbc052636706af6aeee0e3ad105c68b8f/dmipy/core/acquisition_scheme.py#L494

Let me know if it solves the problem!

[villalonreina]

Thanks @rutgerfick. Yes, we have this line in the code already:

scheme = acquisition_scheme_from_bvalues(bvals*1e6, bvecs, small_delta,
                                         big_delta, b0_threshold=50e6,
                                         min_b_shell_distance=50e6)

Our b-zeros are actually b=5, and the code still crashes, unless we modify the bval file. I tracked the error and it took me to the dmipy/signal_models/cylinder_models.py as noted above.

[rutgerfick]

You're right this is a bug.

[rutgerfick]

merging the fix tomorrow morning, it's really small.

[rutgerfick]

Hey @villalonreina , please update to the latest version of dmipy and let me know if this solves this particular problem.