<commit_before># -*- coding: utf-8 -*-
import numpy as np
from dipy.reconst.multi_voxel import multi_voxel_fit
from dipy.reconst.base import ReconstModel, ReconstFit
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(ReconstModel):

    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
    dimensions.
    The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
    extension is provided using a set of three basis functions 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.

    .. [4] Fick et al. "Laplacian-Regularized MAP-MRI: Improving Axonal
           Caliber", ISBI, 2015.

    .. [5] <NAME>., 2014. Estimation and Processing of Ensemble Average
           Propagator and Its Features in Diffusion MRI. Ph.D. Thesis.

    .. [6] Hosseinbor et al. "Bessel fourier orientation reconstruction
           (bfor): An analytical diffusion propagator reconstruction for hybrid
           diffusion imaging and computation of q-space indices. NeuroImage 64,
           2013, 650–670.

    .. [7] Craven et al. "Smoothing Noisy Data with Spline Functions."
           NUMER MATH 31.4 (1978): 377-403.
    """

    def __init__(self,
                 gtab,
                 radial_order=4,
                 laplacian_regularization=True,
                 laplacian_weighting='GCV',
                 positivity_constraint=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
        dimensions.

        The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
        extension is provided using a set of three basis functions 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
        laplacian_regularization: bool,
            Regularize using the Laplacian of the MAP-MRI basis.
        laplacian_weighting: string or scalar or ndarray,
            The string 'GCV' makes it use generalized cross-validation to find
            the regularization weight [3]. A scalar sets the regularization
            weight to that value. A numpy array of values will make the GCV
            function find the optimal value among those values.
        positivity_constraint : 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

        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>. et. 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 analytical ODF.

        >>> from dipy.core.gradients import gradient_table
        >>> from dipy.data import dsi_voxels, get_sphere
        >>> data, gtab = dsi_voxels()
        >>> sphere = get_sphere('symmetric724')
        >>> from dipy.sims.voxel import SticksAndBall
        >>> data, golden_directions = SticksAndBall(gtab, d=0.0015, S0=100,
        ... 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.laplacian_regularization = laplacian_regularization
        if laplacian_weighting == 'GCV' or \
                np.isscalar(laplacian_weighting) or \
                type(laplacian_weighting) == np.ndarray:
            if np.isscalar(laplacian_weighting) and \
                    laplacian_weighting < 0.:
                msg = "Laplacian Regularization weighting must be positive."
                raise ValueError(msg)
            if type(laplacian_weighting) == np.ndarray and \
                    laplacian_weighting.min() < 0.:
                msg = "Laplacian Regularization weighting must be positive."
                raise ValueError(msg)
            else:
                self.laplacian_weighting = laplacian_weighting

        self.R_mat, self.L_mat, self.S_mat = mapmri_RLS_regularization_matrices(
            radial_order)

        self.positivity_constraint = positivity_constraint
        if self.positivity_constraint:
            if not have_cvxopt:
                raise ValueError(
                    'CVXOPT package needed to enforce constraints')
            import cvxopt.solvers

        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.tenmodel = dti.TensorModel(gtab)
        self.ind_mat = mapmri_index_matrix(self.radial_order)
        self.Bm = b_mat(self.ind_mat)

    @multi_voxel_fit
    def fit(self, data):

        tenfit = self.tenmodel.fit(data)
        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)

        if self.laplacian_regularization:
            laplacian_matrix = mapmri_laplacian_reg_matrix(
                self.ind_mat, mu, self.R_mat, self.L_mat, self.S_mat)
            if self.laplacian_weighting == 'GCV':
                lopt = generalized_crossvalidation(data, M, laplacian_matrix)
            elif np.isscalar(self.laplacian_weighting):
                lopt = self.laplacian_weighting
            elif type(self.laplacian_weighting) == np.ndarray:
                lopt = generalized_crossvalidation(data, M, laplacian_matrix,
                                                   self.laplacian_weighting)
        else:
            lopt = 0.
            laplacian_matrix = np.ones((self.ind_mat.shape[0],
                                        self.ind_mat.shape[0]))

        if self.positivity_constraint:
            if not have_cvxopt:
                raise ValueError(<commit_msg>Removed reordering of eigenvalues in the scale factor estimation<commit_after># -*- coding: utf-8 -*-
import numpy as np
from dipy.reconst.multi_voxel import multi_voxel_fit
from dipy.reconst.base import ReconstModel, ReconstFit
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(ReconstModel):

    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
    dimensions.
    The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
    extension is provided using a set of three basis functions 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.

    .. [4] Fick et al. "Laplacian-Regularized MAP-MRI: Improving Axonal
           Caliber", ISBI, 2015.

    .. [5] <NAME>., 2014. Estimation and Processing of Ensemble Average
           Propagator and Its Features in Diffusion MRI. Ph.D. Thesis.

    .. [6] Hosseinbor et al. "Bessel fourier orientation reconstruction
           (bfor): An analytical diffusion propagator reconstruction for hybrid
           diffusion imaging and computation of q-space indices. NeuroImage 64,
           2013, 650–670.

    .. [7] Craven et al. "Smoothing Noisy Data with Spline Functions."
           NUMER MATH 31.4 (1978): 377-403.
    """

    def __init__(self,
                 gtab,
                 radial_order=4,
                 laplacian_regularization=True,
                 laplacian_weighting='GCV',
                 positivity_constraint=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
        dimensions.

        The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
        extension is provided using a set of three basis functions 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
        laplacian_regularization: bool,
            Regularize using the Laplacian of the MAP-MRI basis.
        laplacian_weighting: string or scalar or ndarray,
            The string 'GCV' makes it use generalized cross-validation to find
            the regularization weight [3]. A scalar sets the regularization
            weight to that value. A numpy array of values will make the GCV
            function find the optimal value among those values.
        positivity_constraint : 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

        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>. et. 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 analytical ODF.

        >>> from dipy.core.gradients import gradient_table
        >>> from dipy.data import dsi_voxels, get_sphere
        >>> data, gtab = dsi_voxels()
        >>> sphere = get_sphere('symmetric724')
        >>> from dipy.sims.voxel import SticksAndBall
        >>> data, golden_directions = SticksAndBall(gtab, d=0.0015, S0=100,
        ... 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.laplacian_regularization = laplacian_regularization
        if laplacian_weighting == 'GCV' or \
                np.isscalar(laplacian_weighting) or \
                type(laplacian_weighting) == np.ndarray:
            if np.isscalar(laplacian_weighting) and \
                    laplacian_weighting < 0.:
                msg = "Laplacian Regularization weighting must be positive."
                raise ValueError(msg)
            if type(laplacian_weighting) == np.ndarray and \
                    laplacian_weighting.min() < 0.:
                msg = "Laplacian Regularization weighting must be positive."
                raise ValueError(msg)
            else:
                self.laplacian_weighting = laplacian_weighting

        self.R_mat, self.L_mat, self.S_mat = mapmri_RLS_regularization_matrices(
            radial_order)

        self.positivity_constraint = positivity_constraint
        if self.positivity_constraint:
            if not have_cvxopt:
                raise ValueError(
                    'CVXOPT package needed to enforce constraints')
            import cvxopt.solvers

        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.tenmodel = dti.TensorModel(gtab)
        self.ind_mat = mapmri_index_matrix(self.radial_order)
        self.Bm = b_mat(self.ind_mat)

    @multi_voxel_fit
    def fit(self, data):

        tenfit = self.tenmodel.fit(data)
        evals = tenfit.evals
        R = tenfit.evecs
        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)

        if self.laplacian_regularization:
            laplacian_matrix = mapmri_laplacian_reg_matrix(
                self.ind_mat, mu, self.R_mat, self.L_mat, self.S_mat)
            if self.laplacian_weighting == 'GCV':
                lopt = generalized_crossvalidation(data, M, laplacian_matrix)
            elif np.isscalar(self.laplacian_weighting):
                lopt = self.laplacian_weighting
            elif type(self.laplacian_weighting) == np.ndarray:
                lopt = generalized_crossvalidation(data, M, laplacian_matrix,
                                                   self.laplacian_weighting)
        else:
            lopt = 0.
            laplacian_matrix = np.ones((self.ind_mat.shape[0],
                                        self.ind_mat.shape[0]))

        if self.positivity_constraint:
            if not have_cvxopt:
                raise ValueError(
{'fraction_non_alphanumeric': 0.07574793125397836, 'fraction_numerical': 0.014162953532781668, 'mean_word_length': 2.124979280623239, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 83, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5863079', 'n_tokens_mistral': 5642, 'n_tokens_neox': 5065, 'n_words': 1755}