redo bayesian wrapper for single input testing

Author: oliverchampion

.gitignore

@@ -1,8 +1,9 @@
 .idea*
 .github*
 __pycache__/
-
-
+*.nii.gz
+bvals.txt
+*json
 # Unit test / coverage reports
 .tox/
 .coverage

phantoms/MR_XCAT_qMRI/sim_ivim_sig.py

@@ -8,7 +8,7 @@
 # code adapted from MAtlab code by Eric Schrauben: https://github.com/schrau24/XCAT-ERIC
 # This code generates a 4D IVIM phantom as nifti file
 
-def phantom(bvalue, noise, TR=3000, TE=40, motion=False, rician=False, interleaved=False):
+def phantom(bvalue, noise, TR=8000, TE=80, motion=False, rician=False, interleaved=False):
     np.random.seed(42)
     if motion:
         states = range(1,21)
@@ -367,7 +367,8 @@ def XCAT_to_MR_DCE(XCAT, TR, TE, bvalue, D, f, Ds, b0=3, ivim_cont = True):
     bvalue = np.array([0., 1, 2, 5, 10, 20, 30, 50, 75, 100, 150, 250, 350, 400, 550, 700, 850, 1000])
     noise = 0.0005
     motion = False
-    sig, XCAT, Dim, fim, Dpim, legend = phantom(bvalue, noise, motion=motion, interleaved=False)
+    interleaved = False
+    sig, XCAT, Dim, fim, Dpim, legend = phantom(bvalue, noise, motion=motion, interleaved=interleaved)
     # sig = np.flip(sig,axis=0)
     # sig = np.flip(sig,axis=1)
     res=np.eye(4)
@@ -402,12 +403,16 @@ def XCAT_to_MR_DCE(XCAT, TR, TE, bvalue, D, f, Ds, b0=3, ivim_cont = True):
     nifti_img = nib.Nifti1Image(sig, affine=res)  # Replace affine if necessary
     # Save the NIfTI image to a file
     nifti_img.header.set_data_dtype(np.float64)
-    if motion:
+    if not motion:
+        output_file = 'output.nii.gz'  # Replace with your desired output file name
+    elif interleaved:
         output_file = 'output_resp_int.nii.gz'  # Replace with your desired output file name
     else:
-        output_file = 'output.nii.gz'  # Replace with your desired output file name
+        output_file = 'output_resp.nii.gz'  # Replace with your desired output file name
+
     nib.save(nifti_img, output_file)
 
+
     nifti_img = nib.Nifti1Image(XCAT, affine=res)  # Replace affine if necessary
     # Save the NIfTI image to a file
     output_file = 'output_xcat.nii.gz'  # Replace with your desired output file name

src/original/OGC_AmsterdamUMC/LSQ_fitting.py

@@ -139,7 +139,7 @@ def fit_segmented(bvalues, dw_data, bounds=([0, 0, 0.005],[0.005, 0.7, 0.2]), cu
         Dt, Fp = params[0] / 1000, 1 - params[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], bounds[1][2])
+        bounds2 = (bounds[0][2]*10, bounds[1][2]*10)
         # 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]
@@ -561,19 +561,19 @@ def neg_log_prior(p):
             Dt, Fp, Dp = p[0], p[1], p[2]
         # make D*<D very unlikely
         if (Dp < Dt):
-            return 1e8
+            return 1e3
         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
                     return 0
                 else:
-                    return 1e8
+                    return 1e3
             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 1e8
+                    return 1e3
 
     return neg_log_prior
 
@@ -638,7 +638,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], fitS0=True):
+def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05, 1], fitS0=True):
     '''
     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 +655,7 @@ def fit_bayesian(bvalues, dw_data, neg_log_prior, x0=[0.001, 0.2, 0.05], fitS0=T
     '''
     try:
         # define fit bounds
-        bounds = [(0, 0.005), (0, 1), (0, 2), (0, 2.5)]
+        bounds = [(0, 0.005), (0, 1.5), (0, 2), (0, 2.5)]
         # 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/OGC_AmsterdamUMC_Bayesian_biexp.py

@@ -1,5 +1,5 @@
 from src.wrappers.OsipiBase import OsipiBase
-from src.original.OGC_AmsterdamUMC.LSQ_fitting import flat_neg_log_prior, fit_bayesian
+from src.original.OGC_AmsterdamUMC.LSQ_fitting import flat_neg_log_prior, fit_bayesian, empirical_neg_log_prior, fit_segmented
 import numpy as np
 
 class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
@@ -27,23 +27,32 @@ class OGC_AmsterdamUMC_Bayesian_biexp(OsipiBase):
     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):
+    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, prior_in=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.
+
+            Args:
+                 datain is a 2D array with values of D, f, D* (and S0) that will form the prior.
         """
-        super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues, bounds,initial_guess,fitS0)
+        super(OGC_AmsterdamUMC_Bayesian_biexp, self).__init__(bvalues, bounds,initial_guess,fitS0,prior_in)
         if bounds is None:
-            self.bounds=([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.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][2],self.bounds[1][2]],[self.bounds[0][3],self.bounds[1][3]])
+        if prior_in is None:
+            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:
+            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])
         self.OGC_algorithm = fit_bayesian
         if initial_guess is None:
-            self.initial_guess = [0.001, 0.001, 0.01, 1]
+            self.initial_guess = [0.001, 0.5, 0.1, 1]
         else:
             self.initial_guess = initial_guess
         self.fitS0=fitS0
@@ -59,7 +68,9 @@ def ivim_fit(self, signals, bvalues=None):
             _type_: _description_
         """
         bvalues=np.array(bvalues)
-        fit_results = self.OGC_algorithm(bvalues, signals, self.neg_log_prior, x0=self.initial_guess, fitS0=self.fitS0)
+        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)
 
         D = fit_results[0]
         f = fit_results[1]

src/standardized/OGC_AmsterdamUMC_biexp_segmented.py

@@ -42,7 +42,7 @@ def __init__(self, bvalues=None, thresholds=150, bounds=([0, 0, 0.005],[0.005, 0
         else:
             self.bounds=bounds
         if initial_guess is None:
-            self.initial_guess = [0.001, 0.001, 0.01, 1]
+            self.initial_guess = [0.001, 0.1, 0.03, 1]
         else:
             self.initial_guess = initial_guess
         self.thresholds=thresholds