Merge remote-tracking branch 'upstream/main' into AddGPExperimentScripts

docs/developers.rst

@@ -38,3 +38,4 @@ Information on specific functions, classes, and methods.
    api/eddymotion.registration
    api/eddymotion.testing
    api/eddymotion.utils
+   api/eddymotion.viz

pyproject.toml

@@ -178,6 +178,7 @@ concurrency = ['multiprocessing']
 omit = [
     '*/tests/*',
     '*/testing/*',
+    '*/viz/*',
     '*/__init__.py',
     '*/conftest.py',
     'src/eddymotion/_version.py'

scripts/dwi_estimation_error_analysis.py

@@ -30,14 +30,10 @@
 
 import argparse
 
-import matplotlib.gridspec as gridspec
-
 # import nibabel as nib
 import numpy as np
 from dipy.core.gradients import gradient_table
 from dipy.sims.voxel import single_tensor
-from matplotlib import pyplot as plt
-from scipy.stats import pearsonr
 from sklearn.metrics import root_mean_squared_error
 from sklearn.model_selection import KFold, RepeatedKFold, cross_val_score
 
@@ -280,92 +276,6 @@ def compute_error(
     return np.asarray(mean_rmse), np.asarray(std_dev)
 
 
-def plot_error(
-    kfolds: list[int], mean: np.ndarray, std_dev: np.ndarray, xlabel: str, ylabel: str, title: str
-) -> plt.Figure:
-    """
-    Plot the error and standard deviation.
-
-    Parameters
-    ----------
-    kfolds : :obj:`list`
-        Number of k-folds.
-    mean : :obj:`~numpy.ndarray`
-        Mean RMSE values.
-    std_dev : :obj:`~numpy.ndarray`
-        Standard deviation values.
-    xlabel : :obj:`str`
-        X-axis label.
-    ylabel : :obj:`str`
-        Y-axis label.
-    title : :obj:`str`
-        Plot title.
-
-    Returns
-    -------
-    :obj:`~matplotlib.pyplot.Figure`
-        Matplotlib figure object.
-
-    """
-    fig, ax = plt.subplots()
-    ax.plot(kfolds, mean, c="orange")
-    ax.fill_between(kfolds, mean - std_dev, mean + std_dev, alpha=0.5, color="orange")
-    ax.scatter(kfolds, mean, c="orange")
-    ax.set_xlabel(xlabel)
-    ax.set_ylabel(ylabel)
-    ax.set_xticks(kfolds)
-    ax.set_xticklabels(kfolds)
-    ax.set_title(title)
-    fig.tight_layout()
-    return fig
-
-
-def plot_estimation_carpet(gt_nii, gp_nii, gtab, suptitle, **kwargs):
-    from nireports.reportlets.modality.dwi import nii_to_carpetplot_data
-    from nireports.reportlets.nuisance import plot_carpet
-
-    fig = plt.figure(layout="tight")
-    gs = gridspec.GridSpec(ncols=1, nrows=2, figure=fig)
-    fig.suptitle(suptitle)
-
-    divide_by_b0 = False
-    gt_data, segments = nii_to_carpetplot_data(gt_nii, bvals=gtab.bvals, divide_by_b0=divide_by_b0)
-
-    title = "Ground truth"
-    plot_carpet(gt_data, segments, subplot=gs[0, :], title=title, **kwargs)
-
-    gp_data, segments = nii_to_carpetplot_data(gp_nii, bvals=gtab.bvals, divide_by_b0=divide_by_b0)
-
-    title = "Estimated (GP)"
-    plot_carpet(gt_data, segments, subplot=gs[1, :], title=title, **kwargs)
-
-    return fig
-
-
-def plot_correlation(x, y, title):
-    r = pearsonr(x, y)
-
-    # Fit a linear curve and estimate its y-values and their error
-    a, b = np.polyfit(x, y, deg=1)
-    y_est = a * x + b
-    y_err = x.std() * np.sqrt(1 / len(x) + (x - x.mean()) ** 2 / np.sum((x - x.mean()) ** 2))
-
-    fig, ax = plt.subplots()
-    ax.plot(x, y_est, "-", color="black", label=f"r = {r.correlation:.2f}")
-    ax.fill_between(x, y_est - y_err, y_est + y_err, alpha=0.2, color="lightgray")
-    ax.plot(x, y, marker="o", markersize="4", color="gray")
-
-    ax.set_ylabel("Ground truth")
-    ax.set_xlabel("Estimated")
-
-    plt.title(title)
-    plt.legend(loc="lower right")
-
-    fig.tight_layout()
-
-    return fig, r
-
-
 def _build_arg_parser() -> argparse.ArgumentParser:
     """
     Build argument parser for command-line interface.
@@ -440,7 +350,7 @@ def main() -> None:
     # Compute the error
     # rmse, std_dev = compute_error(data, args.repeats, args.kfold)
 
-    # Plot
+    # Plot - import plot_* from eddymotion.viz.signals
     # xlabel = "N"
     # ylabel = "RMSE"
     # title = f"Gaussian process estimation\n(SNR={args.snr})"

src/eddymotion/model/_dipy.py

@@ -67,7 +67,7 @@ def gp_prediction(
 
     """
 
-    X = gtab.bvecs.T if hasattr("bvecs") else np.asarray(gtab)
+    X = gtab.bvecs.T if hasattr(gtab, "bvecs") else np.asarray(gtab)
 
     # Check it's fitted as they do in sklearn internally
     # https://github.com/scikit-learn/scikit-learn/blob/972e17fe1aa12d481b120ad4a3dc076bae736931/\
@@ -151,7 +151,7 @@ def fit(
 
         # Extract b-vecs: Scikit learn wants (n_samples, n_features)
         # where n_features is 3, and n_samples the different diffusion orientations.
-        X = gtab.bvecs.T if hasattr("bvecs") else np.asarray(gtab)
+        X = gtab.bvecs if hasattr(gtab, "bvecs") else np.asarray(gtab)
 
         # Data must be shapes (n_samples, n_targets) where n_samples is
         # the number of diffusion orientations, and n_targets is number of voxels.

src/eddymotion/viz/signals.py

@@ -0,0 +1,114 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# © The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Visualizing signals and intermediate aspects of models."""
+
+import matplotlib.gridspec as gridspec
+import numpy as np
+from matplotlib import pyplot as plt
+from scipy.stats import pearsonr
+
+
+def plot_error(
+    kfolds: list[int], mean: np.ndarray, std_dev: np.ndarray, xlabel: str, ylabel: str, title: str
+) -> plt.Figure:
+    """
+    Plot the error and standard deviation.
+
+    Parameters
+    ----------
+    kfolds : :obj:`list`
+        Number of k-folds.
+    mean : :obj:`~numpy.ndarray`
+        Mean RMSE values.
+    std_dev : :obj:`~numpy.ndarray`
+        Standard deviation values.
+    xlabel : :obj:`str`
+        X-axis label.
+    ylabel : :obj:`str`
+        Y-axis label.
+    title : :obj:`str`
+        Plot title.
+
+    Returns
+    -------
+    :obj:`~matplotlib.pyplot.Figure`
+        Matplotlib figure object.
+
+    """
+    fig, ax = plt.subplots()
+    ax.plot(kfolds, mean, c="orange")
+    ax.fill_between(kfolds, mean - std_dev, mean + std_dev, alpha=0.5, color="orange")
+    ax.scatter(kfolds, mean, c="orange")
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel(ylabel)
+    ax.set_xticks(kfolds)
+    ax.set_xticklabels(kfolds)
+    ax.set_title(title)
+    fig.tight_layout()
+    return fig
+
+
+def plot_estimation_carpet(gt_nii, gp_nii, gtab, suptitle, **kwargs):
+    from nireports.reportlets.modality.dwi import nii_to_carpetplot_data
+    from nireports.reportlets.nuisance import plot_carpet
+
+    fig = plt.figure(layout="tight")
+    gs = gridspec.GridSpec(ncols=1, nrows=2, figure=fig)
+    fig.suptitle(suptitle)
+
+    divide_by_b0 = False
+    gt_data, segments = nii_to_carpetplot_data(gt_nii, bvals=gtab.bvals, divide_by_b0=divide_by_b0)
+
+    title = "Ground truth"
+    plot_carpet(gt_data, segments, subplot=gs[0, :], title=title, **kwargs)
+
+    gp_data, segments = nii_to_carpetplot_data(gp_nii, bvals=gtab.bvals, divide_by_b0=divide_by_b0)
+
+    title = "Estimated (GP)"
+    plot_carpet(gt_data, segments, subplot=gs[1, :], title=title, **kwargs)
+
+    return fig
+
+
+def plot_correlation(x, y, title):
+    r = pearsonr(x, y)
+
+    # Fit a linear curve and estimate its y-values and their error
+    a, b = np.polyfit(x, y, deg=1)
+    y_est = a * x + b
+    y_err = x.std() * np.sqrt(1 / len(x) + (x - x.mean()) ** 2 / np.sum((x - x.mean()) ** 2))
+
+    fig, ax = plt.subplots()
+    ax.plot(x, y_est, "-", color="black", label=f"r = {r.correlation:.2f}")
+    ax.fill_between(x, y_est - y_err, y_est + y_err, alpha=0.2, color="lightgray")
+    ax.plot(x, y, marker="o", markersize="4", color="gray")
+
+    ax.set_ylabel("Ground truth")
+    ax.set_xlabel("Estimated")
+
+    plt.title(title)
+    plt.legend(loc="lower right")
+
+    fig.tight_layout()
+
+    return fig, r

test/test_model.py

@@ -24,15 +24,15 @@
 
 import numpy as np
 import pytest
-from dipy.core.gradients import gradient_table
-from sklearn.datasets import make_regression
+from dipy.sims.voxel import single_tensor
 
 from eddymotion import model
 from eddymotion.data.dmri import DWI
 from eddymotion.data.splitting import lovo_split
 from eddymotion.exceptions import ModelNotFittedError
 from eddymotion.model._dipy import GaussianProcessModel
 from eddymotion.model.dmri import DEFAULT_MAX_S0, DEFAULT_MIN_S0
+from eddymotion.testing import simulations as _sim
 
 
 def test_trivial_model():
@@ -108,22 +108,39 @@ def test_average_model():
     assert np.all(tmodel_2000.predict([0, 0, 0]) == 1100)
 
 
-def test_gp_model():
-    gp = GaussianProcessModel("test")
-
+@pytest.mark.parametrize(
+    (
+        "bval_shell",
+        "S0",
+        "evals",
+    ),
+    [
+        (
+            1000,
+            100,
+            (0.0015, 0.0003, 0.0003),
+        )
+    ],
+)
+@pytest.mark.parametrize("snr", (10, 20))
+@pytest.mark.parametrize("hsph_dirs", (60, 30))
+def test_gp_model(evals, S0, snr, hsph_dirs, bval_shell):
+    # Simulate signal for a single tensor
+    evecs = _sim.create_single_fiber_evecs()
+    gtab = _sim.create_single_shell_gradient_table(hsph_dirs, bval_shell)
+    signal = single_tensor(gtab, S0=S0, evals=evals, evecs=evecs, snr=snr)
+
+    # Drop the initial b=0
+    gtab = gtab[1:]
+    data = signal[1:]
+
+    gp = GaussianProcessModel(kernel_model="spherical")
     assert isinstance(gp, model._dipy.GaussianProcessModel)
 
-    X, y = make_regression(n_samples=100, n_features=3, noise=0, random_state=0)
-
-    bvecs = X.T / np.linalg.norm(X.T, axis=0)
-    gtab = gradient_table([1000] * bvecs.shape[-1], bvecs)
-
-    gpfit = gp.fit(gtab, y)
-
-    X_qry = bvecs[:, :2].T
-    prediction = gpfit.predict(X_qry)
+    gpfit = gp.fit(data[:-2], gtab[:-2])
+    prediction = gpfit.predict(gtab.bvecs[-2:])
 
-    assert prediction.shape == (X_qry.shape[0],)
+    assert prediction.shape == (2,)
 
 
 def test_two_initialisations(datadir):

test/test_sklearn.py

@@ -24,7 +24,6 @@
 
 import numpy as np
 import pytest
-from dipy.core.gradients import gradient_table
 from dipy.io import read_bvals_bvecs
 
 from eddymotion.model import _sklearn as ems
@@ -287,7 +286,7 @@ def test_kernel(repodata, covariance):
     kernel = KernelType()
     K = kernel(bvecs)
 
-    assert K.shape == (bvecs.shape[0], ) * 2
+    assert K.shape == (bvecs.shape[0],) * 2
 
     assert np.allclose(np.diagonal(K), kernel.diag(bvecs))