diff --git a/.gitignore b/.gitignore index faf3523..977ba28 100644 --- a/.gitignore +++ b/.gitignore @@ -19,8 +19,10 @@ md5sums.txt *.asv src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/test_data src/original/ASD_MemorialSloanKettering/MRI-QAMPER_IVIM/output_files - - +tests/IVIMmodels/unit_tests/*.log +junit/* +ivim_simulation.bval +ivim_simulation.bvec # Unit test / coverage reports .tox/ diff --git a/conftest.py b/conftest.py index 086585c..3c33bda 100644 --- a/conftest.py +++ b/conftest.py @@ -81,6 +81,24 @@ def pytest_addoption(parser): type=str, help="Drop this algorithm from the list" ) + parser.addoption( + "--withmatlab", + action="store_true", + default=False, + help="Run MATLAB-dependent tests" + ) + + +@pytest.fixture(scope="session") +def eng(request): + """Start and return a MATLAB engine session if --withmatlab is set.""" + if not request.config.getoption("--withmatlab"): + return None + import matlab.engine + print("Starting MATLAB engine...") + eng = matlab.engine.start_matlab() + print("MATLAB engine started.") + return eng @pytest.fixture(scope="session") @@ -149,25 +167,20 @@ def use_prior(request): def pytest_generate_tests(metafunc): if "SNR" in metafunc.fixturenames: metafunc.parametrize("SNR", metafunc.config.getoption("SNR")) - if "ivim_algorithm" in metafunc.fixturenames: - algorithms = algorithm_list(metafunc.config.getoption("algorithmFile"), metafunc.config.getoption("selectAlgorithm"), metafunc.config.getoption("dropAlgorithm")) - metafunc.parametrize("ivim_algorithm", algorithms) if "ivim_data" in metafunc.fixturenames: data = data_list(metafunc.config.getoption("dataFile")) metafunc.parametrize("ivim_data", data) + if "data_ivim_fit_saved" in metafunc.fixturenames: + args = data_ivim_fit_saved(metafunc.config.getoption("dataFile"),metafunc.config.getoption("algorithmFile")) + metafunc.parametrize("data_ivim_fit_saved", args) + if "algorithmlist" in metafunc.fixturenames: + args = algorithmlist(metafunc.config.getoption("algorithmFile")) + metafunc.parametrize("algorithmlist", args) + if "bound_input" in metafunc.fixturenames: + args = bound_input(metafunc.config.getoption("dataFile"),metafunc.config.getoption("algorithmFile")) + metafunc.parametrize("bound_input", args) -def algorithm_list(filename, selected, dropped): - current_folder = pathlib.Path.cwd() - algorithm_path = current_folder / filename - with algorithm_path.open() as f: - algorithm_information = json.load(f) - algorithms = set(algorithm_information["algorithms"]) - algorithms = algorithms - set(dropped) - if len(selected) > 0 and selected[0]: - algorithms = algorithms & set(selected) - return list(algorithms) - def data_list(filename): current_folder = pathlib.Path.cwd() data_path = current_folder / filename @@ -178,3 +191,69 @@ def data_list(filename): bvals = bvals['bvalues'] for name, data in all_data.items(): yield name, bvals, data + + +def data_ivim_fit_saved(datafile, algorithmFile): + # Find the algorithms from algorithms.json + current_folder = pathlib.Path.cwd() + algorithm_path = current_folder / algorithmFile + with algorithm_path.open() as f: + algorithm_information = json.load(f) + # Load generic test data generated from the included phantom: phantoms/MR_XCAT_qMRI + generic = current_folder / datafile + with generic.open() as f: + all_data = json.load(f) + algorithms = algorithm_information["algorithms"] + bvals = all_data.pop('config') + bvals = bvals['bvalues'] + for algorithm in algorithms: + first = True + for name, data in all_data.items(): + algorithm_dict = algorithm_information.get(algorithm, {}) + xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}), + "strict": algorithm_dict.get("xfail_names", {}).get(name, True)} + kwargs = algorithm_dict.get("options", {}) + tolerances = algorithm_dict.get("tolerances", {}) + skiptime=False + if first: + if algorithm_dict.get("fail_first_time", False): + skiptime = True + first = False + requires_matlab = algorithm_dict.get("requires_matlab", False) + yield name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime, requires_matlab + +def algorithmlist(algorithmFile): + # Find the algorithms from algorithms.json + current_folder = pathlib.Path.cwd() + algorithm_path = current_folder / algorithmFile + with algorithm_path.open() as f: + algorithm_information = json.load(f) + + algorithms = algorithm_information["algorithms"] + for algorithm in algorithms: + algorithm_dict = algorithm_information.get(algorithm, {}) + requires_matlab = algorithm_dict.get("requires_matlab", False) + yield algorithm, requires_matlab + +def bound_input(datafile,algorithmFile): + # Find the algorithms from algorithms.json + current_folder = pathlib.Path.cwd() + algorithm_path = current_folder / algorithmFile + with algorithm_path.open() as f: + algorithm_information = json.load(f) + # Load generic test data generated from the included phantom: phantoms/MR_XCAT_qMRI + generic = current_folder / datafile + with generic.open() as f: + all_data = json.load(f) + algorithms = algorithm_information["algorithms"] + bvals = all_data.pop('config') + bvals = bvals['bvalues'] + for name, data in all_data.items(): + for algorithm in algorithms: + algorithm_dict = algorithm_information.get(algorithm, {}) + xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}), + "strict": algorithm_dict.get("xfail_names", {}).get(name, True)} + kwargs = algorithm_dict.get("options", {}) + tolerances = algorithm_dict.get("tolerances", {}) + requires_matlab = algorithm_dict.get("requires_matlab", False) + yield name, bvals, data, algorithm, xfail, kwargs, tolerances, requires_matlab \ No newline at end of file diff --git a/pytest.ini b/pytest.ini index e61feaf..ed02be7 100644 --- a/pytest.ini +++ b/pytest.ini @@ -2,4 +2,5 @@ markers = slow: marks tests as slow (deselect with '-m "not slow"') addopts = - -m 'not slow' \ No newline at end of file + -m 'not slow' +testpaths = tests \ No newline at end of file diff --git a/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py new file mode 100644 index 0000000..455e639 --- /dev/null +++ b/src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py @@ -0,0 +1,120 @@ +from src.wrappers.OsipiBase import OsipiBase +import numpy as np +import matlab.engine + + +class ASD_MemorialSloanKettering_QAMPER_IVIM(OsipiBase): + """ + Bi-exponential fitting algorithm by Oliver Gurney-Champion, Amsterdam UMC + """ + + # 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 = "LoCastro, Dr. Ramesh Paudyal, Dr. Amita Shukla-Dave" + id_algorithm_type = "Bi-exponential fit" + id_return_parameters = "f, D*, D, S0" + 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? + + # Supported inputs in the standardized class + supported_bounds = True + supported_initial_guess = True + supported_thresholds = False + + def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, eng=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(OGC_AmsterdamUMC_biexp, self).__init__(bvalues, bounds, initial_guess, fitS0) + super(ASD_MemorialSloanKettering_QAMPER_IVIM, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess) + self.initialize(bounds, initial_guess) + if eng is None: + print('initiating matlab; this may take some time. For repeated testing one could use the optional input eng as an already initiated matlab engine') + self.eng=matlab.engine.start_matlab() + self.keep_alive=False + else: + self.eng = eng + self.keep_alive=True + + def algorithm(self,dwi_arr, bval_arr, LB0, UB0, x0in): + dwi_arr = matlab.double(dwi_arr.tolist()) + bval_arr = matlab.double(bval_arr.tolist()) + LB0 = matlab.double(LB0.tolist()) + UB0 = matlab.double(UB0.tolist()) + x0in = matlab.double(x0in.tolist()) + results = self.eng.IVIM_standard_bcin( + dwi_arr, bval_arr, 0.0, LB0, UB0, x0in, False, 0, 0,nargout=11) + (f_arr, D_arr, Dx_arr, s0_arr, fitted_dwi_arr, RSS, rms_val, chi, AIC, BIC, R_sq) = results + return D_arr/1000, f_arr, Dx_arr/1000, s0_arr + + def initialize(self, bounds, initial_guess): + if bounds is None: + print('warning, no bounds were defined, so algorithm-specific default bounds are used') + self.bounds=([1e-6, 0, 0.004, 0],[0.003, 1.0, 0.2, 5]) + else: + self.bounds=bounds + if initial_guess is None: + print('warning, no initial guesses were defined, so algorithm-specific default initial guess is used') + self.initial_guess = [0.001, 0.2, 0.01, 1] + else: + self.initial_guess = initial_guess + self.use_initial_guess = True + self.use_initial_guess = True + self.use_bounds = True + + def ivim_fit(self, signals, bvalues, **kwargs): + """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) + LB = np.array(self.bounds[0])[[1,0,2,3]] + UB = np.array(self.bounds[1])[[1,0,2,3]] + + fit_results = self.algorithm(np.array(signals)[:,np.newaxis], bvalues, LB, UB, np.array(self.initial_guess)[[1,0,2,3]]) + + results = {} + results["D"] = fit_results[0] + results["f"] = fit_results[1] + results["Dp"] = fit_results[2] + + return results + + def clean(self): + if not self.keep_alive: + if hasattr(self, "eng") and self.eng: + try: + self.eng.quit() + except Exception as e: + print(f"Warning: Failed to quit MATLAB engine cleanly: {e}") + + def __del__(self): + self.clean() + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + self.clean() diff --git a/src/standardized/OGC_AmsterdamUMC_biexp.py b/src/standardized/OGC_AmsterdamUMC_biexp.py index 49a16b2..5afd0cb 100644 --- a/src/standardized/OGC_AmsterdamUMC_biexp.py +++ b/src/standardized/OGC_AmsterdamUMC_biexp.py @@ -55,7 +55,7 @@ def initialize(self, bounds, initial_guess, fitS0): self.bounds=bounds if initial_guess is None: print('warning, no initial guesses were defined, so default bounds are used of [0.001, 0.001, 0.01, 1]') - self.initial_guess = [0.001, 0.001, 0.01, 1] + self.initial_guess = [0.001, 0.1, 0.01, 1] else: self.initial_guess = initial_guess self.use_initial_guess = True diff --git a/tests/IVIMmodels/unit_tests/algorithms.json b/tests/IVIMmodels/unit_tests/algorithms.json index 51a5650..fa90ce5 100644 --- a/tests/IVIMmodels/unit_tests/algorithms.json +++ b/tests/IVIMmodels/unit_tests/algorithms.json @@ -1,5 +1,6 @@ { "algorithms": [ + "ASD_MemorialSloanKettering_QAMPER_IVIM", "ETP_SRI_LinearFitting", "IAR_LU_biexp", "IAR_LU_modified_mix", @@ -14,6 +15,9 @@ "PvH_KB_NKI_IVIMfit", "OJ_GU_seg" ], + "ASD_MemorialSloanKettering_QAMPER_IVIM": { + "requires_matlab": true + }, "IAR_LU_biexp": { "fail_first_time": true }, diff --git a/tests/IVIMmodels/unit_tests/test_ivim_fit.py b/tests/IVIMmodels/unit_tests/test_ivim_fit.py index fb58d13..ea9b336 100644 --- a/tests/IVIMmodels/unit_tests/test_ivim_fit.py +++ b/tests/IVIMmodels/unit_tests/test_ivim_fit.py @@ -7,6 +7,7 @@ from src.wrappers.OsipiBase import OsipiBase #run using python -m pytest from the root folder + def signal_helper(signal): signal = np.asarray(signal) #signal[signal < 0] = 0.00001 @@ -29,45 +30,22 @@ def tolerances_helper(tolerances, data): tolerances["atol"] = tolerances.get("atol", {"f": 2e-1, "D": 5e-4, "Dp": 4e-2}) return tolerances -def data_ivim_fit_saved(): - # Find the algorithms from algorithms.json - file = pathlib.Path(__file__) - algorithm_path = file.with_name('algorithms.json') - with algorithm_path.open() as f: - algorithm_information = json.load(f) - - # Load generic test data generated from the included phantom: phantoms/MR_XCAT_qMRI - generic = file.with_name('generic.json') - with generic.open() as f: - all_data = json.load(f) - - algorithms = algorithm_information["algorithms"] - bvals = all_data.pop('config') - bvals = bvals['bvalues'] - first=True - for name, data in all_data.items(): - for algorithm in algorithms: - algorithm_dict = algorithm_information.get(algorithm, {}) - xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}), - "strict": algorithm_dict.get("xfail_names", {}).get(name, True)} - kwargs = algorithm_dict.get("options", {}) - tolerances = algorithm_dict.get("tolerances", {}) - skiptime=False - if first == True: - if algorithm_dict.get("fail_first_time", {}) == True: - skiptime = True - first = False - yield name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime - -@pytest.mark.parametrize("name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime", data_ivim_fit_saved()) -def test_ivim_fit_saved(name, bvals, data, algorithm, xfail, kwargs, tolerances,skiptime, request, record_property): +def test_ivim_fit_saved(data_ivim_fit_saved, eng, request, record_property): + name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime, requires_matlab = data_ivim_fit_saved + max_time = 0.5 + if requires_matlab: + max_time = 2 + if eng is None: + pytest.skip(reason="Running without matlab; if Matlab is available please run pytest --withmatlab") + else: + kwargs = {**kwargs, 'eng': eng} if xfail["xfail"]: mark = pytest.mark.xfail(reason="xfail", strict=xfail["strict"]) request.node.add_marker(mark) signal = signal_helper(data["data"]) tolerances = tolerances_helper(tolerances, data) - start_time = time.time() # Record the start time fit = OsipiBase(algorithm=algorithm, **kwargs) + start_time = time.time() # Record the start time fit_result = fit.osipi_fit(signal, bvals) elapsed_time = time.time() - start_time # Calculate elapsed time def to_list_if_needed(value): @@ -92,22 +70,18 @@ def to_list_if_needed(value): npt.assert_allclose(fit_result['Dp'],data['Dp'], rtol=tolerances["rtol"]["Dp"], atol=tolerances["atol"]["Dp"]) #assert fit_result['D'] < fit_result['Dp'], f"D {fit_result['D']} is larger than D* {fit_result['Dp']} for {name}" if not skiptime: - assert elapsed_time < 0.5, f"Algorithm {name} took {elapsed_time} seconds, which is longer than 2 second to fit per voxel" #less than 0.5 seconds per voxel - - -def algorithms(): - # Find the algorithms from algorithms.json - file = pathlib.Path(__file__) - algorithm_path = file.with_name('algorithms.json') - with algorithm_path.open() as f: - algorithm_information = json.load(f) - algorithms = algorithm_information["algorithms"] - for algorithm in algorithms: - yield algorithm - -@pytest.mark.parametrize("algorithm", algorithms()) -def test_default_bounds_and_initial_guesses(algorithm): - fit = OsipiBase(algorithm=algorithm) + assert elapsed_time < max_time, f"Algorithm {name} took {elapsed_time} seconds, which is longer than 2 second to fit per voxel" #less than 0.5 seconds per voxel + +def test_default_bounds_and_initial_guesses(algorithmlist,eng): + algorithm, requires_matlab = algorithmlist + if requires_matlab: + if eng is None: + pytest.skip(reason="Running without matlab; if Matlab is available please run pytest --withmatlab") + else: + kwargs = {'eng': eng} + else: + kwargs={} + fit = OsipiBase(algorithm=algorithm,**kwargs) #assert fit.bounds is not None, f"For {algorithm}, there is no default fit boundary" #assert fit.initial_guess is not None, f"For {algorithm}, there is no default fit initial guess" if fit.use_bounds: @@ -127,33 +101,13 @@ def test_default_bounds_and_initial_guesses(algorithm): assert 0.9 <= fit.initial_guess[3] <= 1.1, f"For {algorithm}, the default initial guess for S {fit.initial_guess[3]} is unrealistic; note signal is normalized" -def bound_input(): - # Find the algorithms from algorithms.json - file = pathlib.Path(__file__) - algorithm_path = file.with_name('algorithms.json') - with algorithm_path.open() as f: - algorithm_information = json.load(f) - - # Load generic test data generated from the included phantom: phantoms/MR_XCAT_qMRI - generic = file.with_name('generic.json') - with generic.open() as f: - all_data = json.load(f) - - algorithms = algorithm_information["algorithms"] - bvals = all_data.pop('config') - bvals = bvals['bvalues'] - for name, data in all_data.items(): - for algorithm in algorithms: - algorithm_dict = algorithm_information.get(algorithm, {}) - xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}), - "strict": algorithm_dict.get("xfail_names", {}).get(name, True)} - kwargs = algorithm_dict.get("options", {}) - tolerances = algorithm_dict.get("tolerances", {}) - yield name, bvals, data, algorithm, xfail, kwargs, tolerances - - -@pytest.mark.parametrize("name, bvals, data, algorithm, xfail, kwargs, tolerances", bound_input()) -def test_bounds(name, bvals, data, algorithm, xfail, kwargs, tolerances, request): +def test_bounds(bound_input, eng): + name, bvals, data, algorithm, xfail, kwargs, tolerances, requires_matlab = bound_input + if requires_matlab: + if eng is None: + pytest.skip(reason="Running without matlab; if Matlab is available please run pytest --withmatlab") + else: + kwargs = {**kwargs, 'eng': eng} bounds = ([0.0008, 0.2, 0.01, 1.1], [0.0012, 0.3, 0.02, 1.3]) # deliberately have silly bounds to see whether they are used fit = OsipiBase(algorithm=algorithm, bounds=bounds, initial_guess = [0.001, 0.25, 0.015, 1.2], **kwargs) diff --git a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py index 2074132..6f2cb7c 100644 --- a/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py +++ b/tests/IVIMmodels/unit_tests/test_ivim_synthetic.py @@ -14,8 +14,11 @@ #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, use_prior): +def test_generated(algorithmlist, skip_list, ivim_data, SNR, rtol, atol, fit_count, rician_noise, save_file, save_duration_file, use_prior): # assert save_file == "test" + ivim_algorithm, requires_matlab = algorithmlist + if requires_matlab and eng is None: + pytest.skip(reason="Running without matlab; if Matlab is available please run pytest --withmatlab") rng = np.random.RandomState(42) # random.seed(42) S0 = 1