From 68a9d337aae30fc91fbb18eb4fa2fec1d7840937 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Tue, 26 Nov 2024 21:36:05 -0800 Subject: [PATCH 1/4] Initial sampling commit --- src/original/ETP_SRI/Sampling.py | 130 ++++++++++++++++++ .../simple_test_run_of_algorithm.py | 58 +++++++- 2 files changed, 183 insertions(+), 5 deletions(-) create mode 100644 src/original/ETP_SRI/Sampling.py diff --git a/src/original/ETP_SRI/Sampling.py b/src/original/ETP_SRI/Sampling.py new file mode 100644 index 00000000..5de8b991 --- /dev/null +++ b/src/original/ETP_SRI/Sampling.py @@ -0,0 +1,130 @@ +import numpy as np +import emcee +import tqdm +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 + """ + 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 + ---------- + + """ + 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) + # check this! + # print(f'likelihood {norm.logpdf(self.data, loc=expected/self.data_scale, scale=self.data_scale)}') + return np.sum(norm.logpdf(self.data, loc=expected, scale=self.data_scale), 1) + + # def biexp_loglike_gauss_full(self, f, D, D_star): + # expected = self.signal(f, D, D_star) + # print(f'expected {expected}') + # print(f'data {self.data}') + # return norm.logpdf(self.data, loc=expected, scale=self.data_scale) + + def biexp_loglike_rice(self, f, D, D_star): + expected = self.signal(f, D, D_star) + # print(f'expected {expected}') + 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) + # print(f'bounds params {total}') + neginf = np.isneginf(total) + # print(f'neginf {neginf}') + f = params[~neginf, 0] + D = params[~neginf, 1] + D_star = params[~neginf, 2] + prior = self.prior(params[~neginf, :]) + # print(f'prior {prior}') + likelihood = self.likelihood(f, D, D_star) + # print(f'likelihood {likelihood}') + total[~neginf] += prior + likelihood + return total + + def sample(self, initial_pos): + # f = initial_pos[0] + # D = initial_pos[1] + # D_star = initial_pos[2] + # print(f'initial pos likelihood {self.biexp_loglike_gauss_full(f, D, D_star)}') + 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) + # print(f'pos {pos}') + # print(f'nsteps {self.nsteps}') + 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)}') + # print(f'final pos likelihood {self.biexp_loglike_gauss_full(self.means[0], self.means[1], self.means[2])}') + # print(f'chain {self.chain}') + return self.means, self.stds + + def plot(self, truths=None, labels=('f', 'D', 'D*'), overplot=None): + if truths is None: + truths = self.means + # print(f'chain size {self.chain.shape}') + fig = corner.corner(self.chain, labels=labels, truths=truths) + fig.suptitle("Sampling of the IVIM data", fontsize=16) + if overplot is not None: + corner.overplot_lines(fig, overplot, color='r') + plt.show() + \ No newline at end of file diff --git a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py index 92473bbf..90960729 100644 --- a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py +++ b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py @@ -1,29 +1,66 @@ import numpy as np +import matplotlib.pyplot as plt import os from pathlib import Path #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_biexp import IAR_LU_biexp #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 ## 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) + # 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)) #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[0], results[2], results[1]]) + 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[0], results[1], results[2]) + + + + mcmc = MCMC(signals, bvalues, gaussian_noise=False, data_scale=1e-2) #, priors=((0.07, 1e-1), (0.0135, 1e-1), (0.001, 1e-1))) + means, stds = mcmc.sample(truth) + 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) + #model1 = ETP_SRI_LinearFitting(thresholds=[200]) #model2 = IAR_LU_biexp() @@ -31,3 +68,14 @@ def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001): #dev_test_run(model1, linear_fit_option=True) dev_test_run(model2) + + +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) + + From 02f6690423d23ce7710a508f39f1315e240e5f57 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Mon, 17 Mar 2025 23:03:01 -0700 Subject: [PATCH 2/4] Update standardized sampler --- requirements.txt | 2 ++ src/original/ETP_SRI/Sampling.py | 51 ++++++++++++++++---------------- 2 files changed, 27 insertions(+), 26 deletions(-) diff --git a/requirements.txt b/requirements.txt index 29742e33..ad527754 100644 --- a/requirements.txt +++ b/requirements.txt @@ -15,3 +15,5 @@ pandas sphinx sphinx_rtd_theme pytest-json-report +emcee +corner diff --git a/src/original/ETP_SRI/Sampling.py b/src/original/ETP_SRI/Sampling.py index 5de8b991..87e6b3fa 100644 --- a/src/original/ETP_SRI/Sampling.py +++ b/src/original/ETP_SRI/Sampling.py @@ -1,6 +1,5 @@ import numpy as np import emcee -import tqdm import corner from scipy.stats import rice, norm, uniform from matplotlib import pyplot as plt @@ -9,7 +8,7 @@ class MCMC: """ - Performs sampling of exponential data + Performs sampling of exponential data using MCMC """ def __init__(self, data, @@ -26,6 +25,28 @@ def __init__(self, """ 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)) @@ -68,60 +89,38 @@ def signal(self, f, D, D_star): def biexp_loglike_gauss(self, f, D, D_star): expected = self.signal(f, D, D_star) - # check this! - # print(f'likelihood {norm.logpdf(self.data, loc=expected/self.data_scale, scale=self.data_scale)}') return np.sum(norm.logpdf(self.data, loc=expected, scale=self.data_scale), 1) - # def biexp_loglike_gauss_full(self, f, D, D_star): - # expected = self.signal(f, D, D_star) - # print(f'expected {expected}') - # print(f'data {self.data}') - # return norm.logpdf(self.data, loc=expected, scale=self.data_scale) - def biexp_loglike_rice(self, f, D, D_star): expected = self.signal(f, D, D_star) - # print(f'expected {expected}') 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) - # print(f'bounds params {total}') neginf = np.isneginf(total) - # print(f'neginf {neginf}') f = params[~neginf, 0] D = params[~neginf, 1] D_star = params[~neginf, 2] prior = self.prior(params[~neginf, :]) - # print(f'prior {prior}') likelihood = self.likelihood(f, D, D_star) - # print(f'likelihood {likelihood}') total[~neginf] += prior + likelihood return total def sample(self, initial_pos): - # f = initial_pos[0] - # D = initial_pos[1] - # D_star = initial_pos[2] - # print(f'initial pos likelihood {self.biexp_loglike_gauss_full(f, D, D_star)}') - print(f'initial pos likelihood {self.posterior(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) - # print(f'pos {pos}') - # print(f'nsteps {self.nsteps}') 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)}') - # print(f'final pos likelihood {self.biexp_loglike_gauss_full(self.means[0], self.means[1], self.means[2])}') - # print(f'chain {self.chain}') + # 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): if truths is None: truths = self.means - # print(f'chain size {self.chain.shape}') fig = corner.corner(self.chain, labels=labels, truths=truths) fig.suptitle("Sampling of the IVIM data", fontsize=16) if overplot is not None: From 5265d0b0535cfad86c4f946f38cbb682ecf3cc58 Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Mon, 17 Mar 2025 23:04:22 -0700 Subject: [PATCH 3/4] Adding missing standardized wrapper mcmc file --- src/standardized/ETP_SRI_MCMC.py | 109 +++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 src/standardized/ETP_SRI_MCMC.py diff --git a/src/standardized/ETP_SRI_MCMC.py b/src/standardized/ETP_SRI_MCMC.py new file mode 100644 index 00000000..983f7c13 --- /dev/null +++ b/src/standardized/ETP_SRI_MCMC.py @@ -0,0 +1,109 @@ +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, truths=None, labels=('f', 'D', 'D*'), overplot=None): + """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. + """ + if self.MCMC is None: + raise ValueError('No MCMC fit has been performed. Fit the data first.') + self.MCMC.plot(truths, labels, overplot) + From f5e621cf5cd08e81b83d73661265937e8c685f8e Mon Sep 17 00:00:00 2001 From: Eric Peterson Date: Sun, 6 Apr 2025 15:24:49 -0700 Subject: [PATCH 4/4] Fix test timeout error --- src/original/ETP_SRI/Sampling.py | 12 +++++++--- src/standardized/ETP_SRI_MCMC.py | 7 ++++-- .../simple_test_run_of_algorithm.py | 22 +++++++++++-------- 3 files changed, 27 insertions(+), 14 deletions(-) diff --git a/src/original/ETP_SRI/Sampling.py b/src/original/ETP_SRI/Sampling.py index 87e6b3fa..8e3e2ee2 100644 --- a/src/original/ETP_SRI/Sampling.py +++ b/src/original/ETP_SRI/Sampling.py @@ -118,12 +118,18 @@ def sample(self, initial_pos): # 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): + 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("Sampling of the IVIM data", fontsize=16) + fig.suptitle(title, fontsize=16) if overplot is not None: corner.overplot_lines(fig, overplot, color='r') - plt.show() + return fig + # if save_path is not None: + # fig.savefig(save_path) + # if show: + # fig.show() + # if close: + # plt.close(fig) \ No newline at end of file diff --git a/src/standardized/ETP_SRI_MCMC.py b/src/standardized/ETP_SRI_MCMC.py index 983f7c13..a5bfc146 100644 --- a/src/standardized/ETP_SRI_MCMC.py +++ b/src/standardized/ETP_SRI_MCMC.py @@ -95,15 +95,18 @@ def ivim_fit(self, signals, initial_guess=None, bvalues=None, **kwargs): return {'means': means, 'stds': stds} - def plot(self, truths=None, labels=('f', 'D', 'D*'), overplot=None): + 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.') - self.MCMC.plot(truths, labels, overplot) + return self.MCMC.plot(**kwargs) diff --git a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py index 44ff52ed..e211fa66 100644 --- a/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py +++ b/tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py @@ -8,6 +8,8 @@ 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 @@ -19,12 +21,12 @@ #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.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): @@ -90,14 +92,7 @@ def ivim_model(b, f=0.1, Dstar=0.01, D=0.001, S0=1): mcmc.plot(overplot=truth) - -#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) def run_sampling(): @@ -109,3 +104,12 @@ def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001): signals = ivim_model(bvalues) +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)