Added bi-exp and segmented fit to wrapper

minor changes to origin by author

Author: oliverchampion

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -150,7 +150,7 @@ def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cu
 
 
 def fit_least_squares_array(bvalues, dw_data, S0_output=True, fitS0=True, njobs=4,
-                            bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])):
+                            bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]),p0=[1, 1, 0.1, 1]):
     """
     This is an implementation of the conventional IVIM fit. It is fitted in array form.
     :param bvalues: 1D Array with the b-values
@@ -175,7 +175,7 @@ def fit_least_squares_array(bvalues, dw_data, S0_output=True, fitS0=True, njobs=
             try:
                 # defining parallel function
                 def parfun(i):
-                    return fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0, bounds=bounds)
+                    return fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0, bounds=bounds,p0=p0)
 
                 output = Parallel(n_jobs=njobs)(delayed(parfun)(i) for i in tqdm(range(len(dw_data)), position=0, leave=True))
                 Dt, Fp, Dp, S0 = np.transpose(output)
@@ -192,14 +192,14 @@ def parfun(i):
             # running in a single loop and filling arrays
             for i in tqdm(range(len(dw_data)), position=0, leave=True):
                 Dt[i], Fp[i], Dp[i], S0[i] = fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0,
-                                                               bounds=bounds)
+                                                               bounds=bounds,p0=p0)
         return [Dt, Fp, Dp, S0]
     else:
         # if S0 is not exported
         if njobs > 1:
             try:
                 def parfun(i):
-                    return fit_least_squares(bvalues, dw_data[i, :], fitS0=fitS0, bounds=bounds)
+                    return fit_least_squares(bvalues, dw_data[i, :], fitS0=fitS0, bounds=bounds,p0=p0)
 
                 output = Parallel(n_jobs=njobs)(delayed(parfun)(i) for i in tqdm(range(len(dw_data)), position=0, leave=True))
                 Dt, Fp, Dp = np.transpose(output)
@@ -213,12 +213,12 @@ def parfun(i):
             Fp = np.zeros(len(dw_data))
             for i in range(len(dw_data)):
                 Dt[i], Fp[i], Dp[i] = fit_least_squares(bvalues, dw_data[i, :], S0_output=S0_output, fitS0=fitS0,
-                                                        bounds=bounds)
+                                                        bounds=bounds,p0=p0)
         return [Dt, Fp, Dp]
 
 
 def fit_least_squares(bvalues, dw_data, S0_output=False, fitS0=True,
-                      bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])):
+                      bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]),p0=[1, 1, 0.1, 1]):
     """
     This is an implementation of the conventional IVIM fit. It fits a single curve
     :param bvalues: Array with the b-values
@@ -242,7 +242,7 @@ 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.
             bounds = ([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]])
-            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=[1, 1, 0.1, 1], bounds=bounds)
+            params, _ = curve_fit(ivimN, bvalues, dw_data, p0=p0, bounds=bounds)
             S0 = params[3]
         # correct for the rescaling of parameters
         Dt, Fp, Dp = params[0] / 1000, params[1] / 10, params[2] / 10
@@ -563,7 +563,6 @@ def fit_bayesian_array(bvalues, dw_data, paramslsq, arg):
     :return Dp: Array with Dp in each voxel
     :return S0: Array with S0 in each voxel
     """
-    arg = checkarg_lsq(arg)
     # fill out missing args
     Dt0, Fp0, Dp0, S00 = paramslsq
     # determine prior

src/standardized/OGC_AmsterdamUMC_biexp.py

@@ -24,19 +24,20 @@ class OGC_AmsterdamUMC_biexp(OsipiBase):
     required_bounds = False
     required_bounds_optional = True  # Bounds may not be required but are optional
     required_initial_guess = False
-    required_initial_guess_optional = True
+    required_initial_guess_optional = False
     accepted_dimensions = 1  # Not sure how to define this for the number of accepted dimensions. Perhaps like the thresholds, at least and at most?
 
-    def __init__(self, bvalues, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])):
+    def __init__(self, bvalues=None, thershold=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])):
         """
             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, None, bounds, None)
+        super(OGC_AmsterdamUMC_biexp, self).__init__(bvalues, bounds)
         self.OGC_algorithm = fit_least_squares_array
+        self.bounds=bounds
 
     def ivim_fit(self, signals, bvalues=None):
         """Perform the IVIM fit
@@ -49,7 +50,7 @@ def ivim_fit(self, signals, bvalues=None):
             _type_: _description_
         """
 
-        fit_results = self.OGC_algorithm(bvalues, signals)
+        fit_results = self.OGC_algorithm(bvalues, signals,cutoff=self.thershold, bounds=self.bounds)
 
         D = fit_results[0]
         f = fit_results[1]

src/standardized/OGC_AmsterdamUMC_biexp_segmented.py

@@ -0,0 +1,62 @@
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import fit_segmented
+
+
+class OGC_AmsterdamUMC_segbiexp(OsipiBase):
+    """
+    Segmented 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 = "Oliver Gurney Champion, Amsterdam UMC"
+    id_algorithm_type = "Segmented 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?
+
+    def __init__(self, bvalues=None, thershold=None, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True):
+        """
+            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_segbiexp, self).__init__(bvalues, bounds, initial_guess)
+        self.OGC_algorithm = fit_segmented
+        self.bounds=bounds
+        self.initial_guess=initial_guess
+        self.fitS0=fitS0
+        self.thershold=thershold
+
+    def ivim_fit(self, signals, bvalues=None):
+        """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_
+        """
+
+        fit_results = self.OGC_algorithm(bvalues, signals,bounds=self.bounds,p0=self.initial_guess,fitS0=self.fitS0)
+
+        D = fit_results[0]
+        f = fit_results[1]
+        Dstar = fit_results[2]
+
+        return f, Dstar, D