bounds test automatically check whether bounds are relevant

Author: oliverchampion

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -124,29 +124,30 @@ def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cu
     :return Dp: Fitted Dp
     :return S0: Fitted S0
     """
-    p0 = [p0[0] * 1000, p0[1] * 10, p0[2] * 10, p0[3]]
     try:
         # determine high b-values and data for D
+        dw_data=dw_data/np.mean(dw_data[bvalues==0])
         high_b = bvalues[bvalues >= cutoff]
         high_dw_data = dw_data[bvalues >= cutoff]
         # correct the bounds. Note that S0 bounds determine the max and min of f
-        bounds1 = ([bounds[0][0] * 1000., 0.7 - bounds[1][1]], [bounds[1][0] * 1000., 1.3 - bounds[0][
-            1]])  # By bounding S0 like this, we effectively insert the boundaries of f
-        # fit for S0' and D
-        params, _ = curve_fit(lambda b, Dt, int: int * np.exp(-b * Dt / 1000), high_b, high_dw_data,
-                              p0=(p0[0], p0[3]-p0[1]/10),
+        bounds1 = ([bounds[0][0], 0], [bounds[1][0], 10000000000])
+        params, _ = curve_fit(lambda b, Dt, int: int * np.exp(-b * Dt ), high_b, high_dw_data,
+                              p0=(p0[0], p0[3]-p0[1]),
                               bounds=bounds1)
-        Dt, Fp = params[0] / 1000, 1 - params[1]
+        Dt, Fp = 0+params[0], 1 - params[1]
+        if Fp < bounds[0][1] : Fp = bounds[0][1]
+        if Fp > bounds[1][1] : Fp = bounds[1][1]
+
         # remove the diffusion part to only keep the pseudo-diffusion
         dw_data_remaining = dw_data - (1 - Fp) * np.exp(-bvalues * Dt)
-        bounds2 = (bounds[0][2]*10, bounds[1][2]*10)
+        bounds2 = (bounds[0][2], bounds[1][2])
         # fit for D*
         params, _ = curve_fit(lambda b, Dp: Fp * np.exp(-b * Dp), bvalues, dw_data_remaining, p0=(p0[2]), bounds=bounds2)
-        Dp = params[0]
-        return Dt, Fp, Dp
+        Dp = 0+params[0]
+        return Dt, np.float64(Fp), Dp
     except:
         # if fit fails, return zeros
-        # print('segnetned fit failed')
+        # print('segmented fit failed')
         return 0., 0., 0.
 
 
@@ -235,17 +236,17 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
     try:
         if not fitS0:
             # bounds are rescaled such that each parameter changes at roughly the same rate to help fitting.
-            bounds = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10],
+            bounds2 = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10],
                       [bounds[1][0] * 1000, bounds[1][1] * 10, bounds[1][2] * 10])
             p0=[p0[0]*1000,p0[1]*10,p0[2]*10]
-            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=p0, bounds=bounds)
+            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=p0, bounds=bounds2)
             S0 = 1
         else:
             # bounds are rescaled such that each parameter changes at roughly the same rate to help fitting.
-            bounds = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10, bounds[0][3]],
+            bounds2 = ([bounds[0][0] * 1000, bounds[0][1] * 10, bounds[0][2] * 10, bounds[0][3]],
                       [bounds[1][0] * 1000, bounds[1][1] * 10, bounds[1][2] * 10, bounds[1][3]])
             p0=[p0[0]*1000,p0[1]*10,p0[2]*10,p0[3]]
-            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p0, bounds=bounds)
+            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p0, bounds=bounds2)
             S0 = params[3]
         # correct for the rescaling of parameters
         Dt, Fp, Dp = params[0] / 1000, params[1] / 10, params[2] / 10
@@ -261,7 +262,7 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
             Dt, Fp, Dp = fit_segmented(bvalues, dw_data, bounds=bounds)
             return Dt, Fp, Dp, 1
         else:
-            return fit_segmented(bvalues, dw_data)
+            return fit_segmented(bvalues, dw_data, bounds=bounds)
 
 
 def fit_least_squares_array_tri_exp(bvalues, dw_data, S0_output=True, fitS0=True, njobs=4,
@@ -561,19 +562,19 @@ def neg_log_prior(p):
             Dt, Fp, Dp = p[0], p[1], p[2]
         # make D*<D very unlikely
         if (Dp < Dt):
-            return 1e3
+            return 1e10
         else:
             # determine and return the prior for D, f and D* (and S0)
             if len(p) == 4:
-                if Dt_range[0] < Dt < Dt_range[1] and Fp_range[0] < Fp < Fp_range[1] and Dp_range[0] < Dp < Dp_range[1]: # and S0_range[0] < S0 < S0_range[1]: << not sure whether this helps. Technically it should be here
+                if Dt_range[0] < Dt < Dt_range[1] and Fp_range[0] < Fp < Fp_range[1] and Dp_range[0] < Dp < Dp_range[1] and S0_range[0] < S0 < S0_range[1]: #<< not sure whether this helps. Technically it should be here
                     return 0
                 else:
-                    return 1e3
+                    return 1e10
             else:
                 if Dt_range[0] < Dt < Dt_range[1] and Fp_range[0] < Fp < Fp_range[1] and Dp_range[0] < Dp < Dp_range[1]:
                     return 0
                 else:
-                    return 1e3
+                    return 1e10
 
     return neg_log_prior
 
@@ -638,7 +639,7 @@ def parfun(i):
     return Dt_pred, Fp_pred, Dp_pred, S0_pred
 
 
-def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS0=True):
+def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS0=True, bounds=([0,0,0,0],[0.005,1.5,2,2.5])):
     '''
     This is an implementation of the Bayesian IVIM fit. It returns the Maximum a posterior probability.
     The fit is taken from Barbieri et al. which was initially introduced in http://arxiv.org/10.1002/mrm.25765 and
@@ -655,7 +656,7 @@ def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS
     '''
     try:
         # define fit bounds
-        bounds = [(0, 0.005), (0, 1.5), (0, 2), (0, 2.5)]
+        bounds = [(bounds[0][0], bounds[1][0]), (bounds[0][1], bounds[1][1]), (bounds[0][2], bounds[1][2]), (bounds[0][3], bounds[1][3])]
         # Find the Maximum a posterior probability (MAP) by minimising the negative log of the posterior
         if fitS0:
             params = minimize(neg_log_posterior, x0=x0, args=(bvalues, dw_data, neg_log_prior), bounds=bounds)

src/standardized/ETP_SRI_LinearFitting.py

@@ -37,7 +37,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         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
         # Could be a good idea to have all the submission-specfic variable be 
         # defined with initials?
         self.ETP_weighting = weighting

src/standardized/IAR_LU_biexp.py

@@ -37,7 +37,10 @@ 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
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/IAR_LU_modified_mix.py

@@ -37,7 +37,10 @@ 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
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/IAR_LU_modified_topopro.py

@@ -37,7 +37,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_modified_topopro, 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
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/IAR_LU_segmented_2step.py

@@ -37,7 +37,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_segmented_2step, 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
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/IAR_LU_segmented_3step.py

@@ -37,7 +37,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_segmented_3step, 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
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/IAR_LU_subtracted.py

@@ -37,7 +37,10 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
             the requirements.
         """
         super(IAR_LU_subtracted, 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
         # Check the inputs
 
         # Initialize the algorithm

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -28,7 +28,7 @@ class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
     accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
     accepts_priors = True
 
-    def __init__(self, bvalues=None, thresholds=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, prior_in=None):
+    def __init__(self, bvalues=None, bounds=None, thresholds = None, initial_guess=None, fitS0=True, prior_in=None):
         """
             Everything this algorithm requires should be implemented here.
             Number of segmentation thresholds, bounds, etc.
@@ -42,23 +42,28 @@ def __init__(self, bvalues=None, thresholds=None, bounds=([0, 0, 0.005, 0.7],[0.
         super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues, bounds, initial_guess) #, fitS0, prior_in)
         self.OGC_algorithm = fit_bayesian
         self.initialize(bounds, initial_guess, fitS0, prior_in)
+        self.fit_segmented=fit_segmented
 
-    def initialize(self, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, prior_in=None):
+    def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
         if bounds is None:
-            self.bounds=([0, 0, 0.005, 0.7],[0.005, 1, 0.2, 1.3])
+            self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.1, 1.3])
         else:
             self.bounds=bounds
+            self.use_bounds = True
+        if initial_guess is None:
+            self.initial_guess = [0.001, 0.1, 0.01, 1]
+        else:
+            self.initial_guess = initial_guess
+            self.use_initial_guess = True
         if prior_in is None:
+            print('using a flat prior between bounds')
             self.neg_log_prior=flat_neg_log_prior([self.bounds[0][0],self.bounds[1][0]],[self.bounds[0][1],self.bounds[1][1]],[self.bounds[0][2],self.bounds[1][2]],[self.bounds[0][3],self.bounds[1][3]])
         else:
+            print('warning, bounds are not used, as a prior is used instead')
             if len(prior_in) is 4:
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2],prior_in[3])
             else:
                 self.neg_log_prior = empirical_neg_log_prior(prior_in[0], prior_in[1], prior_in[2])
-        if initial_guess is None:
-            self.initial_guess = [0.001, 0.5, 0.1, 1]
-        else:
-            self.initial_guess = initial_guess
         self.fitS0=fitS0
 
     def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
@@ -71,12 +76,14 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
         Returns:
             _type_: _description_
         """
-        if initial_guess is not None and len(initial_guess) == 4:
-            self.initial_guess = initial_guess
         bvalues=np.array(bvalues)
-        fit_results = fit_segmented(bvalues, signals, bounds=self.bounds, cutoff=150, p0=self.initial_guess)
-        fit_results=fit_results+(1,)
-        fit_results = self.OGC_algorithm(bvalues, signals, self.neg_log_prior, x0=fit_results, fitS0=self.fitS0)
+        epsilon = 0.000001
+        fit_results = self.fit_segmented(bvalues, signals, bounds=self.bounds, cutoff=150, p0=self.initial_guess)
+        fit_results=np.array(fit_results+(1,))
+        for i in range(4):
+            if fit_results[i] < self.bounds[0][i] : fit_results[0] = self.bounds[0][i]+epsilon
+            if fit_results[i] > self.bounds[1][i] : fit_results[0] = self.bounds[1][i]-epsilon
+        fit_results = self.OGC_algorithm(bvalues, signals, self.neg_log_prior, x0=fit_results, fitS0=self.fitS0, bounds=self.bounds)
 
         results = {}
         results["D"] = fit_results[0]

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -38,20 +38,23 @@ def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=Non
         #super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues, bounds, initial_guess, fitS0)
         super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues=bvalues, bounds=bounds, initial_guess=initial_guess)
         self.OGC_algorithm = fit_least_squares
-        #self.initialize(bounds, initial_guess, fitS0)
         self.fitS0=fitS0
+        self.initialize(bounds, initial_guess, fitS0)
 
     def initialize(self, bounds, initial_guess, fitS0):
         if bounds is None:
             self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3])
         else:
             self.bounds=bounds
+            self.use_bounds = True
         if initial_guess is None:
             self.initial_guess = [0.001, 0.1, 0.01, 1]
         else:
             self.initial_guess = initial_guess
+            self.use_initial_guess = True
         self.fitS0=fitS0
 
+
     def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
         """Perform the IVIM fit