Merge pull request #126 from OSIPI/wrapper_dev

Dictionary bounds + standard OSIPI bounds

Author: IvanARashid

conftest.py

@@ -282,8 +282,8 @@ def bound_input(datafile, algorithms):
         for algorithm in algorithms["algorithms"]:
             algorithm_dict = algorithms.get(algorithm, {})
             if not algorithm_dict.get('deep_learning',False):
-                xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}),
-                    "strict": algorithm_dict.get("xfail_names", {}).get(name, True)}
+                xfail = {"xfail": name in algorithm_dict.get("xfail_names", {}) or "bounds" in algorithm_dict.get("xfail_names", {}),
+                    "strict": algorithm_dict.get("xfail_names", {}).get("bounds", 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)

src/original/IAR_LundUniversity/ivim_fit_method_modified_mix.py

@@ -86,9 +86,9 @@ def __init__(self, gtab, bounds=None, maxiter=10, xtol=1e-8, rescale_units=False
                                     (bounds[0][1]*1000, bounds[1][1]*1000), \
                                     (bounds[0][2]*1000, bounds[1][2]*1000)])
         else: # Finally, if units if µm2/ms are already used
-            self.bounds = np.array([(bounds[0][0], bounds[1][0], \
+            self.bounds = np.array([(bounds[0][0], bounds[1][0]), \
                                     (bounds[0][1], bounds[1][1]), \
-                                    (bounds[0][2], bounds[1][2]))])
+                                    (bounds[0][2], bounds[1][2])])
 
     @multi_voxel_fit
     def fit(self, data, bounds_de=None):

src/original/IAR_LundUniversity/ivim_fit_method_modified_topopro.py

@@ -84,9 +84,9 @@ def __init__(self, gtab, bounds=[[0, 0.005, 1e-5], [1, 0.1, 0.004]], \
                                     (bounds[0][1]*1000, bounds[1][1]*1000), \
                                     (bounds[0][2]*1000, bounds[1][2]*1000)])
         else: # Finally, if units if µm2/ms are already used
-            self.bounds = np.array([(bounds[0][0], bounds[1][0], \
+            self.bounds = np.array([(bounds[0][0], bounds[1][0]), \
                                     (bounds[0][1], bounds[1][1]), \
-                                    (bounds[0][2], bounds[1][2]))])
+                                    (bounds[0][2], bounds[1][2])])
 
     @multi_voxel_fit
     def fit(self, data):

src/original/IAR_LundUniversity/ivim_fit_method_segmented_2step.py

@@ -50,7 +50,7 @@ def fit(self, data):
             data = data / data_max
 
         ### Fit the diffusion signal to bvals >= diff_b_threshold_lower
-        diff_bounds = [(self.bounds[0][0], self.bounds[0][3]), \
+        diff_bounds = [(0, self.bounds[0][3]), \
             (self.bounds[1][0], self.bounds[1][3])] # Bounds for S0 and D
 
         diff_bval_indices = np.where(self.bvals >= self.diff_b_threshold_lower)[0]

src/original/IAR_LundUniversity/ivim_fit_method_segmented_3step.py

@@ -54,7 +54,7 @@ def fit(self, data):
             data = data / data_max
 
         ### Fit the diffusion signal to bvals >= diff_b_threshold_lower
-        diff_bounds = [(self.bounds[0][0], self.bounds[0][3]), \
+        diff_bounds = [(0, self.bounds[0][3]), \
             (self.bounds[1][0], self.bounds[1][3])] # Bounds for S0 and D
 
         diff_bval_indices = np.where(self.bvals >= self.diff_b_threshold_lower)[0]
@@ -80,7 +80,8 @@ def fit(self, data):
         f_est = S0_perf_est/(S0_perf_est + S0_diff_est)
 
         # Fit to the full bi-exponential, f estimate as initial guess, D fixed
-        full_initial_guess = np.array([self.initial_guess[0], f_est, self.initial_guess[2]])
+        f_intial_guess = np.min((f_est, self.bounds[0][1])) if f_est > self.bounds[0][1] else np.max((f_est, self.bounds[1][1]))
+        full_initial_guess = np.array([self.initial_guess[0], f_intial_guess, self.initial_guess[2]])
 
         full_bounds_lower = self.bounds[0][:-1]
         full_bounds_upper = self.bounds[1][:-1]

src/original/IAR_LundUniversity/ivim_fit_method_subtracted.py

@@ -52,7 +52,7 @@ def fit(self, data):
             data = data / data_max
 
         ### Fit the diffusion signal to bvals >= diff_b_threshold_lower
-        diff_bounds = [(self.bounds[0][0], self.bounds[0][3]), \
+        diff_bounds = [(0, self.bounds[0][3]), \
             (self.bounds[1][0], self.bounds[1][3])] # Bounds for S0 and D
 
         diff_bval_indices = np.where(self.bvals >= self.diff_b_threshold_lower)[0]
@@ -64,7 +64,7 @@ def fit(self, data):
 
 
         ### Fit the perfusion signal to bvals <= perf_b_threshold_upper
-        perf_bounds = [(self.bounds[0][0], self.bounds[0][2]), \
+        perf_bounds = [(0, self.bounds[0][2]), \
             (self.bounds[1][0], self.bounds[1][2])] # Bounds for S0 and D*
 
         perf_bvals = self.bvals[self.bvals <= self.perf_b_threshold_upper]

src/standardized/ASD_MemorialSloanKettering_QAMPER_IVIM.py

@@ -41,9 +41,11 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             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)
+
+        self.use_bounds = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
+        self.use_initial_guess = {"f" : True, "D" : True, "Dp" : True, "S0" : True}
+
         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()
@@ -63,21 +65,6 @@ def algorithm(self,dwi_arr, bval_arr, LB0, UB0, x0in):
         (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
 
@@ -88,12 +75,16 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         Returns:
             _type_: _description_
         """
+        bounds = ([self.bounds["D"][0], self.bounds["f"][0], self.bounds["Dp"][0], self.bounds["S0"][0]],
+                  [self.bounds["D"][1], self.bounds["f"][1], self.bounds["Dp"][1], self.bounds["S0"][1]])
+
+        initial_guess = [self.initial_guess["D"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["S0"]]
 
         bvalues=np.array(bvalues)
-        LB = np.array(self.bounds[0])[[1,0,2,3]]
-        UB = np.array(self.bounds[1])[[1,0,2,3]]
+        LB = np.array(bounds[0])[[1,0,2,3]]
+        UB = np.array(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]])
+        fit_results = self.algorithm(np.array(signals)[:,np.newaxis], bvalues, LB, UB, np.array(initial_guess)[[1,0,2,3]])
 
         results = {}
         results["D"] = fit_results[0]

src/standardized/ETP_SRI_LinearFitting.py

@@ -48,8 +48,8 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         super(ETP_SRI_LinearFitting, self).__init__(bvalues, thresholds, bounds, initial_guess)
         if bounds is not None:
             print('warning, bounds from wrapper are not (yet) used in this algorithm')
-        self.use_bounds = False
-        self.use_initial_guess = False
+        self.use_bounds = {"f": False, "Dp": False, "D": False, "S0": False}
+        self.use_initial_guess = {"f": False, "Dp": False, "D": False, "S0": False}
 
         # Could be a good idea to have all the submission-specfic variable be 
         # defined with initials?

src/standardized/IAR_LU_biexp.py

@@ -43,12 +43,18 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_biexp, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        if bounds is not None:
-            print('warning, bounds from wrapper are not (yet) used in this algorithm')
-        self.use_bounds = False
-        self.use_initial_guess = False
+        if bounds is None:
+            self.use_bounds = {"f": False, "Dp": False, "D": False}
+        else:
+            self.use_bounds = {"f": True, "Dp": True, "D": True}
+
+        if initial_guess is None:
+            self.use_initial_guess = {"f": False, "Dp": False, "D": False}
+        else:
+            self.use_initial_guess = {"f": True, "Dp": True, "D": True}
+
         # Check the inputs
-        
+
         # Initialize the algorithm
         if self.bvalues is not None:
             bvec = np.zeros((self.bvalues.size, 3))
@@ -70,6 +76,11 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         Returns:
             _type_: _description_
         """
+
+        # Make sure bounds and initial guess conform to the algorithm requirements
+        bounds = [[self.bounds["S0"][0], self.bounds["f"][0], self.bounds["Dp"][0], self.bounds["D"][0]], 
+                       [self.bounds["S0"][1], self.bounds["f"][1], self.bounds["Dp"][1], self.bounds["D"][1]]]
+        initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
 
         if self.IAR_algorithm is None:
             if bvalues is None:
@@ -81,7 +92,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
             bvec[:,2] = 1
             gtab = gradient_table(bvalues, bvecs=bvec, b0_threshold=0)
 
-            self.IAR_algorithm = IvimModelBiExp(gtab, bounds=self.bounds, initial_guess=self.initial_guess)
+            self.IAR_algorithm = IvimModelBiExp(gtab, bounds=bounds, initial_guess=initial_guess)
 
         fit_results = self.IAR_algorithm.fit(signals)
 
@@ -103,7 +114,10 @@ def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
         Returns:
             _type_: _description_
         """
-        
+        # Make sure bounds and initial guess conform to the algorithm requirements
+        bounds = [[self.bounds["S0"][0], self.bounds["f"][0], self.bounds["Dp"][0], self.bounds["D"][0]],
+                       [self.bounds["S0"][1], self.bounds["f"][1], self.bounds["Dp"][1], self.bounds["D"][1]]]
+        initial_guess = [self.initial_guess["S0"], self.initial_guess["f"], self.initial_guess["Dp"], self.initial_guess["D"]]
         if self.IAR_algorithm is None:
             if bvalues is None:
                 bvalues = self.bvalues
@@ -114,7 +128,7 @@ def ivim_fit_full_volume(self, signals, bvalues, **kwargs):
             bvec[:,2] = 1
             gtab = gradient_table(bvalues, bvecs=bvec, b0_threshold=0)
 
-            self.IAR_algorithm = IvimModelBiExp(gtab, bounds=self.bounds, initial_guess=self.initial_guess)
+            self.IAR_algorithm = IvimModelBiExp(gtab, bounds=bounds, initial_guess=initial_guess)
         b0_index = np.where(bvalues == 0)[0][0]
         mask = signals[...,b0_index]>0
         fit_results = self.IAR_algorithm.fit(signals, mask=mask)

src/standardized/IAR_LU_modified_mix.py

@@ -45,23 +45,25 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_modified_mix, self).__init__(bvalues, thresholds, bounds, initial_guess)
-        if bounds is not None:
-            print('warning, bounds from wrapper are not (yet) used in this algorithm')
-        self.use_bounds = False
-        self.use_initial_guess = False
+
+        self.use_bounds = {"f": False, "Dp": False, "D": False} # This algorithm performs intermediate steps that generates initial guesses outside very constrainted bounds
+        self.use_initial_guess = {"f": False, "Dp": False, "D": False} # This algorithm does not use initial guesses
 
         # Additional options
         self.stochastic = True
 
         # Check the inputs
-        
+
         # Initialize the algorithm
         if self.bvalues is not None:
             bvec = np.zeros((self.bvalues.size, 3))
             bvec[:,2] = 1
             gtab = gradient_table(self.bvalues, bvec, b0_threshold=0)
 
-            self.IAR_algorithm = IvimModelVP(gtab, bounds=self.bounds, rescale_results_to_mm2_s=True)
+            bounds = [[self.bounds["f"][0], self.bounds["Dp"][0]*1000, self.bounds["D"][0]*1000], 
+                      [self.bounds["f"][1], self.bounds["Dp"][1]*1000, self.bounds["D"][1]*1000]]
+
+            self.IAR_algorithm = IvimModelVP(gtab, bounds=bounds, rescale_units=False, rescale_results_to_mm2_s=True)
         else:
             self.IAR_algorithm = None
 
@@ -76,6 +78,9 @@ def ivim_fit(self, signals, bvalues, **kwargs):
         Returns:
             _type_: _description_
         """
+
+        bounds = [[self.bounds["f"][0], self.bounds["Dp"][0]*1000, self.bounds["D"][0]*1000], 
+                  [self.bounds["f"][1], self.bounds["Dp"][1]*1000, self.bounds["D"][1]*1000]]
 
         if self.IAR_algorithm is None:
             if bvalues is None:
@@ -87,7 +92,7 @@ def ivim_fit(self, signals, bvalues, **kwargs):
             bvec[:,2] = 1
             gtab = gradient_table(bvalues, bvec, b0_threshold=0)
 
-            self.IAR_algorithm = IvimModelVP(gtab, bounds=self.bounds, rescale_results_to_mm2_s=True)
+            self.IAR_algorithm = IvimModelVP(gtab, bounds=bounds, rescale_results_to_mm2_s=True)
 
         fit_results = self.IAR_algorithm.fit(signals)