Author(s): Paul Miles, Jason Hite | Date Created: August 15, 2018
This example demonstrates how MCMC methods can be used to locate unknown sources of gamma radiation in an urban environment. For this demonstration we consider a simulated 250m x 178m block of downtown Washington D.C.
We utilize a highly simplified radiation transport model which ignores scattering. The model consists of the following components. There are $N_d$ gamma detectors located at positions $\mathbf{r}_i, i = 1, ..., N_d$. A ray tracing algorithm is implemented for a source located at $\mathbf{r}_0 = (x_0, y_0)$ with intensity $I_0$. Buildings affect the readings of the gamma detectors, so we have modeled their impact using a exponential decay function that depends on the building dimensions and average cross-section. For the $i^{th}$ detector, we let $N_i$ denote the number of buildings between $\mathbf{r}_0$ and $\mathbf{r}_i$. For $n_i = 1,...,N_i$, we let $s_{n_i}$ and $\sigma_{n_i}$ denote the length of the building segment and the average cross-section of the $n_i^{th}$ building, respectively. The surface area, efficiency, and measurement duration for the $i^{th}$ detector are denoted by $A_i, \varepsilon_i,$ and $\Delta t_i$, respectively. Finally, $B_i$ represents the expected background count rate at position $\mathbf{r}_i$. $$ \Gamma_i = \frac{\Delta t_i\varepsilon_iA_iI_0}{4\pi|\mathbf{r}_i - \mathbf{r}_0|^2}\exp\Big(-\sum_{n_i=1}^{N_i}\sigma_{n_i}s_{n_i}\Big) + E[B_i\Delta t_i] $$
This model is implemented in the open source Python package gefry3. The input deck for this example can be found here. For additional reading on this area of research, please refer to the following articles:
Acknowledgment: Funding for this research provided by the Consortium for Nonproliferation Enabling Capabilities (CNEC).
Package Installation: All packages can be install via
pip install package_name
with the exception of gefry3, which requires
pip install git+https://github.com/jasonmhite/gefry3
Advanced statistical plots are generated using the mcmcplot package.
# import required packages
import gefry3
import numpy as np
from pymcmcstat.MCMC import MCMC
import mcmcplot.mcmatplot as mcmpl
import mcmcplot.mcseaborn as mcsbn
from descartes.patch import PolygonPatch
import seaborn as sns
import matplotlib.pyplot as plt
np.seterr(over = 'ignore')
NOBS = 1 # number of observations for each detector
NS = int(1e3) # number of MCMC simulations
# Setup Radiation Model - Downtown Urban Environment Simulation
P = gefry3.read_input_problem('g3_deck.yml', 'Simple_Problem')
# Setup background
nominal_background = 300
dwell = np.array([detector.dwell for detector in P.detectors])
B0 = nominal_background * dwell
# True source and intensity
S0 = P.source.R # m
I0 = P.source.I0 # Bq
nominal = P(S0, I0)
ndet = len(P.detectors)
# Define observations
observations = np.random.poisson(
nominal + nominal_background * dwell,
size=(
NOBS,
ndet,
),
)
observations = observations.astype(np.float64)
# Setup parameter bounds
XMIN, YMIN, XMAX, YMAX = P.domain.all.bounds
IMIN, IMAX = I0 * 0.1, I0 * 10
XMIN, YMIN, XMAX, YMAX = P.domain.all.bounds
print(observations)
print('S0 = {} m, I0 = {} Bq'.format(S0, I0))
[[1524. 1602. 1506. 1678. 1734. 1576. 1523. 1850. 1773. 1707.]] S0 = [158. 98.] m, I0 = 3214000000.0 Bq
# Radiation Sum of Squares Function
def radiation_ssfun(theta, data):
x, y, I = theta
model, background = data.user_defined_object[0]
output = model((x, y), I) + background
res = data.ydata[0] - output
ss = (res ** 2).sum(axis = 0)
return ss
Note, we pass the radiation model into our residual function via the data structure. Also note, the background is included as a parameter, but it is not sampled during the simulation.
# Initialize MCMC object
mcstat = MCMC()
# setup data structure for dram
mcstat.data.add_data_set(
x=np.zeros(observations.shape),
y=observations,
user_defined_object=[
P,
B0,
],
)
mcstat.parameters.add_model_parameter(
name='$x_0$',
theta0=140.,
minimum=XMIN,
maximum=XMAX,
)
mcstat.parameters.add_model_parameter(
name='$y_0$',
theta0=85.,
minimum=YMIN,
maximum=YMAX,
)
mcstat.parameters.add_model_parameter(
name='$I_0$',
theta0=3.0e9,
minimum=IMIN,
maximum=IMAX
)
mcstat.simulation_options.define_simulation_options(
nsimu=NS,
updatesigma=False,
method='dram',
adaptint=100,
verbosity=1,
waitbar=1,
save_to_json=False,
save_to_bin=False,
)
mcstat.model_settings.define_model_settings(
sos_function=radiation_ssfun,
sigma2=observations.sum(axis=0),
)
# Run mcmcrun
mcstat.run_simulation()
Sampling these parameters: name start [ min, max] N( mu, sigma^2) $x_0$: 140.00 [ 0.00e+00, 246.62] N( 0.00e+00, inf) $y_0$: 85.00 [ 0.00e+00, 176.33] N( 0.00e+00, inf) $I_0$: 3.00e+09 [ 3.21e+08, 3.21e+10] N( 0.00e+00, inf) [-----------------100%-----------------] 1000 of 1000 complete in 75.2 sec
results = mcstat.simulation_results.results
# specify burnin period
burnin = int(results['nsimu']/2)
# display chain statistics
chain = results['chain']
s2chain = results['s2chain']
sschain = results['sschain']
names = results['names']
mcstat.chainstats(chain[burnin:,:], results)
print('Acceptance rate: {:6.4}%'.format(100*(1 - results['total_rejected'])))
print('Model Evaluations: {}'.format(results['nsimu']
- results['iacce'][0]
+ results['nsimu']))
# generate mcmc plots
sns.set_style('white')
mcmpl.plot_density_panel(chain[burnin:,:], names)
mcmpl.plot_chain_panel(chain[burnin:,:], names)
mcmpl.plot_pairwise_correlation_panel(chain[burnin:, :], names)
--------------------- name : mean std MC_err tau geweke $x_0$ : 155.4023 1.3916 0.1796 10.1524 0.9958 $y_0$ : 100.4459 1.5950 0.2295 13.2855 0.9903 $I_0$ : 3.74e+09 2.983e+08 45261055.0304 17.3932 0.9690 --------------------- Acceptance rate: 70.7% Model Evaluations: 1833
(<Figure size 700x500 with 3 Axes>, {'skip': 1, 'fig': {'figsize': (7, 5), 'dpi': 100}, 'plot': {'color': 'b', 'marker': '.', 'linestyle': 'none'}, 'xlabel': {}, 'ylabel': {}, 'title': {}, 'add_5095_contours': False, 'plot_50': {'color': 'g', 'marker': None, 'linewidth': 2, 'linestyle': '--', 'label': '50%'}, 'plot_95': {'color': 'r', 'marker': None, 'linewidth': 2, 'linestyle': '--', 'label': '95%'}, 'add_legend': False, 'legend': {'loc': 1}})
First we define a function to allow us to add axes within a figure.
def add_subplot_axes(ax,rect,**kwargs):
fig = plt.gcf()
box = ax.get_position()
width = box.width
height = box.height
inax_position = ax.transAxes.transform(rect[0:2])
transFigure = fig.transFigure.inverted()
infig_position = transFigure.transform(inax_position)
x = infig_position[0]
y = infig_position[1]
width *= rect[2]
height *= rect[3] # <= Typo was here
subax = fig.add_axes([x,y,width,height],**kwargs)
x_labelsize = subax.get_xticklabels()[0].get_size()
y_labelsize = subax.get_yticklabels()[0].get_size()
x_labelsize *= rect[2]**0.5
y_labelsize *= rect[3]**0.5
subax.xaxis.set_tick_params(labelsize=x_labelsize)
subax.yaxis.set_tick_params(labelsize=y_labelsize)
return subax
We set the context of our displays to "talk" which will improve their visibility for presentation. To begin, we will plot the joint distribution relative to the true source location.
sns.set_context('talk')
settings = dict(jointplot=dict(kind='kde', color='b'),
marginal_kws=dict(alpha=1))
fig0, us0 = mcsbn.plot_joint_distributions(
chain[burnin:,0:2],
names,
sns_style='white',
settings=settings
);
# add source
fig0[0].ax_joint.plot(P.source.R[0], P.source.R[1], 'or',
linewidth=3, markerfacecolor='r',
markersize=10, label='True Source');
Next, we plot the joint distributions relative to our entire simulated urban environment. We have added an insert to highlight the relative positioning of the posterior distributions and the true source location.
settings = dict(jointplot=dict(kind='kde', color='b'),
marginal_kws=dict(alpha=1))
fig, us = mcsbn.plot_joint_distributions(chain[burnin:, 0:2],
names,
sns_style='white',
settings=settings);
f = fig[0]
f.ax_marg_x.set_xlim([XMIN, XMAX])
f.ax_marg_y.set_ylim([YMIN, YMAX])
ax = f.ax_joint
# add source
ax.plot(P.source.R[0], P.source.R[1], 'or',
linewidth=3, markerfacecolor='r', markersize=3)
# add buildings to plot
for building in P.domain.solids:
patch = PolygonPatch(building.geom, alpha=0.3,
facecolor='g', edgecolor='k')
ax.add_patch(patch)
# Create insert
rect = [0.1, 0.1, 0.3, 0.3]
ax1 = add_subplot_axes(ax, rect, facecolor='w')
sns.kdeplot(chain[:, 0], chain[:, 1], ax=ax1, color='b')
mu = results['mean']
xlimits = [mu[0] - 7.5, mu[0] + 7.5]
ylimits = [mu[1] - 7.5, mu[1] + 7.5]
ax1.set_xlim(xlimits)
ax1.set_ylim(ylimits)
# add source
ax1.plot(P.source.R[0], P.source.R[1], 'or',
linewidth=3, markerfacecolor='r', markersize=10)
# add buildings to plot
for building in P.domain.solids:
patch = PolygonPatch(building.geom, alpha=0.3,
facecolor='g', edgecolor='k')
ax1.add_patch(patch)
# Create lines outlining insert location
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
mark_inset(ax, ax1, loc1=2, loc2=4,
facecolor="none", edgecolor="k", linewidth=2)
ax1.xaxis.set_visible(False)
ax1.yaxis.set_visible(False)
f.savefig('source_distribution.png')
We observe that the poster distributions for the source location are within a couple meters of the true source location (red circle - $R_{true} = [158, 98]$ meters). Furthermore, the true location is within the uncertainty bounds of our posteriors.