<commit_before>""" Classes and functions for fitting ivim model """
from __future__ import division, print_function, absolute_import
import numpy as np

from dipy.core.gradients import gradient_table
from dipy.reconst.base import ReconstModel
from dipy.reconst.dti import apparent_diffusion_coef
from dipy.reconst.dti import TensorModel, mean_diffusivity
from dipy.reconst.multi_voxel import multi_voxel_fit
from dipy.reconst.vec_val_sum import vec_val_vect
from dipy.core.sphere import Sphere

from distutils.version import LooseVersion
import scipy

SCIPY_LESS_0_17 = LooseVersion(scipy.version.short_version) < '0.17'

if not SCIPY_LESS_0_17:
    least_squares = scipy.optimize.least_squares
else:
    leastsq = scipy.optimize.leastsq


def ivim_prediction(params, gtab, S0=1.):
    """The Intravoxel incoherent motion (IVIM) model function.

    Parameters
    ----------
    params : array
        An array of IVIM parameters - S0, f, D_star, D

    gtab : GradientTable class instance
        Gradient directions and bvalues

    S0 : float, optional
        This has been added just for consistency with the existing
        API. Unlike other models, IVIM predicts S0 and this is over written
        by the S0 value in params.

    References
    ----------
    .. [1] <NAME>, Denis, et al. "Separation of diffusion
               and perfusion in intravoxel incoherent motion MR
               imaging." Radiology 168.2 (1988): 497-505.
    .. [2] Federau, Christian, et al. "Quantitative measurement
               of brain perfusion with intravoxel incoherent motion
               MR imaging." Radiology 265.3 (2012): 874-881.
    """
    S0, f, D_star, D = params
    b = gtab.bvals
    S = S0 * (f * np.exp(-b * D_star) + (1 - f) * np.exp(-b * D))
    return S


def _ivim_error(params, gtab, signal):
    """Error function to be used in fitting the IVIM model

    Parameters
    ----------
    params : array
        An array of IVIM parameters. [S0, f, D_star, D]

    gtab : GradientTable class instance
        Gradient directions and bvalues.

    signal : array
        Array containing the actual signal values.

    """
    return (signal - ivim_prediction(params, gtab))


class IvimModel(ReconstModel):
    """Ivim model
    """

    def __init__(self, gtab, split_b=200.0,
                 bounds=None, tol=1e-7,
                 options={'gtol': 1e-7, 'ftol': 1e-7,
                          'eps': 1e-7, 'maxiter': 1000}):
        """
        Initialize an IVIM model.

        The IVIM model assumes that biological tissue includes a volume
        fraction 'f' of water flowing in perfused capillaries, with a
        perfusion coefficient D* and a fraction (1-f) of static (diffusion
        only), intra and extracellular water, with a diffusion coefficient
        D. In this model the echo attenuation of a signal in a single voxel
        can be written as

            .. math::

            S(b) = S_0[f*e^{(-b*D\*)} + (1-f)e^{(-b*D)}]

            Where:
            .. math::

            S_0, f, D\* and D are the IVIM parameters.

        Parameters
        ----------
        gtab : GradientTable class instance
            Gradient directions and bvalues

        split_b : float, optional
            The b-value to split the data on for two-stage fit.
            default : 200.

        bounds : tuple of arrays with 4 elements, optional
            Bounds to constrain the fitted model parameters. This is only supported for
            Scipy version > 0.17. When using a older scipy version, this function will raise
            an error if bounds are different from None.
            default : ([0., 0., 0., 0.], [np.inf, 1., 1., 1.])

        tol : float, optional
            Tolerance for convergence of minimization.
            default : 1e-7

        options : dict, optional
            Dictionary containing gtol, ftol, eps and maxiter. This is passed
            to leastsq.
            default : options={'gtol': 1e-7, 'ftol': 1e-7, 'eps': 1e-7, 'maxiter': 1000}

        References
        ----------
        .. [1] <NAME>, Denis, et al. "Separation of diffusion
                   and perfusion in intravoxel incoherent motion MR
                   imaging." Radiology 168.2 (1988): 497-505.
        .. [2] Federau, Christian, et al. "Quantitative measurement
                   of brain perfusion with intravoxel incoherent motion
                   MR imaging." Radiology 265.3 (2012): 874-881.
        """
        ReconstModel.__init__(self, gtab)
        self.split_b = split_b
        self.bounds = bounds
        self.tol = tol
        self.options = options

        if SCIPY_LESS_0_17 and self.bounds is not None:
            e_s = "Scipy versions less than 0.17 do not support "
            e_s += "bounds. Please update to Scipy 0.17 to use bounds"
            raise ValueError(e_s)
        else:
            self.bounds = (np.array([0., 0., 0., 0.]),
                           np.array([np.inf, 1., 0.1, 0.1]))

    def estimate_x0(self, data):
        """
        Fit the ivim params using a two stage fit.

        In the two stage fitting routine, initially, we fit the signal
        at bvals higher than the specified `split_b` using `TensorModel`
        and get an intial guess for f and D. Then, using these parameters
        we fit the entire data for all bvalues.

        Parameters
        ----------
        data : array
            The measured signal from one voxel. A multi voxel decorator
            will be applied to this fit method to scale it and apply it
            to multiple voxels.

        Returns
        -------
        x0_guess : array
            An array with initial values of S0, f, D_star, D for each voxel.
        """
        S0_hat, D_guess = self._estimate_S0_D(data)
        f_guess = 1 - S0_hat / data[0]
        # We set the S0 guess as the signal value at b=0
        # The D* guess is roughly 10 times the D value
        x0 = np.array([data[0], f_guess, 10 * D_guess, D_guess])
        # The API does not allow bounds for Scipy < 0.17. While estimating x0,
        # if there is noise in the data the estimated x0 might not be feasible.
        # In such a case we will use the values for the parameters<commit_msg>replaced data[0] with mean of data[b0_mask]<commit_after>""" Classes and functions for fitting ivim model """
from __future__ import division, print_function, absolute_import
import numpy as np

from dipy.core.gradients import gradient_table
from dipy.reconst.base import ReconstModel
from dipy.reconst.dti import apparent_diffusion_coef
from dipy.reconst.dti import TensorModel, mean_diffusivity
from dipy.reconst.multi_voxel import multi_voxel_fit
from dipy.reconst.vec_val_sum import vec_val_vect
from dipy.core.sphere import Sphere

from distutils.version import LooseVersion
import scipy

SCIPY_LESS_0_17 = LooseVersion(scipy.version.short_version) < '0.17'

if not SCIPY_LESS_0_17:
    least_squares = scipy.optimize.least_squares
else:
    leastsq = scipy.optimize.leastsq


def ivim_prediction(params, gtab, S0=1.):
    """The Intravoxel incoherent motion (IVIM) model function.

    Parameters
    ----------
    params : array
        An array of IVIM parameters - S0, f, D_star, D

    gtab : GradientTable class instance
        Gradient directions and bvalues

    S0 : float, optional
        This has been added just for consistency with the existing
        API. Unlike other models, IVIM predicts S0 and this is over written
        by the S0 value in params.

    References
    ----------
    .. [1] <NAME>, Denis, et al. "Separation of diffusion
               and perfusion in intravoxel incoherent motion MR
               imaging." Radiology 168.2 (1988): 497-505.
    .. [2] Federau, Christian, et al. "Quantitative measurement
               of brain perfusion with intravoxel incoherent motion
               MR imaging." Radiology 265.3 (2012): 874-881.
    """
    S0, f, D_star, D = params
    b = gtab.bvals
    S = S0 * (f * np.exp(-b * D_star) + (1 - f) * np.exp(-b * D))
    return S


def _ivim_error(params, gtab, signal):
    """Error function to be used in fitting the IVIM model

    Parameters
    ----------
    params : array
        An array of IVIM parameters. [S0, f, D_star, D]

    gtab : GradientTable class instance
        Gradient directions and bvalues.

    signal : array
        Array containing the actual signal values.

    """
    return (signal - ivim_prediction(params, gtab))


class IvimModel(ReconstModel):
    """Ivim model
    """

    def __init__(self, gtab, split_b=200.0,
                 bounds=None, tol=1e-7,
                 options={'gtol': 1e-7, 'ftol': 1e-7,
                          'eps': 1e-7, 'maxiter': 1000}):
        """
        Initialize an IVIM model.

        The IVIM model assumes that biological tissue includes a volume
        fraction 'f' of water flowing in perfused capillaries, with a
        perfusion coefficient D* and a fraction (1-f) of static (diffusion
        only), intra and extracellular water, with a diffusion coefficient
        D. In this model the echo attenuation of a signal in a single voxel
        can be written as

            .. math::

            S(b) = S_0[f*e^{(-b*D\*)} + (1-f)e^{(-b*D)}]

            Where:
            .. math::

            S_0, f, D\* and D are the IVIM parameters.

        Parameters
        ----------
        gtab : GradientTable class instance
            Gradient directions and bvalues

        split_b : float, optional
            The b-value to split the data on for two-stage fit.
            default : 200.

        bounds : tuple of arrays with 4 elements, optional
            Bounds to constrain the fitted model parameters. This is only supported for
            Scipy version > 0.17. When using a older scipy version, this function will raise
            an error if bounds are different from None.
            default : ([0., 0., 0., 0.], [np.inf, 1., 1., 1.])

        tol : float, optional
            Tolerance for convergence of minimization.
            default : 1e-7

        options : dict, optional
            Dictionary containing gtol, ftol, eps and maxiter. This is passed
            to leastsq.
            default : options={'gtol': 1e-7, 'ftol': 1e-7, 'eps': 1e-7, 'maxiter': 1000}

        References
        ----------
        .. [1] <NAME>, Denis, et al. "Separation of diffusion
                   and perfusion in intravoxel incoherent motion MR
                   imaging." Radiology 168.2 (1988): 497-505.
        .. [2] <NAME>, et al. "Quantitative measurement
                   of brain perfusion with intravoxel incoherent motion
                   MR imaging." Radiology 265.3 (2012): 874-881.
        """
        ReconstModel.__init__(self, gtab)
        self.split_b = split_b
        self.bounds = bounds
        self.tol = tol
        self.options = options

        if SCIPY_LESS_0_17 and self.bounds is not None:
            e_s = "Scipy versions less than 0.17 do not support "
            e_s += "bounds. Please update to Scipy 0.17 to use bounds"
            raise ValueError(e_s)
        else:
            self.bounds = (np.array([0., 0., 0., 0.]),
                           np.array([np.inf, 1., 0.1, 0.1]))

    def estimate_x0(self, data):
        """
        Fit the ivim params using a two stage fit.

        In the two stage fitting routine, initially, we fit the signal
        at bvals higher than the specified `split_b` using `TensorModel`
        and get an intial guess for f and D. Then, using these parameters
        we fit the entire data for all bvalues.

        Parameters
        ----------
        data : array
            The measured signal from one voxel. A multi voxel decorator
            will be applied to this fit method to scale it and apply it
            to multiple voxels.

        Returns
        -------
        x0_guess : array
            An array with initial values of S0, f, D_star, D for each voxel.
        """
        S0_hat, D_guess = self._estimate_S0_D(data)
        f_guess = 1 - S0_hat / np.mean(data[self.gtab.b0s_mask])
        # We set the S0 guess as the signal value at b=0
        # The D* guess is roughly 10 times the D value
        x0 = np.array([np.mean(data[self.gtab.b0s_mask]), f_guess, 10 * D_guess, D_guess])
        # The API does not allow bounds for Scipy < 0.17. While estimating x0,
        # if there is noise in the data the estimated x0 might not be feasible.
        # In such a case we will use the values for the parameters
{'fraction_non_alphanumeric': 0.08906075586030936, 'fraction_numerical': 0.026311593047360866, 'mean_word_length': 1.9842969307637401, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '311312', 'n_tokens_mistral': 3944, 'n_tokens_neox': 3583, 'n_words': 1387}