diff --git a/.github/workflows/analysis.yml b/.github/workflows/analysis.yml new file mode 100644 index 0000000..ba1d26f --- /dev/null +++ b/.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 diff --git a/README.md b/README.md index 5cb4185..6f9c976 100644 --- a/README.md +++ b/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) diff --git a/conftest.py b/conftest.py index ce16f10..fd6e0ea 100644 --- a/conftest.py +++ b/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: diff --git a/src/original/PV_MUMC/two_step_IVIM_fit.py b/src/original/PV_MUMC/two_step_IVIM_fit.py index db9d6d6..7007a60 100644 --- a/src/original/PV_MUMC/two_step_IVIM_fit.py +++ b/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 diff --git a/src/standardized/PV_MUMC_biexp.py b/src/standardized/PV_MUMC_biexp.py new file mode 100644 index 0000000..1425aee --- /dev/null +++ b/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 diff --git a/src/standardized/PvH_KB_NKI_IVIMfit.py b/src/standardized/PvH_KB_NKI_IVIMfit.py new file mode 100644 index 0000000..3d21ad5 --- /dev/null +++ b/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 diff --git a/tests/IVIMmodels/unit_tests/algorithms.json b/tests/IVIMmodels/unit_tests/algorithms.json index 06fe41b..e4320f5 100644 --- a/tests/IVIMmodels/unit_tests/algorithms.json +++ b/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": { diff --git a/tests/IVIMmodels/unit_tests/analyze.r b/tests/IVIMmodels/unit_tests/analyze.r index 51bdbf2..27025a6 100644 --- a/tests/IVIMmodels/unit_tests/analyze.r +++ b/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) { diff --git a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py index fe7f339..92473bb 100644 --- a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py +++ b/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) diff --git a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py index 37ca940..cf4b09a 100644 --- a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py +++ b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py @@ -14,7 +14,7 @@ #run using pytest --saveFileName test_output.txt --SNR 50 100 200 #e.g. pytest -m slow tests/IVIMmodels/unit_tests/test_ivim_synthetic.py --saveFileName test_output.csv --SNR 10 50 100 200 --fitCount 20 @pytest.mark.slow -def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, rician_noise, save_file, save_duration_file): +def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, rician_noise, save_file, save_duration_file, use_prior): # assert save_file == "test" random.seed(42) S0 = 1 @@ -24,8 +24,8 @@ def test_generated(ivim_algorithm, ivim_data, SNR, rtol, atol, fit_count, rician f = data["f"] Dp = data["Dp"] fit = OsipiBase(algorithm=ivim_algorithm) - # here try a prior, but it's not seeing it, why? - if hasattr(fit, "accepts_priors") and fit.accepts_priors: + # here is a prior + if use_prior and hasattr(fit, "accepts_priors") and fit.accepts_priors: prior = [np.random.normal(D, D/3, 10), np.random.normal(f, f/3, 10), np.random.normal(Dp, Dp/3, 10), np.random.normal(1, 1/3, 10)] # prior = [np.repeat(D, 10)+np.random.normal(0,D/3,np.shape(np.repeat(D, 10))), np.repeat(f, 10)+np.random.normal(0,f/3,np.shape(np.repeat(D, 10))), np.repeat(Dp, 10)+np.random.normal(0,Dp/3,np.shape(np.repeat(D, 10))),np.repeat(np.ones_like(Dp), 10)+np.random.normal(0,1/3,np.shape(np.repeat(D, 10)))] # D, f, D* fit.initialize(prior_in=prior)