Analysis/workflow

.github/workflows/analysis.yml

@@ -0,0 +1,63 @@
+name: Algorithm Analysis
+
+on:
+  push:
+    branches:
+      - 'main'
+      - 'analysis/**'
+
+jobs:
+  build:
+
+    runs-on: ubuntu-latest
+    continue-on-error: false
+    strategy:
+      fail-fast: false
+    steps:
+      - uses: actions/checkout@v3
+      - name: Set up Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: "3.11"
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -r requirements.txt
+      - name: Set up R
+        uses: r-lib/actions/setup-r@v2
+        with: 
+          use-public-rspm: true
+      - name: Install R dependencies
+        uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          packages: |
+            any::plyr
+            any::dplyr
+            any::tidyverse
+            any::data.table
+            any::ggplot2
+      - name: Generate fitting data
+        run: |
+          pip install pytest
+          python -m pytest -m slow --saveFileName test_output.csv --SNR 10 30 50 100 200 --fitCount 300 --saveDurationFileName test_duration.csv
+      - name: Generate figures
+        run: Rscript --vanilla tests/IVIMmodels/unit_tests/analyze.r test_output.csv test_duration.csv
+      - name: Upload raw data
+        uses: actions/upload-artifact@v3
+        with:
+            name: Raw data
+            path: |
+              test_output.csv
+              test_duration.csv
+      - name: Upload figures
+        uses: actions/upload-artifact@v3
+        with:
+            name: Fit figures
+            path: |
+              D.pdf
+              f.pdf
+              Dp.pdf
+              D_limited.pdf
+              f_limited.pdf
+              Dp_limited.pdf
+              durations.pdf

README.md

@@ -30,4 +30,4 @@ The **test** folder contains the test files corresponding to the contributed cod
 The **utils** folder contains various helpful tools.
 
 ## View Testing Reports
-*to be added*
+[![Unit tests](https://github.com/OSIPI/TF2.4_IVIM-MRI_CodeCollection/actions/workflows/unit_test.yml/badge.svg?branch=main)](https://github.com/OSIPI/TF2.4_IVIM-MRI_CodeCollection/actions/workflows/unit_test.yml)

conftest.py

@@ -19,6 +19,12 @@ def pytest_addoption(parser):
         type=bool,
         help="Use Rician noise, non-rician is gaussian",
     )
+    parser.addoption(
+        "--usePrior",
+        default=False,
+        type=bool,
+        help="Use a prior where accepted",
+    )
     parser.addoption(
         "--algorithmFile",
         default="tests/IVIMmodels/unit_tests/algorithms.json",
@@ -110,6 +116,10 @@ def fit_count(request):
 def rician_noise(request):
     return request.config.getoption("--ricianNoise")
 
+@pytest.fixture(scope="session")
+def use_prior(request):
+    return request.config.getoption("--usePrior")
+
 
 def pytest_generate_tests(metafunc):
     if "SNR" in metafunc.fixturenames:

src/original/PV_MUMC/two_step_IVIM_fit.py

@@ -7,36 +7,34 @@
 numpy
 tqdm
 scipy
-joblib
 """
 
 # load relevant libraries
-from scipy.optimize import curve_fit, nnls
+from scipy.optimize import curve_fit
 import numpy as np
-from joblib import Parallel, delayed
 import tqdm
 
 
 
 
 def two_exp_noS0(bvalues, Dpar, Fmv, Dmv):
-    """ tri-exponential IVIM function, and S0 set to 1"""
+    """ bi-exponential IVIM function, and S0 set to 1"""
     return Fmv * np.exp(-bvalues * Dmv) + (1 - Fmv ) * np.exp(-bvalues * Dpar)
 
 def two_exp(bvalues, S0, Dpar, Fmv, Dmv):
-    """ tri-exponential IVIM function"""
+    """ bi-exponential IVIM function"""
     return S0 * (Fmv * np.exp(-bvalues * Dmv) + (1 - Fmv ) * np.exp(-bvalues * Dpar))
 
 
 
-def fit_least_squares_array(bvalues, dw_data, fitS0=True, bounds=([0.9, 0.0001, 0.0, 0.0025], [1.1, 0.0025, 0.2, 0.2]), cutoff=200):
+def fit_least_squares_array(bvalues, dw_data, fitS0=False, bounds=([0.9, 0.0001, 0.0, 0.0025], [1.1, 0.0025, 0.5, 0.2]), cutoff=200):
     """
     This is the LSQ implementation, in which we first estimate Dpar using a curve fit to b-values>=cutoff;
     Second, we fit the other parameters using all b-values, while fixing Dpar from step 1. This fit
     is done on an array.
     :param bvalues: 1D Array with the b-values
     :param dw_data: 2D Array with diffusion-weighted signal in different voxels at different b-values
-    :param bounds: Array with fit bounds ([S0min, Dparmin, Fintmin, Dintmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fintmax, Dintmax, Fmvmax, Dmvmax]). default: ([0.9, 0.0001, 0.0, 0.0015, 0.0, 0.004], [1.1, 0.0015, 0.4, 0.004, 0.2, 0.2])
+    :param bounds: Array with fit bounds ([S0min, Dparmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fmvmax, Dmvmax]). 
     :param cutoff: cutoff b-value used in step 1 
     :return Dpar: 1D Array with Dpar in each voxel
     :return Fmv: 1D Array with Fmv in each voxel
@@ -54,69 +52,56 @@ def fit_least_squares_array(bvalues, dw_data, fitS0=True, bounds=([0.9, 0.0001,
     return [Dpar, Fmv, Dmv, S0]
 
 
-def fit_least_squares(bvalues, dw_data, IR=False, S0_output=False, fitS0=True,
-                      bounds=([0.9, 0.0001, 0.0, 0.0025], [1.1, 0.0025, 0.2, 0.2]), cutoff=200):
+def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=False,
+                      bounds=([0.9, 0.0001, 0.0, 0.0025], [1.1, 0.003, 1, 0.2]), cutoff=200):
     """
    This is the LSQ implementation, in which we first estimate Dpar using a curve fit to b-values>=cutoff;
    Second, we fit the other parameters using all b-values, while fixing Dpar from step 1. This fit
    is done on an array. It fits a single curve
     :param bvalues: 1D Array with the b-values
     :param dw_data: 1D Array with diffusion-weighted signal in different voxels at different b-values
-    :param IR: Boolean; True will fit the IVIM accounting for inversion recovery, False will fit IVIM without IR; default = True
-    :param S0_output: Boolean determining whether to output (often a dummy) variable S0; default = False
-    :param fix_S0: Boolean determining whether to fix S0 to 1; default = True
-    :param bounds: Array with fit bounds ([S0min, Dparmin, Fintmin, Dintmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fintmax, Dintmax, Fmvmax, Dmvmax]). Default: ([0, 0, 0, 0.005, 0, 0.06], [2.5, 0.005, 1, 0.06, 1, 0.5])
+    :param S0_output: Boolean determining whether to output (often a dummy) variable S0; 
+    :param fitS0: Boolean determining whether to fix S0 to 1;
+    :param bounds: Array with fit bounds ([S0min, Dparmin, Fmvmin, Dmvmin],[S0max, Dparmax, Fmvmax, Dmvmax]). 
     :param cutoff: cutoff b-value used in step 1 
     :return S0: optional 1D Array with S0 in each voxel
     :return Dpar: scalar with Dpar of the specific voxel
-    :return Fint: scalar with Fint of the specific voxel
-    :return Dint: scalar with Dint of the specific voxel
     :return Fmv: scalar with Fmv of the specific voxel
     :return Dmv: scalar with Dmv of the specific voxel
     """
 
     #try:
-    def monofit(bvalues, Dpar):
-            return np.exp(-bvalues * Dpar)
+    def monofit(bvalues, Dpar, Fmv):
+            return (1-Fmv)*np.exp(-bvalues * Dpar)
 
     high_b = bvalues[bvalues >= cutoff]
     high_dw_data = dw_data[bvalues >= cutoff]
-    boundspar = ([bounds[0][1]], [bounds[1][1]])
-    params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2], bounds=boundspar)
-    Dpar = params[0]
-
+    boundsmonoexp = ([bounds[0][1], bounds[0][2]], [bounds[1][1], bounds[1][2]])
+    params, _ = curve_fit(monofit, high_b, high_dw_data, p0=[(bounds[1][1]-bounds[0][1])/2,(bounds[1][2]-bounds[0][2])/2], bounds=boundsmonoexp)
+    Dpar1 = params[0]
     if not fitS0:
-        boundsupdated=([Dpar1 , bounds[0][2] , bounds[0][3] ],
-                    [Dpar1 , bounds[1][2] , bounds[1][3] ])    
-        params, _ = curve_fit(two_exp_noS0, bvalues, dw_data, p0=[Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
-        Dpar1, Fmv, Dmv = params[0], params[1], params[2]
+        boundsupdated=([bounds[0][2] , bounds[0][3] ],
+                    [bounds[1][2] , bounds[1][3] ])    
+        params, _ = curve_fit(lambda b, Fmv, Dmv: two_exp_noS0(b, Dpar1, Fmv, Dmv), bvalues, dw_data, p0=[(bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
+        Fmv, Dmv = params[0], params[1]
         #when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN)
-        if Fmv < 1e-4:
-            Dmv = float("NaN")
+        #if Fmv < 1e-4:
+        #    Dmv = float("NaN")
 
     else: 
-        #boundsupdated = ([bounds[0][0] , Dpar1 , bounds[0][2] , bounds[0][3] ],
-        #            [bounds[1][0] , Dpar1, bounds[1][2] , bounds[1][3] ])   
-        #params, _ = curve_fit(two_exp, bvalues, dw_data, p0=[1, Dpar1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
         boundsupdated = ([bounds[0][0] , bounds[0][2] , bounds[0][3] ],
                     [bounds[1][0] , bounds[1][2] , bounds[1][3] ])   
-        params, _ = curve_fit(lambda bvalues, S0, Fmv, Dmv: two_exp(bvalues, S0, Dpar, Fmv, Dmv), bvalues, dw_data, p0=[1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
+        params, _ = curve_fit(lambda b, S0, Fmv, Dmv: two_exp(b, S0, Dpar1, Fmv, Dmv), bvalues, dw_data, p0=[1, (bounds[0][2]+bounds[1][2])/2, (bounds[0][3]+bounds[1][3])/2], bounds=boundsupdated)
         S0 = params[0]
         Fmv, Dmv = params[1] , params[2]
         #when the fraction of a compartment equals zero (or very very small), the corresponding diffusivity is non-existing (=NaN)
-        if Fmv < 1e-4:
-            Dmv = float("NaN")     
+        #if Fmv < 1e-4:
+        #    Dmv = float("NaN")     
 
     if S0_output:
-        return Dpar, Fmv, Dmv, S0
-    else:
-        return Dpar, Fmv, Dmv
-    #except:
-
-    if S0_output:
-        return 0, 0, 0, 0, 0, 0
+        return Dpar1, Fmv, Dmv, S0
     else:
-        return 0, 0, 0, 0, 0
+        return Dpar1, Fmv, Dmv
 
 
 

src/standardized/PV_MUMC_biexp.py

@@ -0,0 +1,56 @@
+import numpy as np
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.PV_MUMC.two_step_IVIM_fit import fit_least_squares_array, fit_least_squares
+
+
+class PV_MUMC_biexp(OsipiBase):
+    """
+    Bi-exponential fitting algorithm by Paulien Voorter, Maastricht University
+    """
+
+    # Some basic stuff that identifies the algorithm
+    id_author = "Paulien Voorter MUMC"
+    id_algorithm_type = "Bi-exponential fit"
+    id_return_parameters = "f, D*, D"
+    id_units = "seconds per milli metre squared or milliseconds per micro metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0,0] # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds
+    required_bounds = False
+    required_bounds_optional = True # Bounds may not be required but are optional
+    required_initial_guess = False
+    required_initial_guess_optional = True
+    accepted_dimensions = 1 # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+
+    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, weighting=None, stats=False):
+        """
+            Everything this algorithm requires should be implemented here.
+            Number of segmentation thresholds, bounds, etc.
+
+            Our OsipiBase object could contain functions that compare the inputs with
+            the requirements.
+        """
+        super(PV_MUMC_biexp, self).__init__(bvalues, None, bounds, None)
+        self.PV_algorithm = fit_least_squares
+
+
+    def ivim_fit(self, signals, bvalues=None):
+        """Perform the IVIM fit
+
+        Args:
+            signals (array-like)
+            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
+
+        Returns:
+            _type_: _description_
+        """
+
+
+        fit_results = self.PV_algorithm(bvalues, signals)
+
+        f = fit_results[1]
+        Dstar = fit_results[2]
+        D = fit_results[0]
+
+        return f, Dstar, D

src/standardized/PvH_KB_NKI_IVIMfit.py

@@ -0,0 +1,62 @@
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.PvH_KB_NKI.DWI_functions_standalone import generate_IVIMmaps_standalone, generate_ADC_standalone
+import numpy as np
+
+class PvH_KB_NKI_IVIMfit(OsipiBase):
+    """
+    Bi-exponential fitting algorithm by Petra van Houdt and Koen Baas, NKI
+    """
+
+    # I'm thinking that we define default attributes for each submission like this
+    # And in __init__, we can call the OsipiBase control functions to check whether
+    # the user inputs fulfil the requirements
+
+    # Some basic stuff that identifies the algorithm
+    id_author = "Group Uulke van der Heide, NKI"
+    id_algorithm_type = "Bi-exponential fit"
+    id_return_parameters = "f, D*, D"
+    id_units = "seconds per milli metre squared or milliseconds per micro metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 4
+    required_thresholds = [0,
+                           0]  # Interval from "at least" to "at most", in case submissions allow a custom number of thresholds
+    required_bounds = False
+    required_bounds_optional = False  # Bounds may not be required but are optional
+    required_initial_guess = False
+    required_initial_guess_optional =False
+    accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
+
+    def __init__(self, bvalues=None, thresholds=None,bounds=None,initial_guess=None):
+        """
+            Everything this algorithm requires should be implemented here.
+            Number of segmentation thresholds, bounds, etc.
+
+            Our OsipiBase object could contain functions that compare the inputs with
+            the requirements.
+        """
+        super(PvH_KB_NKI_IVIMfit, self).__init__(bvalues, thresholds,bounds,initial_guess)
+        self.NKI_algorithm = generate_IVIMmaps_standalone
+
+
+    def ivim_fit(self, signals, bvalues=None):
+        """Perform the IVIM fit
+
+        Args:
+            signals (array-like)
+            bvalues (array-like, optional): b-values for the signals. If None, self.bvalues will be used. Default is None.
+
+        Returns:
+            _type_: _description_
+        """
+        #bvalues = np.array(bvalues)
+        bvalues = bvalues.tolist() #NKI code expects a list instead of nparray
+        # reshape signal as the NKI code expects a 4D array
+        signals = np.reshape(signals, (1, 1, 1, len(signals)))  # assuming that in this test the signals are always single voxel
+        fit_results = self.NKI_algorithm(signals,bvalues)
+
+        D = fit_results[0][0,0,0]/1000
+        f = fit_results[1][0,0,0]
+        Dstar = fit_results[2][0,0,0]/1000
+
+        return f, Dstar, D

tests/IVIMmodels/unit_tests/algorithms.json

@@ -9,7 +9,9 @@
         "IAR_LU_subtracted",
         "OGC_AmsterdamUMC_Bayesian_biexp",
         "OGC_AmsterdamUMC_biexp_segmented",
-        "OGC_AmsterdamUMC_biexp"
+        "OGC_AmsterdamUMC_biexp",
+        "PV_MUMC_biexp",
+        "PvH_KB_NKI_IVIMfit"
     ],
     "IAR_LU_biexp": {
         "tolerances": {

tests/IVIMmodels/unit_tests/analyze.r

@@ -25,7 +25,7 @@ plot_ivim <- function(data, fileExtension) {
     f_plot <- ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=f_fitted)) + geom_boxplot(color="red", aes(y=f)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0, 1) + ggtitle("Perfusion fraction grid") + ylab("Perfusion fraction")
     print(f_plot)
     ggsave(paste("f", fileExtension, sep=""), plot=f_plot, width = 50, height = 50, units = "cm")
-    D_plot <- ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=D_fitted)) + geom_boxplot(color="red", aes(y=D)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Diffusion grid") + ylab("Diffusion")
+    D_plot <- ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=D_fitted)) + geom_boxplot(color="red", aes(y=D)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0, 0.005) + ggtitle("Diffusion grid") + ylab("Diffusion")
     print(D_plot)
     ggsave(paste("D", fileExtension, sep=""), plot=D_plot, width = 50, height = 50, units = "cm")
     Dp_plot <- ggplot(data, aes(x=Algorithm)) + geom_boxplot(aes(y=Dp_fitted)) + geom_boxplot(color="red", aes(y=Dp)) + facet_grid(SNR ~ Region) + scale_x_discrete(guide = guide_axis(angle = 90)) + ylim(0, 0.25) + ggtitle("Perfusion grid") + ylab("Perfusion")
@@ -43,8 +43,8 @@ plot_ivim(data_restricted, "_limited.pdf")
 data_duration <- read.csv(duration_name)
 data_duration <- data_duration %>% mutate_if(is.character, as.factor)
 data_duration$ms <- data_duration$Duration..us./data_duration$Count/1000
-ggplot(data_duration, aes(x=Algorithm, y=ms)) + geom_boxplot() + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Fit Duration") + ylab("Time (ms)")
-ggsave("durations.pdf", width = 20, height = 20, units = "cm")
+duration_plot <- ggplot(data_duration, aes(x=Algorithm, y=ms)) + geom_boxplot() + scale_x_discrete(guide = guide_axis(angle = 90)) + ggtitle("Fit Duration") + ylab("Time (ms)")
+ggsave("durations.pdf", plot=duration_plot, width = 20, height = 20, units = "cm")
 
 
 if (runPrediction) {

tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py

@@ -3,7 +3,8 @@
 from pathlib import Path
 #from src.standardized.ETP_SRI_LinearFitting import ETP_SRI_LinearFitting
 from src.standardized.IAR_LU_biexp import IAR_LU_biexp
-from src.standardized.IAR_LU_segmented_2step import IAR_LU_segmented_2step
+#from src.standardized.IAR_LU_segmented_2step import IAR_LU_segmented_2step
+from src.standardized.PvH_KB_NKI_IVIMfit import PvH_KB_NKI_IVIMfit
 #from src.standardized.PV_MUMC_biexp import PV_MUMC_biexp
 
 ## Simple test code... 
@@ -26,7 +27,7 @@ def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
 
 #model1 = ETP_SRI_LinearFitting(thresholds=[200])
 #model2 = IAR_LU_biexp()
-model2 = IAR_LU_segmented_2step()
+model2 = PvH_KB_NKI_IVIMfit()
 
 #dev_test_run(model1, linear_fit_option=True)
-dev_test_run(model2, thresholds="Arbitrary kwarg for ivim_fit")
+dev_test_run(model2)