diff --git a/scripts/dwi_estimation_error_analysis.py b/scripts/dwi_estimation_error_analysis.py index d5070d69..fab59963 100644 --- a/scripts/dwi_estimation_error_analysis.py +++ b/scripts/dwi_estimation_error_analysis.py @@ -73,7 +73,7 @@ def cross_validate( """ gpm = EddyMotionGPR( - kernel=SphericalKriging(a=1.8, lambda_s=1000), + kernel=SphericalKriging(beta_a=1.8, beta_l=1000), alpha=100, disp=True, # optimizer="Nelder-Mead", diff --git a/src/eddymotion/model/_dipy.py b/src/eddymotion/model/_dipy.py index e64b9e27..1e9185b5 100644 --- a/src/eddymotion/model/_dipy.py +++ b/src/eddymotion/model/_dipy.py @@ -90,8 +90,8 @@ class GaussianProcessModel(ReconstModel): def __init__( self, kernel_model: str = "spherical", - lambda_s: float = 2.0, - a: float = 0.1, + beta_l: float = 2.0, + beta_a: float = 0.1, sigma_sq: float = 1.0, *args, **kwargs, @@ -119,8 +119,8 @@ def __init__( KernelType = SphericalKriging if kernel_model == "spherical" else ExponentialKriging self.kernel = KernelType( - a=a, - lambda_s=lambda_s, + beta_a=beta_a, + beta_l=beta_l, ) def fit( diff --git a/src/eddymotion/model/_sklearn.py b/src/eddymotion/model/_sklearn.py index db8ed0d2..3baa8a00 100644 --- a/src/eddymotion/model/_sklearn.py +++ b/src/eddymotion/model/_sklearn.py @@ -103,7 +103,7 @@ BOUNDS_A: tuple[float, float] = (0.1, np.pi) """The limits for the parameter *a*.""" -BOUNDS_LAMBDA_S: tuple[float, float] = (1.0, 1000) +BOUNDS_LAMBDA: tuple[float, float] = (1.0, 1000) """The limits for the parameter lambda.""" THETA_EPSILON: float = 1e-5 """Minimum nonzero angle.""" @@ -246,38 +246,38 @@ class ExponentialKriging(Kernel): def __init__( self, - a: float = 0.01, - lambda_s: float = 2.0, + beta_a: float = 0.01, + beta_l: float = 2.0, a_bounds: tuple[float, float] = BOUNDS_A, - lambda_s_bounds: tuple[float, float] = BOUNDS_LAMBDA_S, + l_bounds: tuple[float, float] = BOUNDS_LAMBDA, ): r""" Initialize an exponential Kriging kernel. Parameters ---------- - a : :obj:`float` + beta_a : :obj:`float` Minimum angle in rads. - lambda_s : :obj:`float` - The :math:`\lambda_s` hyperparameter. + beta_l : :obj:`float` + The :math:`\lambda` hyperparameter. a_bounds : :obj:`tuple` Bounds for the a parameter. - lambda_s_bounds : :obj:`tuple` - Bounds for the :math:`\lambda_s` hyperparameter. + l_bounds : :obj:`tuple` + Bounds for the :math:`\lambda` hyperparameter. """ - self.a = a - self.lambda_s = lambda_s + self.beta_a = beta_a + self.beta_l = beta_l self.a_bounds = a_bounds - self.lambda_s_bounds = lambda_s_bounds + self.l_bounds = l_bounds @property def hyperparameter_a(self) -> Hyperparameter: - return Hyperparameter("a", "numeric", self.a_bounds) + return Hyperparameter("beta_a", "numeric", self.a_bounds) @property - def hyperparameter_lambda_s(self) -> Hyperparameter: - return Hyperparameter("lambda_s", "numeric", self.lambda_s_bounds) + def hyperparameter_beta_l(self) -> Hyperparameter: + return Hyperparameter("beta_l", "numeric", self.l_bounds) def __call__( self, X: np.ndarray, Y: np.ndarray | None = None, eval_gradient: bool = False @@ -310,16 +310,16 @@ def __call__( """ thetas = compute_pairwise_angles(X, Y) thetas[np.abs(thetas) < THETA_EPSILON] = 0.0 - C_theta = np.exp(-thetas / self.a) + C_theta = np.exp(-thetas / self.beta_a) if not eval_gradient: - return self.lambda_s * C_theta + return self.beta_l * C_theta K_gradient = np.zeros((*thetas.shape, 2)) - K_gradient[..., 0] = self.lambda_s * C_theta * thetas / self.a**2 # Derivative w.r.t. a + K_gradient[..., 0] = self.beta_l * C_theta * thetas / self.beta_a**2 # Derivative w.r.t. a K_gradient[..., 1] = C_theta - return self.lambda_s * C_theta, K_gradient + return self.beta_l * C_theta, K_gradient def diag(self, X: np.ndarray) -> np.ndarray: """Returns the diagonal of the kernel k(X, X). @@ -338,14 +338,14 @@ def diag(self, X: np.ndarray) -> np.ndarray: K_diag : ndarray of shape (n_samples_X,) Diagonal of kernel k(X, X) """ - return self.lambda_s * np.ones(X.shape[0]) + return self.beta_l * np.ones(X.shape[0]) def is_stationary(self) -> bool: """Returns whether the kernel is stationary.""" return True def __repr__(self) -> str: - return f"ExponentialKriging (a={self.a}, λₛ={self.lambda_s})" + return f"ExponentialKriging (a={self.beta_a}, λ={self.beta_l})" class SphericalKriging(Kernel): @@ -353,38 +353,38 @@ class SphericalKriging(Kernel): def __init__( self, - a: float = 0.01, - lambda_s: float = 2.0, + beta_a: float = 0.01, + beta_l: float = 2.0, a_bounds: tuple[float, float] = BOUNDS_A, - lambda_s_bounds: tuple[float, float] = BOUNDS_LAMBDA_S, + l_bounds: tuple[float, float] = BOUNDS_LAMBDA, ): r""" Initialize a spherical Kriging kernel. Parameters ---------- - a : :obj:`float` + beta_a : :obj:`float` Minimum angle in rads. - lambda_s : :obj:`float` - The :math:`\lambda_s` hyperparameter. + beta_l : :obj:`float` + The :math:`\lambda` hyperparameter. a_bounds : :obj:`tuple` Bounds for the a parameter. - lambda_s_bounds : :obj:`tuple` - Bounds for the :math:`\lambda_s` hyperparameter. + l_bounds : :obj:`tuple` + Bounds for the :math:`\lambda` hyperparameter. """ - self.a = a - self.lambda_s = lambda_s + self.beta_a = beta_a + self.beta_l = beta_l self.a_bounds = a_bounds - self.lambda_s_bounds = lambda_s_bounds + self.l_bounds = l_bounds @property def hyperparameter_a(self) -> Hyperparameter: - return Hyperparameter("a", "numeric", self.a_bounds) + return Hyperparameter("beta_a", "numeric", self.a_bounds) @property - def hyperparameter_lambda_s(self) -> Hyperparameter: - return Hyperparameter("lambda_s", "numeric", self.lambda_s_bounds) + def hyperparameter_beta_l(self) -> Hyperparameter: + return Hyperparameter("beta_l", "numeric", self.l_bounds) def __call__( self, X: np.ndarray, Y: np.ndarray | None = None, eval_gradient: bool = False @@ -418,23 +418,25 @@ def __call__( thetas = compute_pairwise_angles(X, Y) thetas[np.abs(thetas) < THETA_EPSILON] = 0.0 - nonzero = thetas <= self.a + nonzero = thetas <= self.beta_a C_theta = np.zeros_like(thetas) C_theta[nonzero] = ( - 1 - 1.5 * thetas[nonzero] / self.a + 0.5 * thetas[nonzero] ** 3 / self.a**3 + 1 - 1.5 * thetas[nonzero] / self.beta_a + 0.5 * thetas[nonzero] ** 3 / self.beta_a**3 ) if not eval_gradient: - return self.lambda_s * C_theta + return self.beta_l * C_theta deriv_a = np.zeros_like(thetas) deriv_a[nonzero] = ( - 1.5 * self.lambda_s * (thetas[nonzero] / self.a**2 - thetas[nonzero] ** 3 / self.a**4) + 1.5 + * self.beta_l + * (thetas[nonzero] / self.beta_a**2 - thetas[nonzero] ** 3 / self.beta_a**4) ) K_gradient = np.dstack((deriv_a, C_theta)) - return self.lambda_s * C_theta, K_gradient + return self.beta_l * C_theta, K_gradient def diag(self, X: np.ndarray) -> np.ndarray: """Returns the diagonal of the kernel k(X, X). @@ -453,14 +455,14 @@ def diag(self, X: np.ndarray) -> np.ndarray: K_diag : ndarray of shape (n_samples_X,) Diagonal of kernel k(X, X) """ - return self.lambda_s * np.ones(X.shape[0]) + return self.beta_l * np.ones(X.shape[0]) def is_stationary(self) -> bool: """Returns whether the kernel is stationary.""" return True def __repr__(self) -> str: - return f"SphericalKriging (a={self.a}, λₛ={self.lambda_s})" + return f"SphericalKriging (a={self.beta_a}, λ={self.beta_l})" def compute_pairwise_angles(