Sampling/initial

requirements.txt

@@ -15,3 +15,5 @@ pandas
 sphinx
 sphinx_rtd_theme
 pytest-json-report
+emcee
+corner

src/original/ETP_SRI/Sampling.py

@@ -0,0 +1,135 @@
+import numpy as np
+import emcee
+import corner
+from scipy.stats import rice, norm, uniform
+from matplotlib import pyplot as plt
+
+from utilities.data_simulation.GenerateData import GenerateData
+
+class MCMC:
+    """
+    Performs sampling of exponential data using MCMC
+    """
+    def __init__(self,
+                 data,
+                 b,
+                 data_scale=1e-5,
+                 parameter_scale=(1e-7, 1e-11, 1e-9),
+                 bounds=((0, 1), (0, 1), (0, 1)),
+                 priors=None,
+                 gaussian_noise=False,
+                 nwalkers=16,
+                 nsteps=10000,
+                 burn_in=2000,
+                 progress=True):
+        """
+        Parameters
+        ----------
+        data : array_like
+            The data to be fitted
+        b : array_like
+            The b-values for the data
+        data_scale : float, optional
+            The scale of the data, by default 1e-5
+        parameter_scale : tuple, optional
+            The scale of the parameters, by default (1e-7, 1e-11, 1e-9)
+        bounds : tuple, optional
+            The bounds for the parameters, by default ((0, 1), (0, 1), (0, 1))
+        priors : tuple, optional
+            The priors for the parameters, by default None
+        gaussian_noise : bool, optional
+            Whether the noise is gaussian, by default False
+        nwalkers : int, optional
+            The number of walkers, by default 16
+        nsteps : int, optional
+            The number of steps, by default 10000
+        burn_in : int, optional
+            The burn in, by default 2000
+        progress : bool, optional
+            Whether to show progress, by default True
+
+        """
+        self.data = np.atleast_2d(np.asarray(data))
+        self.b = np.atleast_2d(np.asarray(b)).T
+        self.data_scale = np.asarray(data_scale)
+        self.parameter_scale = np.asarray(parameter_scale)
+        self.bounds = np.asarray(bounds)
+        self.priors = np.asarray(priors)
+        self.nwalkers = nwalkers
+        self.nsteps = nsteps
+        self.burn_in = burn_in
+        self.progress = progress
+        if priors is None:
+            self.prior = self.zero_prior
+        else:
+            self.prior = self.loglike_gauss_prior
+        if gaussian_noise:
+            self.likelihood = self.biexp_loglike_gauss
+        else:
+            self.likelihood = self.biexp_loglike_rice
+        self.ndim = 3
+        self.chain = None
+        self.means = None
+        self.stds = None
+
+    def accepted_dimensions(self):
+        return (1, 1)
+
+    def bounds_prior(self, params):
+        return np.sum(uniform.logpdf(params, loc=self.bounds[:, 0], scale=self.bounds[:, 1] - self.bounds[:, 0]), 1)
+
+    def loglike_gauss_prior(self, params):
+        return np.sum(norm.logpdf(params, loc=self.priors[:, 0], scale=self.priors[:, 1]), 1)
+
+    def zero_prior(self, params):
+        return 0
+
+    def signal(self, f, D, D_star):
+        return (f * np.exp(-self.b * D_star) + (1 - f) * np.exp(-self.b * D)).T
+
+    def biexp_loglike_gauss(self, f, D, D_star):
+        expected = self.signal(f, D, D_star)
+        return np.sum(norm.logpdf(self.data, loc=expected, scale=self.data_scale), 1)
+
+    def biexp_loglike_rice(self, f, D, D_star):
+        expected = self.signal(f, D, D_star)
+        return np.sum(rice.logpdf(self.data, b=expected/self.data_scale, scale=self.data_scale), 1)
+
+    def posterior(self, params):
+        params = np.atleast_2d(params)
+        total = self.bounds_prior(params)
+        neginf = np.isneginf(total)
+        f = params[~neginf, 0]
+        D = params[~neginf, 1]
+        D_star = params[~neginf, 2]
+        prior = self.prior(params[~neginf, :])
+        likelihood = self.likelihood(f, D, D_star)
+        total[~neginf] += prior + likelihood
+        return total
+
+    def sample(self, initial_pos):
+        # print(f'initial pos likelihood {self.posterior(initial_pos)}')
+        sampler = emcee.EnsembleSampler(self.nwalkers, 3, self.posterior, vectorize=True)
+        pos = initial_pos + self.parameter_scale * np.random.randn(self.nwalkers, self.ndim)
+        sampler.run_mcmc(pos, self.nsteps, progress=True)
+        self.chain = sampler.get_chain(discard=self.burn_in, flat=True)
+        self.means = np.mean(self.chain, 0)
+        self.stds = np.std(self.chain, 0)
+        # print(f'final pos likelihood {self.posterior(self.means)}')
+        return self.means, self.stds
+
+    def plot(self, truths=None, labels=('f', 'D', 'D*'), overplot=None, title='Sampling of the IVIM data'): #, show=True, save_path=None, close=False):
+        if truths is None:
+            truths = self.means
+        fig = corner.corner(self.chain, labels=labels, truths=truths)
+        fig.suptitle(title, fontsize=16)
+        if overplot is not None:
+            corner.overplot_lines(fig, overplot, color='r')
+        return fig
+        # if save_path is not None:
+        #     fig.savefig(save_path)
+        # if show:
+        #     fig.show()
+        # if close:
+        #     plt.close(fig)
+

src/standardized/ETP_SRI_MCMC.py

@@ -0,0 +1,112 @@
+import numpy as np
+# import emcee
+from src.wrappers.OsipiBase import OsipiBase
+from src.original.ETP_SRI.Sampling import MCMC
+
+
+class ETP_SRI_MCMC(OsipiBase):
+    """WIP
+    Implementation and execution of the submitted algorithm
+    """
+
+    # 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 = "Eric T. Peterson, SRI"
+    id_algorithm_type = "Markov Chain Monte Carlo"
+    id_return_parameters = "Mean, Standard deviation"
+    id_units = "seconds per milli metre squared"
+
+    # Algorithm requirements
+    required_bvalues = 3
+    required_thresholds = [0,1] # Interval from 1 to 1, 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 = 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=None,
+                 data_scale=1e-2,
+                 parameter_scale=(1e-7, 1e-11, 1e-9),
+                 bounds=((0, 1), (0, 1), (0, 1)),
+                 priors=None,
+                 gaussian_noise=False,
+                 nwalkers=16,
+                 nsteps=10000,
+                 burn_in=2000,
+                 progress=True,
+                 initial_guess=None):
+        """
+            Everything this method 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(ETP_SRI_MCMC, self).__init__(bvalues)
+        self.result_keys = ['means', 'stds']
+
+        # Could be a good idea to have all the submission-specfic variable be 
+        # defined with initials?
+        self.bvalues = bvalues
+        self.data_scale = data_scale
+        self.parameter_scale = parameter_scale
+        self.bounds = bounds
+        self.priors = priors
+        self.gaussian_noise = gaussian_noise
+        self.nwalkers = nwalkers
+        self.nsteps = nsteps
+        self.burn_in = burn_in
+        self.progress = progress
+        self.initial_guess = initial_guess
+
+        self.MCMC = None
+
+        # Check the inputs
+
+
+    def ivim_fit(self, signals, initial_guess=None, bvalues=None, **kwargs):
+        """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.
+            linear_fit_option (bool, optional): This fit has an option to only run a linear fit. Defaults to False.
+
+        Returns:
+            _type_: _description_
+        """
+        if bvalues is None:
+            bvalues = self.bvalues
+
+        if initial_guess is not None:
+            self.initial_guess = initial_guess
+
+        self.__dict__.update(kwargs)
+
+        self.MCMC = MCMC(signals, bvalues, self.data_scale, self.parameter_scale, self.bounds, self.priors, self.gaussian_noise, self.nwalkers, self.nsteps, self.burn_in, self.progress)
+
+        means, stds = self.MCMC.sample(self.initial_guess)
+
+        return {'means': means, 'stds': stds}
+
+    def plot(self, **kwargs):
+        """Plot the results of the MCMC fit
+
+        Args:
+            truths (array-like, optional): True values. Defaults to None.
+            labels (tuple, optional): Labels for the parameters. Defaults to ('f', 'D', 'D*').
+            overplot (array-like, optional): Overplot the true values. Defaults to None.
+            title (str, optional): Title of the plot. Defaults to 'MCMC Fit'.
+            show (bool, optional): Show the plot. Defaults to True.
+            save_path (str, optional): Path to save the plot. Defaults to None.
+        """
+        if self.MCMC is None:
+            raise ValueError('No MCMC fit has been performed. Fit the data first.')
+        return self.MCMC.plot(**kwargs)
+

tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py

@@ -1,43 +1,115 @@
+
 import sys
-sys.path.append(r"C:\Users\ivan5\Box\OSIPI\TF2.4_IVIM-MRI_CodeCollection")
+# sys.path.append(r"C:\Users\ivan5\Box\OSIPI\TF2.4_IVIM-MRI_CodeCollection")
+sys.path.append(f".")
 
 
 import os
 import numpy as np
+import matplotlib.pyplot as plt
 from pathlib import Path
+import nibabel as nib
+
 #from src.standardized.ETP_SRI_LinearFitting import ETP_SRI_LinearFitting
-from src.standardized.IAR_LU_biexp import IAR_LU_biexp
-from src.standardized.IAR_LU_modified_mix import IAR_LU_modified_mix
+
+# from src.standardized.IAR_LU_biexp import IAR_LU_biexp
+
+# from src.standardized.IAR_LU_biexp import IAR_LU_biexp
+# from src.standardized.IAR_LU_modified_mix import IAR_LU_modified_mix
+
 #from src.standardized.IAR_LU_segmented_2step import IAR_LU_segmented_2step
 #from src.standardized.PvH_KB_NKI_IVIMfit import PvH_KB_NKI_IVIMfit
 #from src.standardized.PV_MUMC_biexp import PV_MUMC_biexp
+
+# from src.original.ETP_SRI.Sampling import MCMC
+from src.standardized.ETP_SRI_MCMC import ETP_SRI_MCMC
+
 from src.standardized.OGC_AmsterdamUMC_biexp import OGC_AmsterdamUMC_biexp
 
+# /Users/e29007/Stanford/sherlock/sub-3900670/dwi/sub-3900670_dir-AP_desc-denoise_dwi.nii.gz
 ## Simple test code... 
 # Used to just do a test run of an algorithm during development
 def dev_test_run(model, **kwargs):
-    bvalues = np.array([0, 50, 100, 150, 200, 500, 800])
+    bvalues = np.array([0, 20, 50, 75, 100, 150, 200, 300, 400, 500, 800, 1000, 1500])
 
-    def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
+    def ivim_model(b, f=0.1, Dstar=0.01, D=0.001, S0=1):
+        # print(f'S0 {S0}')
+        # print(f'Dstar {f*np.exp(-b*Dstar)}')
+        # print(f'D {(1-f)*np.exp(-b*D)}')
+        # print(f'sum {f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D)}')
+        # print(f'S0 {(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))}')
         return S0*(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))
 
-    signals = ivim_model(bvalues)
-    data = np.array([signals, signals, signals])
-    #print(data)
-    signals = data
+
+    # TODO: add S0 fitting!
+    true_f = 0.4
+    true_Dstar = 0.01
+    true_D = 0.001
+    truth = [true_f, true_D, true_Dstar]
+    signals_noiseless = ivim_model(bvalues, true_f, true_Dstar, true_D)
+    print(f'noiselss {signals_noiseless}')
+    signals = signals_noiseless + np.abs(1e-1 * (np.random.randn(len(bvalues)) + 1j * np.random.randn(len(bvalues))) / np.sqrt(2))
+
+    # signals = ivim_model(bvalues)
+    # data = np.array([signals, signals, signals])
+    # #print(data)
+    # signals = data
+
 
     #model = ETP_SRI_LinearFitting(thresholds=[200])
     if kwargs:
         results = model.osipi_fit(signals, bvalues, **kwargs)
     else:
         results = model.osipi_fit(signals, bvalues)
-    print(results)
+    # print(results) # f, D*, D
+    results_reordered = np.asarray([results['f'], results['Dp'], results['D']])
+    print(truth)
+    print(results_reordered)
     #test = model.osipi_simple_bias_and_RMSE_test(SNR=20, bvalues=bvalues, f=0.1, Dstar=0.03, D=0.001, noise_realizations=10)
+    signal_results = ivim_model(bvalues, results['f'], results['Dp'], results['D'])
+
+
+
+    # mcmc = MCMC(signals, bvalues, gaussian_noise=True, data_scale=1e-2) #, priors=((0.07, 1e-1), (0.0135, 1e-1), (0.001, 1e-1)))
+    # means, stds = mcmc.sample(truth)
+
+    mcmc = ETP_SRI_MCMC(bvalues, gaussian_noise=True, data_scale=1e-2) #, priors=((0.07, 1e-1), (0.0135, 1e-1), (0.001, 1e-1)))
+    mcmc_results = mcmc.ivim_fit(signals, initial_guess=truth)
+    means = mcmc_results['means']
+    stds = mcmc_results['stds']
+
+    print(f'means {means} stds {stds}')
+    print(f'expected {results_reordered}')
+    print(f'truth {truth}')
+
+    signal_means = ivim_model(bvalues, means[0], means[2], means[1])
+    plt.plot(bvalues, signals_noiseless, 'g', label='Noiseless Signal')
+    plt.plot(bvalues, signals, '.g', label='Noisy Signal')
+    plt.plot(bvalues, signal_results, 'r', label='Results Signal')
+    plt.plot(bvalues, signal_means, 'b', label='Means Signal')
+    plt.legend()
+    # plt.show()
+
+    mcmc.plot(overplot=truth)
+
+
+
+
+def run_sampling():
+    bvalues = np.array([0, 50, 100, 150, 200, 500, 800])
+
+    def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
+        return S0*(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))
+
+    signals = ivim_model(bvalues)
+
 
-#model1 = ETP_SRI_LinearFitting(thresholds=[200])
-#model2 = IAR_LU_biexp(bounds=([0,0,0,0], [1,1,1,1]))
-#model2 = IAR_LU_modified_mix()
-model2 = OGC_AmsterdamUMC_biexp()
+if __name__ == "__main__":
+
+    #model1 = ETP_SRI_LinearFitting(thresholds=[200])
+    #model2 = IAR_LU_biexp(bounds=([0,0,0,0], [1,1,1,1]))
+    #model2 = IAR_LU_modified_mix()
+    model2 = OGC_AmsterdamUMC_biexp()
 
-#dev_test_run(model1)
-dev_test_run(model2)
+    #dev_test_run(model1)
+    dev_test_run(model2)