Automated testing compatible with Wrapper

some work on initial guess testing, but did not work robust.

Author: oliverchampion

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -130,6 +130,8 @@ def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cu
         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
+        assert bounds[0][1] < p0[1]
+        assert bounds[1][1] > p0[1]
         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]),
@@ -238,15 +240,15 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
             # bounds are rescaled such that each parameter changes at roughly the same rate to help fitting.
             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=bounds2)
+            p1=[p0[0]*1000,p0[1]*10,p0[2]*10]
+            params, _ = curve_fit(ivimN_noS0, bvalues, dw_data, p0=p1, bounds=bounds2)
             S0 = 1
         else:
             # bounds are rescaled such that each parameter changes at roughly the same rate to help fitting.
             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=bounds2)
+            p1=[p0[0]*1000,p0[1]*10,p0[2]*10,p0[3]]
+            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p1, bounds=bounds2)
             S0 = params[3]
         # correct for the rescaling of parameters
         Dt, Fp, Dp = params[0] / 1000, params[1] / 10, params[2] / 10
@@ -259,10 +261,10 @@ def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
         # if fit fails, then do a segmented fit instead
         # print('lsq fit failed, trying segmented')
         if S0_output:
-            Dt, Fp, Dp = fit_segmented(bvalues, dw_data, bounds=bounds)
+            Dt, Fp, Dp = fit_segmented(bvalues, dw_data, bounds=bounds,p0=p0)
             return Dt, Fp, Dp, 1
         else:
-            return fit_segmented(bvalues, dw_data, bounds=bounds)
+            return fit_segmented(bvalues, dw_data, bounds=bounds,p0=p0)
 
 
 def fit_least_squares_array_tri_exp(bvalues, dw_data, S0_output=True, fitS0=True, njobs=4,

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -34,7 +34,7 @@ class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
     supported_initial_guess = True
     supported_thresholds = 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=[0.001, 0.1, 0.01, 1], fitS0=True, prior_in=None):
+    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
 
         """
             Everything this algorithm requires should be implemented here.
@@ -51,11 +51,20 @@ def __init__(self, bvalues=None, thresholds=None, bounds=([0, 0, 0.005, 0.7],[0.
         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=[0.001, 0.1, 0.01, 1], fitS0=True, prior_in=None):
-        self.bounds=bounds
-        self.use_bounds = True
-        self.initial_guess = initial_guess
+    def initialize(self, bounds=None, initial_guess=None, fitS0=True, prior_in=None):
+        if bounds is None:
+            print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]')
+            self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3])
+        else:
+            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]
+        else:
+            self.initial_guess = initial_guess
         self.use_initial_guess = True
+        self.use_bounds = True
+        self.thresholds = 150
         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]])
@@ -81,7 +90,7 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
 
         epsilon = 0.000001
         fit_results = fit_segmented(bvalues, signals, bounds=self.bounds, cutoff=self.thresholds, p0=self.initial_guess)
-        fit_results=fit_results+(1,)
+        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

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -33,7 +33,7 @@ class OGC_AmsterdamUMC_biexp(OsipiBase):
     supported_initial_guess = True
     supported_thresholds = False
 
-    def __init__(self, bvalues=None, thresholds=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=[0.001, 0.1, 0.01, 1], fitS0=True):
+    def __init__(self, bvalues=None, thresholds=None, bounds=None, initial_guess=None, fitS0=True):
         """
             Everything this algorithm requires should be implemented here.
             Number of segmentation thresholds, bounds, etc.
@@ -48,14 +48,22 @@ def __init__(self, bvalues=None, thresholds=None, bounds=([0, 0, 0.005, 0.7],[0.
         self.initialize(bounds, initial_guess, fitS0)
 
     def initialize(self, bounds, initial_guess, fitS0):
-        self.bounds=bounds
-        self.use_bounds = True
-        self.initial_guess = initial_guess
-        self.use_initial_guess = True
+        if bounds is None:
+            print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]')
+            self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3])
+        else:
+            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]
+        else:
+            self.initial_guess = initial_guess
+            self.use_initial_guess = True
         self.fitS0=fitS0
+        self.use_initial_guess = True
+        self.use_bounds = True
 
-
-    def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
+    def ivim_fit(self, signals, bvalues, **kwargs):
         """Perform the IVIM fit
 
         Args:
@@ -66,8 +74,6 @@ def ivim_fit(self, signals, bvalues, initial_guess=None, **kwargs):
             _type_: _description_
         """
 
-        if initial_guess is not None and len(initial_guess) == 4:
-            self.initial_guess = initial_guess
         bvalues=np.array(bvalues)
         fit_results = self.OGC_algorithm(bvalues, signals, p0=self.initial_guess, bounds=self.bounds, fitS0=self.fitS0)
 

src/standardized/OGC_AmsterdamUMC_biexp_segmented.py

@@ -33,7 +33,7 @@ class OGC_AmsterdamUMC_biexp_segmented(OsipiBase):
     supported_initial_guess = True
     supported_thresholds = True
 
-    def __init__(self, bvalues=None, thresholds=150, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), initial_guess=[0.001, 0.01, 0.01,1]):
+    def __init__(self, bvalues=None, thresholds=150, bounds=None, initial_guess=None):
         """
             Everything this algorithm requires should be implemented here.
             Number of segmentation thresholds, bounds, etc.
@@ -47,12 +47,23 @@ def __init__(self, bvalues=None, thresholds=150, bounds=([0, 0, 0.005],[0.005, 0
 
 
     def initialize(self, bounds, initial_guess, thresholds):
-        self.bounds=bounds
-        self.use_bounds=True
-        self.initial_guess = initial_guess
-        self.use_initial_guess=True
-        self.thresholds = thresholds
-
+        if bounds is None:
+            print('warning, no bounds were defined, so default bounds are used of [0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3]')
+            self.bounds=([0, 0, 0.005, 0.7],[0.005, 1.0, 0.2, 1.3])
+        else:
+            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]
+        else:
+            self.initial_guess = initial_guess
+        self.use_initial_guess = True
+        self.use_bounds = True
+        if thresholds is None:
+            self.thresholds = 150
+            print('warning, no thresholds were defined, so default bounds are used of  150')
+        else:
+            self.thresholds = thresholds
     def ivim_fit(self, signals, bvalues, **kwargs):
         """Perform the IVIM fit
 

tests/IVIMmodels/unit_tests/test_ivim_fit.py

@@ -6,6 +6,8 @@
 import os
 import time
 
+from exceptiongroup import catch
+
 from src.wrappers.OsipiBase import OsipiBase
 from utilities.data_simulation.GenerateData import GenerateData
 #run using python -m pytest from the root folder
@@ -70,7 +72,7 @@ def test_ivim_fit_saved(name, bvals, data, algorithm, xfail, kwargs, tolerances,
     signal = signal_helper(data["data"])
     tolerances = tolerances_helper(tolerances, data)
     start_time = time.time()  # Record the start time
-    fit = OsipiBase(algorithm=algorithm, **kwargs)
+    fit = OsipiBase(algorithm=algorithm,  **kwargs)
     fit_result = fit.osipi_fit(signal, bvals)
     elapsed_time = time.time() - start_time  # Calculate elapsed time
     def to_list_if_needed(value):
@@ -113,7 +115,7 @@ def test_default_bounds_and_initial_guesses(algorithm):
     fit = OsipiBase(algorithm=algorithm)
     #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.bounds is not None:
+    if fit.use_bounds:
         assert 0 <= fit.bounds[0][0] <= 0.003, f"For {algorithm}, the default lower bound of D {fit.bounds[0][0]} is unrealistic"
         assert 0 <= fit.bounds[1][0] <= 0.01, f"For {algorithm}, the default upper bound of D {fit.bounds[1][0]} is unrealistic"
         assert 0 <= fit.bounds[0][1] <= 1, f"For {algorithm}, the default lower bound of f {fit.bounds[0][1]} is unrealistic"
@@ -123,7 +125,7 @@ def test_default_bounds_and_initial_guesses(algorithm):
         assert 0 <= fit.bounds[0][3] <= 1, f"For {algorithm}, the default lower bound of S {fit.bounds[0][3]} is unrealistic; note data is normaized"
         assert 1 <= fit.bounds[1][3] <= 1000, f"For {algorithm}, the default upper bound of S {fit.bounds[1][3]} is unrealistic; note data is normaized"
         assert fit.bounds[1][0] <= fit.bounds[0][2], f"For {algorithm}, the default upper bound of D {fit.bounds[1][0]} is higher than lower bound of D* {fit.bounds[0][2]}"
-    if fit.initial_guess is not None:
+    if fit.use_initial_guess:
         assert 0.0008 <= fit.initial_guess[0] <= 0.002, f"For {algorithm}, the default initial guess for D {fit.initial_guess[0]} is unrealistic"
         assert 0 <= fit.initial_guess[1] <= 0.5, f"For {algorithm}, the default initial guess for f {fit.initial_guess[1]} is unrealistic"
         assert 0.003 <= fit.initial_guess[2] <= 0.1, f"For {algorithm}, the default initial guess for Ds {fit.initial_guess[2]} is unrealistic"
@@ -169,5 +171,37 @@ def test_bounds(name, bvals, data, algorithm, xfail, kwargs, tolerances, request
         assert bounds[0][2] <= fit_result['Dp'] <= bounds[1][2], f"Result {fit_result['Dp']} out of bounds for data: {name}"
         # S0 is not returned as argument...
         #assert bounds[0][3] <= fit_result['S0'] <= bounds[1][3], f"Result {fit_result['S0']} out of bounds for data: {name}"
+        '''if fit.use_initial_guess:
+            passD=False
+            passf=False
+            passDp=False
+            fit = OsipiBase(algorithm=algorithm, bounds=bounds, initial_guess=[0.0007, 0.25, 0.015, 1.2], **kwargs)
+            try:
+                fit_result = fit.osipi_fit(signal, bvals)
+                if fit_result['D'] == 0:
+                    passD=True
+            except:
+                passD=True
+            fit = OsipiBase(algorithm=algorithm, bounds=bounds, initial_guess=[0.001, 0.19, 0.015, 1.2], **kwargs)
+            try:
+                fit_result = fit.osipi_fit(signal, bvals)
+                if fit_result['f'] == 0:
+                    passf=True
+            except:
+                passf = True
+            try:
+                fit = OsipiBase(algorithm=algorithm, bounds=bounds, initial_guess=[0.001, 0.25, 0.009, 1.2], **kwargs)
+                if fit_result['Dp'] == 0:
+                    passDp = True
+            except:
+                passDp = True
+            assert passD, f"Fit still passes when initial guess D is out of fit bounds; potentially initial guesses not respected for: {name}"
+            assert passf, f"Fit still passes when initial guess f is out of fit bounds; potentially initial guesses not respected for: {name}"
+            assert passDp, f"Fit still passes when initial guess Ds is out of fit bounds; potentially initial guesses not respected for: {name}" '''
+
+
+
+
+