Added Bayesian fit

With flat prior, as hard to get emperical prior from single input

Author: oliverchampion

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -519,7 +519,37 @@ def neg_log_prior(p):
 
     return neg_log_prior
 
+def flat_neg_log_prior(Dt_range, Fp_range, Dp_range, S0_range=None):
+    """
+    This function determines the negative of the log of the empirical prior probability of the IVIM parameters
+    :param Dt0: 1D Array with the initial D estimates
+    :param Dt0: 1D Array with the initial f estimates
+    :param Dt0: 1D Array with the initial D* estimates
+    :param Dt0: 1D Array with the initial S0 estimates (optional)
+    """
+    def neg_log_prior(p):
+        # depends on whether S0 is fitted or not
+        if len(p) == 4:
+            Dt, Fp, Dp, S0 = p[0], p[1], p[2], p[3]
+        else:
+            Dt, Fp, Dp = p[0], p[1], p[2]
+        # make D*<D very unlikely
+        if (Dp < Dt):
+            return 1e8
+        else:
+            eps = 1e-8
+            Dp_prior = stats.loguniform.pdf(Dp, Dp_range[0], scale=Dp_range[1]-Dp_range[0])
+            Dt_prior = stats.loguniform.pdf(Dt, Dt_range[0], scale=Dt_range[1]-Dt_range[0])
+            Fp_prior = stats.loguniform.pdf(Fp, Fp_range[0], scale=Fp_range[1]-Fp_range[0])
+            # determine and return the prior for D, f and D* (and S0)
+            if len(p) == 4:
+                S0_prior = stats.loguniform.pdf(S0, S0_range[0], scale=S0_range[1] - S0_range[0])
+                return -np.log(Dp_prior + eps) - np.log(Dt_prior + eps) - np.log(Fp_prior + eps) - np.log(
+                    S0_prior + eps)
+            else:
+                return -np.log(Dp_prior + eps) - np.log(Dt_prior + eps) - np.log(Fp_prior + eps)
 
+    return neg_log_prior
 def neg_log_likelihood(p, bvalues, dw_data):
     """
     This function determines the negative of the log of the likelihood of parameters p, given the data dw_data for the Bayesian fit

src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -0,0 +1,66 @@
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import flat_neg_log_prior, fit_bayesian
+import numpy as np
+
+class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
+    """
+    Bayesian 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 = "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, bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3]), initial_guess=None, fitS0=True, thresholds=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_Bayesian_biexp, self).__init__(bvalues, bounds,initial_guess,fitS0)
+        if bounds is None:
+            self.bounds=([0, 0, 0.005, 0.7],[0.005, 0.7, 0.2, 1.3])
+        else:
+            self.bounds=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][1],self.bounds[1][1]],[self.bounds[0][2],self.bounds[1][2]])
+        self.OGC_algorithm = fit_bayesian
+        self.bounds=bounds
+        self.initial_guess=initial_guess
+        self.fitS0=fitS0
+
+    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_
+        """
+        bvalues=np.array(bvalues)
+        fit_results = self.OGC_algorithm(bvalues, signals, self.neg_log_prior, x0=self.initial_guess, fitS0=self.fitS0)
+
+        D = fit_results[0]
+        f = fit_results[1]
+        Dstar = fit_results[2]
+
+        return f, Dstar, D