DIPY MSMT Fit Not Converging

,
1

Hello,

I am trying to migrate my MrTrix workflow to DIPY (1.4.1). I have multishell data (1 average b0, 64 directions at b=1000 and 2000) preprocessed with QSIPrep that I want to run MSMT CSD on. The only thing that has happened with the data between preprocessing and this script is running MRtrix dwigradcheck to make sure the gradient table is correct. I am working with the tutorial linked here. My code is below.

# DEFINE PATHS TO DATA
dwi_fodf_path ='/om4/group/gablab/data/HBN/derivatives/sbci/sub-NDARAA306NT2/dwi/dwi.nii.gz'
bvals_path = '/om4/group/gablab/data/HBN/derivatives/sbci/sub-NDARAA306NT2/dwi/dwi.bval'
bvecs_path = '/om4/group/gablab/data/HBN/derivatives/sbci/sub-NDARAA306NT2/dwi/dwi.bvec'

# LOAD DATA
data, affine = load_nifti(dwi_fodf_path)
bvals, bvecs = read_bvals_bvecs(bvals_path, bvecs_path)
bvals=bvals.round(1) # round b-vals to whole numbers for easier indexing
gtab = gradient_table(bvals,bvecs)

# INDEX B-VALS < 1200 FOR BETTER TENSOR FITTING
low_b = np.logical_or(bvals == 0, bvals == 1000)
data_low_b = data[..., low_b]
gtab_low_b = gradient_table(bvals[low_b], bvecs[low_b])

# GET TISSUE MASKS
mask_wm, mask_gm, mask_csf = mask_for_response_msmt(gtab_low_b, data_low_b, roi_radii=10,
                                                    wm_fa_thr=0.7,
                                                    gm_fa_thr=0.3,
                                                    csf_fa_thr=0.15,
                                                    gm_md_thr=0.001,
                                                    csf_md_thr=0.0032)

# GET TISSUE RESPONSES
response_wm, response_gm, response_csf = response_from_mask_msmt(gtab, data,
                                                                 mask_wm,
                                                                 mask_gm,
                                                                 mask_csf)
# GET MSMT FIBNER RESPONSE
ubvals = unique_bvals_tolerance(gtab.bvals)
response_mcsd = multi_shell_fiber_response(sh_order=8,
                                           bvals=ubvals,
                                           wm_rf=response_wm,
                                           gm_rf=response_gm,
                                           csf_rf=response_csf)
response = np.array([response_wm, response_gm, response_csf])

# DEFINE AND FIT MULTI SHELL SPHERICAL DECONVOLUTION MODEL
mcsd_model_simple_response = MultiShellDeconvModel(gtab, response, sh_order=8)
mcsd_model = MultiShellDeconvModel(gtab, response_mcsd)
mcsd_fit = mcsd_model.fit(data)

I get the following error message at the last step of the process:
---------------------------------------------------------------------------
ArpackNoConvergence Traceback (most recent call last)
/om2/user/smeisler/anaconda3/envs/nipype/lib/python3.8/site-packages/cvxpy/expressions/constants/constant.py in is_psd_within_tol(A, tol)
260 try:
→ 261 ev = SA_eigsh(-tol) # might return np.NaN, or raise exception
262 finally:

    /om2/user/smeisler/anaconda3/envs/nipype/lib/python3.8/site-packages/cvxpy/expressions/constants/constant.py in SA_eigsh(sigma)
        244     def SA_eigsh(sigma):
    --> 245         return eigsh(A, k=1, sigma=sigma, which='SA', return_eigenvectors=False)
        246         # Returns the eigenvalue w[i] of A where 1/(w[i] - sigma) is minimized.

    /om2/user/smeisler/anaconda3/envs/nipype/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py in eigsh(A, k, M, sigma, which, v0, ncv, maxiter, tol, return_eigenvectors, Minv, OPinv, mode)
       1689         while not params.converged:
    -> 1690             params.iterate()
       1691 

    /om2/user/smeisler/anaconda3/envs/nipype/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py in iterate(self)
        569             elif self.info == 1:
    --> 570                 self._raise_no_convergence()
        571             else:

    /om2/user/smeisler/anaconda3/envs/nipype/lib/python3.8/site-packages/scipy/sparse/linalg/eigen/arpack/arpack.py in _raise_no_convergence(self)
        375             k_ok = 0
    --> 376         raise ArpackNoConvergence(msg % (num_iter, k_ok, self.k), ev, vec)
        377 

    ArpackNoConvergence: ARPACK error -1: No convergence (471 iterations, 0/1 eigenvectors converged)

MSMT CSD has worked on MrTrix for this subject. Any suggestions on what I can do here?

Thanks,
Steven

2

Hi Steven,

I can help you with this! Can you tell me your CVXPY version? This typically has happened in the past when CVXPY changed something in their system causing the quadratic programming solver within MSMT to not converge. The version of CVXPY for which we’ve tested is: 1.1.10!

Quick side question: are you using a mask? Sometimes having 0s or NaNs (outside the brain) can affect CVXPY convergence.

Thanks,
Shreyas

3

Ah, thanks for the heads up! Current version I have is 1.1.15, I can try to roll back. The DWI data has that is being loaded in the first few lines has not been masked.

Best,
Steven

1 Like
4

Perfect! Let us know if you need any further help :slight_smile: I would suggest also using Gitter and GitHub discussions of DIPY. You’ll get almost instant responses on these forums and would help other users running into similar issues!

Thanks,
Shreyas

1 Like