{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Radiation Source Localization\n", "Author(s): Paul Miles, Jason Hite | Date Created: August 15, 2018\n", "\n", "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.\n", "\n", "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$. \n", "$$\n", "\\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]\n", "$$\n", "\n", "This model is implemented in the open source Python package [gefry3](https://github.com/jasonmhite/gefry3). The input deck for this example can be found [here](https://github.com/jasonmhite/gefry3/blob/master/examples/g3_deck.yml). For additional reading on this area of research, please refer to the following articles:\n", "\n", "- Hite, J., & Mattingly, J. (2018). Bayesian Metropolis Methods for Source Localization in an Urban Environment. Radiation Physics and Chemistry. [https://doi.org/10.1016/j.radphyschem.2018.06.024](https://doi.org/10.1016/j.radphyschem.2018.06.024)\n", "- Ştefănescu, R., Schmidt, K., Hite, J., Smith, R. C., & Mattingly, J. (2017). Hybrid optimization and Bayesian inference techniques for a non‐smooth radiation detection problem. International Journal for Numerical Methods in Engineering, 111(10), 955-982. [https://doi.org/10.1002/nme.5491](https://doi.org/10.1002/nme.5491)\n", "- Hite, J. M., Mattingly, J. K., Schmidt, K. L., Ştefănescu, R., & Smith, R. (2016, September). Bayesian metropolis methods applied to sensor networks for radiation source localization. In Multisensor Fusion and Integration for Intelligent Systems (MFI), 2016 IEEE International Conference on (pp. 389-393). IEEE. [https://doi.org/10.1109/MFI.2016.7849519](https://doi.org/10.1109/MFI.2016.7849519)\n", "\n", "__Acknowledgment:__\n", "Funding for this research provided by the [Consortium for Nonproliferation Enabling Capabilities (CNEC)](https://cnec.ncsu.edu/).\n", "\n", "__Package Installation:__\n", "All packages can be install via\n", "```python\n", "pip install package_name\n", "```\n", "with the exception of [gefry3](https://github.com/jasonmhite/gefry3), which requires\n", "```python\n", "pip install git+https://github.com/jasonmhite/gefry3\n", "```\n", "Advanced statistical plots are generated using the [mcmcplot](https://prmiles.wordpress.ncsu.edu/codes/python-packages/mcmcplot/) package." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# import required packages\n", "import gefry3\n", "import numpy as np\n", "from pymcmcstat.MCMC import MCMC\n", "import mcmcplot.mcmatplot as mcmpl\n", "import mcmcplot.mcseaborn as mcsbn\n", "from descartes.patch import PolygonPatch\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "np.seterr(over = 'ignore')\n", "NOBS = 1 # number of observations for each detector\n", "NS = int(1e3) # number of MCMC simulations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Radiation Model, Problem Domain, and Background Radiation" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1524. 1602. 1506. 1678. 1734. 1576. 1523. 1850. 1773. 1707.]]\n", "S0 = [158. 98.] m, I0 = 3214000000.0 Bq\n" ] } ], "source": [ "# Setup Radiation Model - Downtown Urban Environment Simulation\n", "P = gefry3.read_input_problem('g3_deck.yml', 'Simple_Problem')\n", "# Setup background\n", "nominal_background = 300\n", "dwell = np.array([detector.dwell for detector in P.detectors])\n", "B0 = nominal_background * dwell\n", "# True source and intensity\n", "S0 = P.source.R # m\n", "I0 = P.source.I0 # Bq\n", "nominal = P(S0, I0)\n", "ndet = len(P.detectors)\n", "# Define observations\n", "observations = np.random.poisson(\n", " nominal + nominal_background * dwell,\n", " size=(\n", " NOBS,\n", " ndet,\n", " ),\n", ")\n", "observations = observations.astype(np.float64)\n", "# Setup parameter bounds\n", "XMIN, YMIN, XMAX, YMAX = P.domain.all.bounds\n", "IMIN, IMAX = I0 * 0.1, I0 * 10\n", "XMIN, YMIN, XMAX, YMAX = P.domain.all.bounds\n", "print(observations)\n", "print('S0 = {} m, I0 = {} Bq'.format(S0, I0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define Sum-of-Squares Error Function." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Radiation Sum of Squares Function\n", "def radiation_ssfun(theta, data):\n", " x, y, I = theta\n", " model, background = data.user_defined_object[0]\n", " output = model((x, y), I) + background\n", "\n", " res = data.ydata[0] - output\n", " ss = (res ** 2).sum(axis = 0)\n", " return ss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Initialize MCMC Object and Setup Simulation\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Initialize MCMC object\n", "mcstat = MCMC()\n", "# setup data structure for dram\n", "mcstat.data.add_data_set(\n", " x=np.zeros(observations.shape),\n", " y=observations,\n", " user_defined_object=[\n", " P,\n", " B0,\n", " ],\n", ")\n", "mcstat.parameters.add_model_parameter(\n", " name='$x_0$',\n", " theta0=140.,\n", " minimum=XMIN,\n", " maximum=XMAX,\n", ")\n", "mcstat.parameters.add_model_parameter(\n", " name='$y_0$',\n", " theta0=85.,\n", " minimum=YMIN,\n", " maximum=YMAX,\n", ")\n", "mcstat.parameters.add_model_parameter(\n", " name='$I_0$',\n", " theta0=3.0e9,\n", " minimum=IMIN,\n", " maximum=IMAX\n", ")\n", "\n", "mcstat.simulation_options.define_simulation_options(\n", " nsimu=NS,\n", " updatesigma=False,\n", " method='dram',\n", " adaptint=100,\n", " verbosity=1,\n", " waitbar=1,\n", " save_to_json=False,\n", " save_to_bin=False,\n", ")\n", "\n", "mcstat.model_settings.define_model_settings(\n", " sos_function=radiation_ssfun,\n", " sigma2=observations.sum(axis=0),\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Run Simulation" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Sampling these parameters:\n", " name start [ min, max] N( mu, sigma^2)\n", " $x_0$: 140.00 [ 0.00e+00, 246.62] N( 0.00e+00, inf)\n", " $y_0$: 85.00 [ 0.00e+00, 176.33] N( 0.00e+00, inf)\n", " $I_0$: 3.00e+09 [ 3.21e+08, 3.21e+10] N( 0.00e+00, inf)\n", " [-----------------100%-----------------] 1000 of 1000 complete in 75.2 sec" ] } ], "source": [ "# Run mcmcrun\n", "mcstat.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Process Results and Display Statistics" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "---------------------\n", "name : mean std MC_err tau geweke\n", "$x_0$ : 155.4023 1.3916 0.1796 10.1524 0.9958\n", "$y_0$ : 100.4459 1.5950 0.2295 13.2855 0.9903\n", "$I_0$ : 3.74e+09 2.983e+08 45261055.0304 17.3932 0.9690\n", "---------------------\n", "Acceptance rate: 70.7%\n", "Model Evaluations: 1833\n" ] }, { "data": { "text/plain": [ "(