<commit_before>import numpy as np
from dipy.reconst.multi_voxel import multi_voxel_fit
from scipy.special import hermite, gamma
from scipy.misc import factorial, factorial2
import dipy.reconst.dti as dti
from warnings import warn
from dipy.core.gradients import gradient_table
from ..utils.optpkg import optional_package

cvxopt, have_cvxopt, _ = optional_package("cvxopt")


class MapmriModel():
    r"""Mean Apparent Propagator MRI (MAPMRI) [1]_ of the diffusion signal.

    The main idea is to model the diffusion signal as a linear combination of
    the continuous functions presented in [2]_ but extending it in three
    dimension.
    The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
    extension is provided using a set of three basis function for the radial
    part, one for the signal along x, one for y and one for z, while [3]_ 
    uses one basis function to model the radial part and real Spherical
    Harmonics to model the angular part.
    From the MAPMRI coefficients is possible to use the analytical formulae 
    to estimate the ODF.

    References
    ----------
    .. [1] <NAME>. et. al, "Mean apparent propagator (MAP) MRI: A novel
           diffusion imaging method for mapping tissue microstructure",
           NeuroImage, 2013.

    .. [2] <NAME>. et. al, "Simple harmonic oscillator based reconstruction
           and estimation for one-dimensional q-space magnetic resonance
           1D-SHORE)", eapoc Intl Soc Mag Reson Med, vol. 16, p. 35., 2008.

    .. [3] <NAME>. et. al, "Continuous diffusion signal, EAP and ODF
           estimation via Compressive Sensing in diffusion MRI", Medical
           Image Analysis, 2013.
    """

    def __init__(self, gtab, radial_order=4, lambd=1e-09,  eap_cons = False, anisotropic_scaling = True):
        r""" Analytical and continuous modeling of the diffusion signal with
        respect to the MAPMRI basis [1]_.

        The main idea is to model the diffusion signal as a linear combination of
        the continuous functions presented in [2]_ but extending it in three
        dimension.

        The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
        extension is provided using a set of three basis function for the radial
        part, one for the signal along x, one for y and one for z, while [3]_ 
        uses one basis function to model the radial part and real Spherical
        Harmonics to model the angular part.
        From the MAPMRI coefficients is possible to use the analytical formulae 
        to estimate the ODF.


        Parameters
        ----------
        gtab : GradientTable,
            gradient directions and bvalues container class
        radial_order : unsigned int,
            an even integer that represent the order of the basis
        lambd : float,
            radial regularisation constant
        eap_cons : bool,
            Constrain the propagator to be positive.
        anisotropic_scaling : bool,
            If false, force the basis function to be identical in the three
            dimensions (SHORE like).

        References
        ----------
        .. [1] <NAME>. et. al, "Mean apparent propagator (MAP) MRI: A novel
               diffusion imaging method for mapping tissue microstructure",
               NeuroImage, 2013.

        .. [2] <NAME>. et. al, "Simple harmonic oscillator based reconstruction
               and estimation for one-dimensional q-space magnetic resonance
               1D-SHORE)", eapoc Intl Soc Mag Reson Med, vol. 16, p. 35., 2008.

        .. [3] <NAME>. et. al, "Continuous diffusion signal, EAP and ODF
               estimation via Compressive Sensing in diffusion MRI", Medical
               Image Analysis, 2013.

        Examples
        --------
        In this example, where the data, gradient table and sphere tessellation
        used for reconstruction are provided, we model the diffusion signal
        with respect to the MAPMRI model and compute the real and analytical
        ODF.

        from dipy.data import get_data,get_sphere
        sphere = get_sphere('symmetric724')
        fimg, fbvals, fbvecs = get_data('ISBI_testing_2shells_table')
        bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
        gtab = gradient_table(bvals, bvecs)
        from dipy.sims.voxel import SticksAndBall
        data, golden_directions = SticksAndBall(gtab, d=0.0015,
                                                S0=1, angles=[(0, 0), (90, 0)],
                                                fractions=[50, 50], snr=None)
        from dipy.reconst.mapmri import MapmriModel
        radial_order = 4
        map_model = MapmriModel(gtab, radial_order=radial_order)
        mapfit = map_model.fit(data)
        odf= mapfit.odf(sphere)
        """

        self.bvals = gtab.bvals
        self.bvecs = gtab.bvecs
        self.gtab = gtab
        self.radial_order = radial_order
        self.lambd = lambd
        self.eap_cons = eap_cons
        self.anisotropic_scaling=anisotropic_scaling
        if (gtab.big_delta is None) or (gtab.small_delta is None):
            self.tau = 1 / (4 * np.pi ** 2)
        else:
            self.tau = gtab.big_delta - gtab.small_delta / 3.0


    @multi_voxel_fit
    def fit(self, data):
        ind=self.gtab.bvals<=2000
        gtab2 = gradient_table(self.gtab.bvals[ind], self.gtab.bvecs[ind,:])
        tenmodel=dti.TensorModel(gtab2)
        tenfit=tenmodel.fit(data[...,ind])
        
        evals = tenfit.evals
        R=tenfit.evecs

        ind_evals = np.argsort(evals)[::-1]
        evals = evals[ind_evals]
        R = R[ind_evals,:]
        
        evals = np.clip(evals,1e-04,evals.max())
        
        if self.anisotropic_scaling:
            mu = np.sqrt(evals*2*self.tau)
            
        else:
            mumean=np.sqrt(evals.mean()*2*self.tau)
            mu=np.array([mumean,mumean,mumean])
        
        qvals=np.sqrt(self.gtab.bvals/self.tau) / (2 * np.pi)
        qvecs=np.dot(self.gtab.bvecs,R)
        q=qvecs*qvals[:,None]

        M = mapmri_phi_matrix(self.radial_order, mu, q.T)
        
        ind_mat = mapmri_index_matrix(self.radial_order)

        # This is a simple empirical regularization, to be replaced
        I = np.diag(ind_mat.sum(1)**2)

        if self.eap_cons:
            if not have_cvxopt:
                raise ValueError(
                    'CVXOPT package needed to enforce constraints')
            import cvxopt.solvers
            rmax = 2* np.sqrt(10 * evals.max()*self.tau)
            r_index, r_grad = create_rspace(11, rmax)
            K = mapmri_psi_matrix(self.radial_order,  mu, r_grad[0:len(r_grad)/2,:])<commit_msg>fix regularization with identity matrix<commit_after>import numpy as np
from dipy.reconst.multi_voxel import multi_voxel_fit
from scipy.special import hermite, gamma
from scipy.misc import factorial, factorial2
import dipy.reconst.dti as dti
from warnings import warn
from dipy.core.gradients import gradient_table
from ..utils.optpkg import optional_package

cvxopt, have_cvxopt, _ = optional_package("cvxopt")


class MapmriModel():
    r"""Mean Apparent Propagator MRI (MAPMRI) [1]_ of the diffusion signal.

    The main idea is to model the diffusion signal as a linear combination of
    the continuous functions presented in [2]_ but extending it in three
    dimension.
    The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
    extension is provided using a set of three basis function for the radial
    part, one for the signal along x, one for y and one for z, while [3]_ 
    uses one basis function to model the radial part and real Spherical
    Harmonics to model the angular part.
    From the MAPMRI coefficients is possible to use the analytical formulae 
    to estimate the ODF.

    References
    ----------
    .. [1] <NAME>. al, "Mean apparent propagator (MAP) MRI: A novel
           diffusion imaging method for mapping tissue microstructure",
           NeuroImage, 2013.

    .. [2] <NAME>. al, "Simple harmonic oscillator based reconstruction
           and estimation for one-dimensional q-space magnetic resonance
           1D-SHORE)", eapoc Intl Soc Mag Reson Med, vol. 16, p. 35., 2008.

    .. [3] <NAME>. al, "Continuous diffusion signal, EAP and ODF
           estimation via Compressive Sensing in diffusion MRI", Medical
           Image Analysis, 2013.
    """

    def __init__(self,
                 gtab,
                 radial_order=4,
                 lambd=1e-16,
                 eap_cons = False,
                 anisotropic_scaling = True,
                 eigenvalue_threshold = 1e-04,
                 bmax_threshold = 2000):
        r""" Analytical and continuous modeling of the diffusion signal with
        respect to the MAPMRI basis [1]_.

        The main idea is to model the diffusion signal as a linear combination of
        the continuous functions presented in [2]_ but extending it in three
        dimension.

        The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
        extension is provided using a set of three basis function for the radial
        part, one for the signal along x, one for y and one for z, while [3]_ 
        uses one basis function to model the radial part and real Spherical
        Harmonics to model the angular part.
        From the MAPMRI coefficients is possible to use the analytical formulae 
        to estimate the ODF.


        Parameters
        ----------
        gtab : GradientTable,
            gradient directions and bvalues container class
        radial_order : unsigned int,
            an even integer that represent the order of the basis
        lambd : float,
            radial regularisation constant
        eap_cons : bool,
            Constrain the propagator to be positive.
        anisotropic_scaling : bool,
            If false, force the basis function to be identical in the three
            dimensions (SHORE like).
        eigenvalue_threshold : float,
            set the minimum of the tensor eigenvalues in order to avoid 
            stability problem
        bmax_threshold : float,
            set the maximum b-value for the tensor estimation

        References
        ----------
        .. [1] <NAME>. al, "Mean apparent propagator (MAP) MRI: A novel
               diffusion imaging method for mapping tissue microstructure",
               NeuroImage, 2013.

        .. [2] <NAME>. et. al, "Simple harmonic oscillator based reconstruction
               and estimation for one-dimensional q-space magnetic resonance
               1D-SHORE)", eapoc Intl Soc Mag Reson Med, vol. 16, p. 35., 2008.

        .. [3] <NAME>. al, "Simple harmonic oscillator based reconstruction
               and estimation for three-dimensional q-space mri", ISMRM 2009.

        Examples
        --------
        In this example, where the data, gradient table and sphere tessellation
        used for reconstruction are provided, we model the diffusion signal
        with respect to the MAPMRI model and compute the real and analytical
        ODF.

        from dipy.data import get_data,get_sphere
        sphere = get_sphere('symmetric724')
        fimg, fbvals, fbvecs = get_data('ISBI_testing_2shells_table')
        bvals, bvecs = read_bvals_bvecs(fbvals, fbvecs)
        gtab = gradient_table(bvals, bvecs)
        from dipy.sims.voxel import SticksAndBall
        data, golden_directions = SticksAndBall(gtab, d=0.0015,
                                                S0=1, angles=[(0, 0), (90, 0)],
                                                fractions=[50, 50], snr=None)
        from dipy.reconst.mapmri import MapmriModel
        radial_order = 4
        map_model = MapmriModel(gtab, radial_order=radial_order)
        mapfit = map_model.fit(data)
        odf= mapfit.odf(sphere)
        """

        self.bvals = gtab.bvals
        self.bvecs = gtab.bvecs
        self.gtab = gtab
        self.radial_order = radial_order
        self.lambd = lambd
        self.eap_cons = eap_cons
        self.anisotropic_scaling=anisotropic_scaling
        if (gtab.big_delta is None) or (gtab.small_delta is None):
            self.tau = 1 / (4 * np.pi ** 2)
        else:
            self.tau = gtab.big_delta - gtab.small_delta / 3.0
        self.eigenvalue_threshold = eigenvalue_threshold
        self.bmax_threshold = bmax_threshold

    @multi_voxel_fit
    def fit(self, data):
        ind=self.gtab.bvals<=self.bmax_threshold
        gtab2 = gradient_table(self.gtab.bvals[ind], self.gtab.bvecs[ind,:])
        tenmodel=dti.TensorModel(gtab2)
        tenfit=tenmodel.fit(data[...,ind])
        
        evals = tenfit.evals
        R=tenfit.evecs

        ind_evals = np.argsort(evals)[::-1]
        evals = evals[ind_evals]
        R = R[ind_evals,:]
        
        evals = np.clip(evals,self.eigenvalue_threshold,evals.max())
        
        if self.anisotropic_scaling:
            mu = np.sqrt(evals*2*self.tau)
            
        else:
            mumean=np.sqrt(evals.mean()*2*self.tau)
            mu=np.array([mumean,mumean,mumean])
        
        qvals=np.sqrt(self.gtab.bvals/self.tau) / (2 * np.pi)
        qvecs=np.dot(self.gtab.bvecs,R)
        q=qvecs*qvals[:,None]

        M = mapmri_phi_matrix(self.radial_order, mu, q.T)
        
        ind_mat = mapmri_index_matrix(self.radial_order)

        # This is a simple empirical regularization, to be replaced
        I = np.eye(ind_mat.shape[0])

        if self.eap_cons:
            if not have_cvxopt:
                raise ValueError(
                    'CVXOPT package needed to enforce constraints')
            import cvxopt.solvers
            rmax = 2 * np.sqrt(10 * evals.max()*self.tau)
            r_index, r_grad = create_rspace(11, rmax)
            K = mapmri_psi_matrix(self.radial_order,  mu, r_grad[0:len(r_grad)/2,:])
{'fraction_non_alphanumeric': 0.07857554786620531, 'fraction_numerical': 0.014777970011534025, 'mean_word_length': 2.122439792932703, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6306537', 'n_tokens_mistral': 4266, 'n_tokens_neox': 3943, 'n_words': 1341}