Matlab installation + adding Amita's code to wrapper

.github/workflows/unit_test.yml

@@ -30,6 +30,24 @@ jobs:
         run: |
           python -m pip install --upgrade pip
           pip install -r requirements.txt
+      - name: Check out repository
+        uses: actions/checkout@v4
+      - name: Set environment variables for MATLAB Engine
+	if: ${{ runner.os == 'Linux' || runner.os == 'macOS' }}
+        run: |
+          if [[ "${{ runner.os }}" == "macOS" ]]; then
+            echo "DYLD_LIBRARY_PATH=/Users/runner/hostedtoolcache/MATLAB/2025.1.999/arm64/MATLAB.app/bin/maca64:$DYLD_LIBRARY_PATH" >> $GITHUB_ENV
+          elif [[ "${{ runner.os }}" == "Linux" ]]; then
+            echo "LD_LIBRARY_PATH=/opt/hostedtoolcache/MATLAB/2025.1.999/x64/bin/glnxa64:$LD_LIBRARY_PATH" >> $GITHUB_ENV
+          fi
+      - name: Set up MATLAB
+        uses: matlab-actions/setup-matlab@v2
+        with:
+          release: R2025a 
+      - name: Install MATLAB Engine for Python
+        run: |
+          python -m pip install matlabengine==25.1.2
+
       - name: Test with pytest
         run: |
           pip install pytest pytest-cov

.gitignore

@@ -19,7 +19,7 @@ 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
 
 
 # Unit test / coverage reports

src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py

@@ -0,0 +1,127 @@
+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 __del__(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 __enter__(self):
+        return self
+
+
+    def __exit__(self, exc_type, exc_val, exc_tb):
+        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}")

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

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,10 @@
         "PvH_KB_NKI_IVIMfit",
         "OJ_GU_seg"
     ],
+    "ASD_MemorialSloanKettering_QAMPER_IVIM": {
+        "requieres_matlab": true,
+        "fail_first_time": true
+    },
     "IAR_LU_biexp": {
         "fail_first_time": true
     },

tests/IVIMmodels/unit_tests/test_ivim_fit.py

@@ -6,6 +6,9 @@
 import time
 from src.wrappers.OsipiBase import OsipiBase
 #run using python -m pytest from the root folder
+import matlab.engine
+
+eng=matlab.engine.start_matlab()
 
 def signal_helper(signal):
     signal = np.asarray(signal)
@@ -44,9 +47,9 @@ def data_ivim_fit_saved():
     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:
+    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)}
@@ -56,7 +59,9 @@ def data_ivim_fit_saved():
             if first == True:
                 if algorithm_dict.get("fail_first_time", {}) == True:
                     skiptime = True
-                    first = False
+                    first = False # this will only work for 1 algorithm... If two algorithms use fail_first_time, the second will not see it is the first time.
+            if algorithm_dict.get("requieres_matlab", {}) == True:
+                kwargs={**kwargs,'eng': eng}
             yield name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime
 
 @pytest.mark.parametrize("name, bvals, data, algorithm, xfail, kwargs, tolerances, skiptime", data_ivim_fit_saved())
@@ -103,11 +108,15 @@ def algorithms():
         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)
+        algorithm_dict = algorithm_information.get(algorithm, {})
+        args={}
+        if algorithm_dict.get("requieres_matlab", {}) == True:
+            args['eng'] = eng
+        yield algorithm, args
+
+@pytest.mark.parametrize("algorithm, args", algorithms())
+def test_default_bounds_and_initial_guesses(algorithm, args):
+    fit = OsipiBase(algorithm=algorithm,**args)
     #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:
@@ -149,6 +158,8 @@ def bound_input():
                 "strict": algorithm_dict.get("xfail_names", {}).get(name, True)}
             kwargs = algorithm_dict.get("options", {})
             tolerances = algorithm_dict.get("tolerances", {})
+            if algorithm_dict.get("requieres_matlab", {}) == True:
+                kwargs={**kwargs,'eng': eng}
             yield name, bvals, data, algorithm, xfail, kwargs, tolerances