Skip to content

Commit 68a9d33

Browse files
etpetersonetpeterson
committed
Initial sampling commit
1 parent 5c217fc commit 68a9d33

File tree

2 files changed

+183
-5
lines changed

2 files changed

+183
-5
lines changed

‎src/original/ETP_SRI/Sampling.py‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
import numpy as np
2+
import emcee
3+
import tqdm
4+
import corner
5+
from scipy.stats import rice, norm, uniform
6+
from matplotlib import pyplot as plt
7+
8+
from utilities.data_simulation.GenerateData import GenerateData
9+
10+
class MCMC:
11+
"""
12+
Performs sampling of exponential data
13+
"""
14+
def __init__(self,
15+
data,
16+
b,
17+
data_scale=1e-5,
18+
parameter_scale=(1e-7, 1e-11, 1e-9),
19+
bounds=((0, 1), (0, 1), (0, 1)),
20+
priors=None,
21+
gaussian_noise=False,
22+
nwalkers=16,
23+
nsteps=10000,
24+
burn_in=2000,
25+
progress=True):
26+
"""
27+
Parameters
28+
----------
29+
30+
"""
31+
self.data = np.atleast_2d(np.asarray(data))
32+
self.b = np.atleast_2d(np.asarray(b)).T
33+
self.data_scale = np.asarray(data_scale)
34+
self.parameter_scale = np.asarray(parameter_scale)
35+
self.bounds = np.asarray(bounds)
36+
self.priors = np.asarray(priors)
37+
self.nwalkers = nwalkers
38+
self.nsteps = nsteps
39+
self.burn_in = burn_in
40+
self.progress = progress
41+
if priors is None:
42+
self.prior = self.zero_prior
43+
else:
44+
self.prior = self.loglike_gauss_prior
45+
if gaussian_noise:
46+
self.likelihood = self.biexp_loglike_gauss
47+
else:
48+
self.likelihood = self.biexp_loglike_rice
49+
self.ndim = 3
50+
self.chain = None
51+
self.means = None
52+
self.stds = None
53+
54+
def accepted_dimensions(self):
55+
return (1, 1)
56+
57+
def bounds_prior(self, params):
58+
return np.sum(uniform.logpdf(params, loc=self.bounds[:, 0], scale=self.bounds[:, 1] - self.bounds[:, 0]), 1)
59+
60+
def loglike_gauss_prior(self, params):
61+
return np.sum(norm.logpdf(params, loc=self.priors[:, 0], scale=self.priors[:, 1]), 1)
62+
63+
def zero_prior(self, params):
64+
return 0
65+
66+
def signal(self, f, D, D_star):
67+
return (f * np.exp(-self.b * D_star) + (1 - f) * np.exp(-self.b * D)).T
68+
69+
def biexp_loglike_gauss(self, f, D, D_star):
70+
expected = self.signal(f, D, D_star)
71+
# check this!
72+
# print(f'likelihood {norm.logpdf(self.data, loc=expected/self.data_scale, scale=self.data_scale)}')
73+
return np.sum(norm.logpdf(self.data, loc=expected, scale=self.data_scale), 1)
74+
75+
# def biexp_loglike_gauss_full(self, f, D, D_star):
76+
# expected = self.signal(f, D, D_star)
77+
# print(f'expected {expected}')
78+
# print(f'data {self.data}')
79+
# return norm.logpdf(self.data, loc=expected, scale=self.data_scale)
80+
81+
def biexp_loglike_rice(self, f, D, D_star):
82+
expected = self.signal(f, D, D_star)
83+
# print(f'expected {expected}')
84+
return np.sum(rice.logpdf(self.data, b=expected/self.data_scale, scale=self.data_scale), 1)
85+
86+
def posterior(self, params):
87+
params = np.atleast_2d(params)
88+
total = self.bounds_prior(params)
89+
# print(f'bounds params {total}')
90+
neginf = np.isneginf(total)
91+
# print(f'neginf {neginf}')
92+
f = params[~neginf, 0]
93+
D = params[~neginf, 1]
94+
D_star = params[~neginf, 2]
95+
prior = self.prior(params[~neginf, :])
96+
# print(f'prior {prior}')
97+
likelihood = self.likelihood(f, D, D_star)
98+
# print(f'likelihood {likelihood}')
99+
total[~neginf] += prior + likelihood
100+
return total
101+
102+
def sample(self, initial_pos):
103+
# f = initial_pos[0]
104+
# D = initial_pos[1]
105+
# D_star = initial_pos[2]
106+
# print(f'initial pos likelihood {self.biexp_loglike_gauss_full(f, D, D_star)}')
107+
print(f'initial pos likelihood {self.posterior(initial_pos)}')
108+
sampler = emcee.EnsembleSampler(self.nwalkers, 3, self.posterior, vectorize=True)
109+
pos = initial_pos + self.parameter_scale * np.random.randn(self.nwalkers, self.ndim)
110+
# print(f'pos {pos}')
111+
# print(f'nsteps {self.nsteps}')
112+
sampler.run_mcmc(pos, self.nsteps, progress=True)
113+
self.chain = sampler.get_chain(discard=self.burn_in, flat=True)
114+
self.means = np.mean(self.chain, 0)
115+
self.stds = np.std(self.chain, 0)
116+
print(f'final pos likelihood {self.posterior(self.means)}')
117+
# print(f'final pos likelihood {self.biexp_loglike_gauss_full(self.means[0], self.means[1], self.means[2])}')
118+
# print(f'chain {self.chain}')
119+
return self.means, self.stds
120+
121+
def plot(self, truths=None, labels=('f', 'D', 'D*'), overplot=None):
122+
if truths is None:
123+
truths = self.means
124+
# print(f'chain size {self.chain.shape}')
125+
fig = corner.corner(self.chain, labels=labels, truths=truths)
126+
fig.suptitle("Sampling of the IVIM data", fontsize=16)
127+
if overplot is not None:
128+
corner.overplot_lines(fig, overplot, color='r')
129+
plt.show()
130+

‎tests/IVIMmodels/unit_tests/simple_test_run_of_algorithm.py‎

Lines changed: 53 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,33 +1,81 @@
11
import numpy as np
2+
import matplotlib.pyplot as plt
23
import os
34
from pathlib import Path
45
#from src.standardized.ETP_SRI_LinearFitting import ETP_SRI_LinearFitting
5-
from src.standardized.IAR_LU_biexp import IAR_LU_biexp
6+
# from src.standardized.IAR_LU_biexp import IAR_LU_biexp
67
#from src.standardized.IAR_LU_segmented_2step import IAR_LU_segmented_2step
78
from src.standardized.PvH_KB_NKI_IVIMfit import PvH_KB_NKI_IVIMfit
89
#from src.standardized.PV_MUMC_biexp import PV_MUMC_biexp
10+
from src.original.ETP_SRI.Sampling import MCMC
911

1012
## Simple test code...
1113
# Used to just do a test run of an algorithm during development
1214
def dev_test_run(model, **kwargs):
13-
bvalues = np.array([0, 50, 100, 150, 200, 500, 800])
15+
bvalues = np.array([0, 20, 50, 75, 100, 150, 200, 300, 400, 500, 800, 1000, 1500])
1416

15-
def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
17+
def ivim_model(b, f=0.1, Dstar=0.01, D=0.001, S0=1):
18+
# print(f'S0 {S0}')
19+
# print(f'Dstar {f*np.exp(-b*Dstar)}')
20+
# print(f'D {(1-f)*np.exp(-b*D)}')
21+
# print(f'sum {f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D)}')
22+
# print(f'S0 {(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))}')
1623
return S0*(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))
1724

18-
signals = ivim_model(bvalues)
25+
# TODO: add S0 fitting!
26+
true_f = 0.4
27+
true_Dstar = 0.01
28+
true_D = 0.001
29+
truth = [true_f, true_D, true_Dstar]
30+
signals_noiseless = ivim_model(bvalues, true_f, true_Dstar, true_D)
31+
print(f'noiselss {signals_noiseless}')
32+
signals = signals_noiseless + np.abs(1e-1 * (np.random.randn(len(bvalues)) + 1j * np.random.randn(len(bvalues))) / np.sqrt(2))
1933

2034
#model = ETP_SRI_LinearFitting(thresholds=[200])
2135
if kwargs:
2236
results = model.osipi_fit(signals, bvalues, **kwargs)
2337
else:
2438
results = model.osipi_fit(signals, bvalues)
25-
print(results)
39+
# print(results) # f, D*, D
40+
results_reordered = np.asarray([results[0], results[2], results[1]])
41+
print(truth)
42+
print(results_reordered)
2643
#test = model.osipi_simple_bias_and_RMSE_test(SNR=20, bvalues=bvalues, f=0.1, Dstar=0.03, D=0.001, noise_realizations=10)
44+
signal_results = ivim_model(bvalues, results[0], results[1], results[2])
45+
46+
47+
48+
mcmc = MCMC(signals, bvalues, gaussian_noise=False, data_scale=1e-2) #, priors=((0.07, 1e-1), (0.0135, 1e-1), (0.001, 1e-1)))
49+
means, stds = mcmc.sample(truth)
50+
print(f'means {means} stds {stds}')
51+
print(f'expected {results_reordered}')
52+
print(f'truth {truth}')
53+
54+
signal_means = ivim_model(bvalues, means[0], means[2], means[1])
55+
plt.plot(bvalues, signals_noiseless, 'g', label='Noiseless Signal')
56+
plt.plot(bvalues, signals, '.g', label='Noisy Signal')
57+
plt.plot(bvalues, signal_results, 'r', label='Results Signal')
58+
plt.plot(bvalues, signal_means, 'b', label='Means Signal')
59+
plt.legend()
60+
# plt.show()
61+
62+
mcmc.plot(overplot=truth)
63+
2764

2865
#model1 = ETP_SRI_LinearFitting(thresholds=[200])
2966
#model2 = IAR_LU_biexp()
3067
model2 = PvH_KB_NKI_IVIMfit()
3168

3269
#dev_test_run(model1, linear_fit_option=True)
3370
dev_test_run(model2)
71+
72+
73+
def run_sampling():
74+
bvalues = np.array([0, 50, 100, 150, 200, 500, 800])
75+
76+
def ivim_model(b, S0=1, f=0.1, Dstar=0.01, D=0.001):
77+
return S0*(f*np.exp(-b*Dstar) + (1-f)*np.exp(-b*D))
78+
79+
signals = ivim_model(bvalues)
80+
81+

0 commit comments

Comments
 (0)