{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU",
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "![MCX Banner](https://mcx.space/img/mcx18_banner.png)\n",
        "\n",
        "# » Getting Started with pmcx\n",
        "\n",
        "- MCX Website: https://mcx.space\n",
        "- License: GNU Public License V3 or later\n",
        "- URL: https://pypi.org/project/pmcx/\n",
        "- Github: https://github.com/fangq/mcx\n",
        "\n",
        "Monte Carlo eXtreme (MCX) is a fast photon transport simulation software for 3D heterogeneous turbid media. By taking advantage of the massively parallel threads and extremely low memory latency in a modern graphics processing unit (GPU), MCX is capable of performing Monte Carlo (MC) photon simulations at a blazing speed, typically hundreds to a thousand times faster than a fully optimized CPU-based MC implementation."
      ],
      "metadata": {
        "id": "bWGJE5zOuMEh"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Step 0: install necessary Python modules\n",
        "\n",
        "> Only need to run this once at the beginning\n",
        "\n",
        "\n",
        "The only required dependency for `pmcx` is `numpy`. Additional (and optional) modules can be installed, including\n",
        "- [`jdata`](https://pypi.org/project/jdata/), for saving simulations to human-readable JSON/[NeuroJSON](https://neurojson.org) files, also for importing MCX .jnii/.jdat/.jdt output files\n",
        "- [`bjdata`](https://pypi.org/project/bjdata/), for saving simulations to more compact binary JSON files, also for importing MCX .bnii/.jdb output files\n",
        "- [`matplotlib`](https://pypi.org/project/matplotlib/), for plotting and visualization of the data\n",
        "\n",
        "If you intend to install pmcx on your local computer (assuming you have an NVIDIA GPU), please first install `python` and run the below `pip` command in a command line window."
      ],
      "metadata": {
        "id": "GX0JSvizv2H9"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mU_sMdR6rNh0"
      },
      "outputs": [],
      "source": [
        "!pip install numpy pmcx jdata bjdata matplotlib"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Step 1: Import modules\n",
        "\n",
        "---\n",
        "\n",
        "> Only need to run this once at the beginning\n"
      ],
      "metadata": {
        "id": "3FuiSDH_w6Dh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import pmcx\n",
        "import jdata as jd\n",
        "from matplotlib import pyplot as plt\n",
        "\n",
        "pmcx.__version__   # print imported pmcx version number"
      ],
      "metadata": {
        "id": "BdyGGnVPr75P",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 36
        },
        "outputId": "1d5c91e9-124f-42ad-b79e-286e4dfc9402"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'0.2.2'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 2
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Step 2: Verify you have access to NVIDIA GPUs\n",
        "\n",
        "---\n",
        "\n",
        "If you run this on Google Colab, please go to menu **Edit\\Notebook Settings\\Hardware accelerator**, and select \"GPU\". By default, Colab does not enable GPU and the below command will fail."
      ],
      "metadata": {
        "id": "oPiJ3kcF4MK2"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "pmcx.gpuinfo()"
      ],
      "metadata": {
        "id": "LnrnJOcJsQzW",
        "outputId": "2e670bef-0f58-4af6-aff8-2a2770303aef",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "[{'name': 'Tesla T4',\n",
              "  'id': 1,\n",
              "  'devcount': 1,\n",
              "  'major': 7,\n",
              "  'minor': 5,\n",
              "  'globalmem': 15835398144,\n",
              "  'constmem': 65536,\n",
              "  'sharedmem': 49152,\n",
              "  'regcount': 65536,\n",
              "  'clock': 1590000,\n",
              "  'sm': 40,\n",
              "  'core': 2560,\n",
              "  'autoblock': 64,\n",
              "  'autothread': 81920,\n",
              "  'maxgate': 0}]"
            ]
          },
          "metadata": {},
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Step 3: Define simulation (as a dictionary variable)\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "SPu7A3Mm4gBh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cfg = {\n",
        "       'nphoton': 1000000,\n",
        "       'vol':np.ones([60,60,60],dtype='uint8'),\n",
        "       'tstart':0,\n",
        "       'tend':5e-9,\n",
        "       'tstep':5e-9,\n",
        "       'srcpos': [30,30,0],\n",
        "       'srcdir':[0,0,1],\n",
        "       'prop':[[0,0,1,1],[0.005,0.1,0.01,1.37]]\n",
        "       }"
      ],
      "metadata": {
        "id": "lY6fw-SnsWxZ"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "##Step 4: Run pmcx simulation\n",
        "\n",
        "---\n"
      ],
      "metadata": {
        "id": "dAT-8UBA4t62"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "res = pmcx.mcxlab(cfg)\n",
        "# or res=pmcx.run(cfg)"
      ],
      "metadata": {
        "id": "fqo7qGqG4tHv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Step 5: Visualize results\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "kSkaE1Uw40Xr"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.imshow(np.log10(res['flux'][30,:, :]))\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 432
        },
        "id": "7IxvdrdWtyMT",
        "outputId": "0e0f80c1-cc72-4dfc-a26b-d5c95c80c458"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "[encoded data, 22536 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# » Customize how to run pmcx"
      ],
      "metadata": {
        "id": "I7V8Gygw-8AM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Use positional parameters\n",
        "\n",
        "---\n",
        "\n",
        "You can call `pmcx.run` in a more compact way (as a one-liner) by directly passing on each setting as a \"positional parameter\""
      ],
      "metadata": {
        "id": "9dRn-1oy_Erl"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "res = pmcx.run(nphoton=1000000, vol=np.ones([60, 60, 60], dtype='uint8'),\n",
        "               tstart=0, tend=5e-9, tstep=5e-9, srcpos=[30,30,0], srcdir=[0,0,1],\n",
        "               prop=np.array([[0, 0, 1, 1], [0.005, 1, 0.01, 1.37]]))\n",
        "\n",
        "res['flux'].shape"
      ],
      "metadata": {
        "id": "qglKOlSE_Jz4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Use itemized dictionary format\n",
        "\n",
        "---\n",
        "\n",
        "Alternatively, you can also be explicity and define each input field in separate lines. This provides better readability."
      ],
      "metadata": {
        "id": "3raVltQ1_wMF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cfg={}\n",
        "cfg['nphoton']=1e7\n",
        "cfg['vol']=np.ones([60, 60, 60], dtype='uint8')\n",
        "cfg['vol'][20:40, 30:40, 20:30]=2\n",
        "cfg['tstart']=0\n",
        "cfg['tend']=5e-9\n",
        "cfg['tstep']=5e-9\n",
        "cfg['srcpos']=[30,30,0]\n",
        "cfg['srcdir']=[0,0,1]\n",
        "cfg['prop']=[[0, 0, 1, 1], [0.005, 1, 0.01, 1.37], [0.1, 10, 0.9, 1]]\n",
        "\n",
        "res=pmcx.run(cfg)\n",
        "\n",
        "res['stat']"
      ],
      "metadata": {
        "id": "-A3ar2ne_7rF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "When using a `dict` cfg input, it is recommended to call `pmcx.mcxlab(cfg)` instead of `pmcx.run()` because `pmcx.mcxlab()` contains the pre- and post-processing code to make the output `res` similar to those returned by `mcxlab()` in MATLAB/Octave.\n",
        "\n",
        "Specificically, `pmcx.mcxlab()` further processes the `res['detp']` data into directionary with keys `{'detid', 'ppath', 'nscat', 'p', 'v', 'mom', 'iquv', 'w0'}` depends on `cfg['savedetflag']` and related input settings. In addition, `pmcx.mcxlab()` also process `res['traj']` trajectory data into directionary-formatted output for easy processing, including keys like `{'id', 'pos'}`. Please run `help(pmcx.mcxlab)` for more details."
      ],
      "metadata": {
        "id": "6liEDDD-RXSy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "res=pmcx.mcxlab(cfg)\n",
        "res.keys()"
      ],
      "metadata": {
        "id": "IYigemuKRVb6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# » Import and export"
      ],
      "metadata": {
        "id": "r85o-0ZJHRnq"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Export data to human-readable JSON/NeuroJSON formats for easy sharing\n",
        "\n",
        "---\n",
        "The US NIH-funded NeuroJSON Project (https://neurojson.org), also developed and led by MCX's author Dr. Fang, is aimed at breaking the barriers for sharing and reusing scientific data between diverse software and programming environments (such as MATLAB, Python, web etc). The NeuroJSON project envisons that the key for making valuable scientific data long-lasting and easily reusable is **human-readability**, thus has adopted JSON and binary JSON as the primary universal format to store data files. JSON is not only directly human-readable, but also directly supported across almost every programming language and software."
      ],
      "metadata": {
        "id": "hx9jSjo3CK9u"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To read/write JSON/NeuroJSON encoded universal data files, one just need to install and import a tiny module called `jdata` (https://pypi.org/project/jdata/)"
      ],
      "metadata": {
        "id": "z1xhSCuoFA3r"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import jdata as jd"
      ],
      "metadata": {
        "id": "CpZcE964CTzQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Once you imported `jdata`, you can then call the `jdata.encode()` function to serialize arbitrary Python data structures, including lists, dicts, numpy arrays etc, into easy-to-parse JSON-like formats. The binary data is losslessly stored in JSON based annotations.\n",
        "\n",
        "It is highly recommended to take advantage of the built-in data compression capability to save your data in a more compact format without losing precision. A compression codec can be pased by adding a key named `compression` in the 2nd parameter of `encode`. Supported compression methods include: `zlib`, `gzip`, `lzma`, `lz4` etc."
      ],
      "metadata": {
        "id": "9rdk1hE3FUM7"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "jd.encode(res['flux'], {'compression':'zlib','base64':1})"
      ],
      "metadata": {
        "id": "6snCnUoMCZJB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Similarly, one can use `lzma` compression to achieve smaller output size at the costof longer compression time. One can compress any numpy array such as `cfg['vol']` to easy-to-exchange JSON data"
      ],
      "metadata": {
        "id": "e29K7DBxGvnJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "jd.encode(cfg['vol'], {'compression':'lzma','base64':1})"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "g4NFlqiHGmM2",
        "outputId": "8199fe1e-7861-49bc-833c-197b6285ad06"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'_ArrayType_': 'uint8',\n",
              " '_ArraySize_': [60, 60, 60],\n",
              " '_ArrayZipType_': 'lzma',\n",
              " '_ArrayZipSize_': [1, 216000],\n",
              " '_ArrayZipData_': b'XQAAgAD//////////wAA7/v//6O3/0c+SBVyOWFRuJIo5qOGB/nu5B6C0y/FOjwBS7F+yYqKTS+jDdl/puOMIxFT4FkYxXWK4nf4tpR/DGrA3nRJZOLpXFOyBNj3RAyrXw1tRunlw3aIt5ZXp4XZZkuf//9alAAA'}"
            ]
          },
          "metadata": {},
          "execution_count": 33
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "JSON/`jdata` encoded JSON-formatted data can be readily decoded to losslessly restore the original data. One can call `jdata.decode()` to achieve this."
      ],
      "metadata": {
        "id": "BeJQxj4jIEg_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "newflux = jd.decode(jd.encode(res['flux'], {'compression':'zlib'}))\n",
        "\n",
        "print([newflux.shape, newflux.dtype])"
      ],
      "metadata": {
        "id": "so8F5PU_IRXz"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Save mcx output data to JSON/NeuroJSON files\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "of3C2OrRHLBV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "jd.save(res['flux'], 'mcx_flux.json', {'compression':'zlib'})"
      ],
      "metadata": {
        "id": "4l9KiVfD43Aa"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Download .json file to local disk\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "zO4kC-t26OVV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import files\n",
        "files.download('mcx_flux.json')"
      ],
      "metadata": {
        "id": "CYBL-Hf55Kpy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# » Use detectors"
      ],
      "metadata": {
        "id": "6Yx9hNng7I7h"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Setup detectors and detection flags"
      ],
      "metadata": {
        "id": "GfYd007cP51D"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cfg['detpos']=[[30,27,0,1], [30,25,0,1]]   # to detect photons, one must first define detectors\n",
        "cfg['issavedet']=1                         # cfg.issavedet must be set to 1 or True in order to save detected photons\n",
        "cfg['issrcfrom0']=1                        # set this flag to ensure src/det coordinates align with voxel space\n",
        "cfg.keys()"
      ],
      "metadata": {
        "id": "NMAqz61S7Nx3",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c9c6d4cf-9ed0-45a9-e6ba-fe4a2cb38867"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "dict_keys(['nphoton', 'vol', 'tstart', 'tend', 'tstep', 'srcpos', 'srcdir', 'prop', 'detpos', 'issavedet', 'issrcfrom0'])"
            ]
          },
          "metadata": {},
          "execution_count": 9
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Run simulation, detected photons are stored in `res['detp']`"
      ],
      "metadata": {
        "id": "dJ7LO-ixP8l5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "res=pmcx.mcxlab(cfg)\n",
        "res['detp'].keys()"
      ],
      "metadata": {
        "id": "bT-OArmNKBHw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Display histogram of partial-paths\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "ULL6n6huPvMH"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plt.hist(res['detp']['ppath'][:,0], bins=100, range=[0,200]);"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 430
        },
        "id": "K2w3S9BdKBnS",
        "outputId": "616b45e9-c2b0-4b63-f0dd-d97dea000c4f"
      },
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "[encoded data, 13532 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Use `savedetflag` to save multiple output data types\n",
        "\n",
        "---\n",
        "in the below example, we ask mcx to output detected photon specific data according to the letters presented in the `savedetflag` string:\n",
        "\n",
        "- `d`: detected photon ID (starting from 1), width: 1\n",
        "- `p`: partial path for each medium type, width: #media\n",
        "- `x`: detected photon existing position, width: 3\n",
        "\n",
        "other possible flags include\n",
        "- `s`: partial scattering count per medium, width: #media\n",
        "- `v`: existing direction vector, width: 3\n",
        "- `m`: momentum transfer per medium, width: #media\n",
        "- `w`: photon initial weight, width: 1"
      ],
      "metadata": {
        "id": "6RWYId8WQo14"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cfg['savedetflag']='dpx'"
      ],
      "metadata": {
        "id": "PClnQkheOZUY"
      },
      "execution_count": 15,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "res=pmcx.run(cfg)\n",
        "\n",
        "res['detp'].shape"
      ],
      "metadata": {
        "id": "8j6dvd2LPj0v"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "In this case, **it is recommended to call `pmcx.mcxlab(cfg)`** instead of `pmcx.run()` as the `pmcx.mcxlab` function splits the detected photon data into detailed named dictionary fields"
      ],
      "metadata": {
        "id": "ppZhI53VTdl2"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "res=pmcx.mcxlab(cfg)\n",
        "\n",
        "res['detp'].keys()"
      ],
      "metadata": {
        "id": "pkRmeLRATP68"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Plot detected photon existing positions\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "dUsLGjGGUM1E"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# plot photon existing position for Det#1 in red\n",
        "plt.scatter(res['detp']['p'][res['detp']['detid']==1,0], res['detp']['p'][res['detp']['detid']==1,1],\n",
        "            marker='.',color='red');\n",
        "\n",
        "# plot photon existing position for Det#2 in blue\n",
        "plt.scatter(res['detp']['p'][res['detp']['detid']==2,0], res['detp']['p'][res['detp']['detid']==2,1],\n",
        "            marker='.',color='blue');\n",
        "plt.axis('equal');"
      ],
      "metadata": {
        "id": "POapUtwoQmWx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Save both fluence and detected photon data to disk\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "T3Ks52zvs70g"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "jd.save(res, 'mcx_flux_dept.json', {'compression':'zlib'})"
      ],
      "metadata": {
        "id": "UNjjTvk0stWa"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Download output file\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "ZAav2LOPtHge"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import files\n",
        "files.download('mcx_flux_dept.json')"
      ],
      "metadata": {
        "id": "dbKqbo6wtK1U"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# » Save diffuse reflectance"
      ],
      "metadata": {
        "id": "sOwcFCxftUvD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Padding volume by zeros to save diffuse reflectance\n",
        "\n",
        "---\n",
        "\n",
        "mcx can accumulate a photon's remaining weight upon exiting the domain surface. The accumulated distribution can be used to compute \"diffuse reflectance\" (`dref`). However, mcx only accumulates `dref` in 0-valued voxels (i.e. background).\n",
        "\n",
        "That means in order to save `dref` along any medium surface, one must pad a layer of zero voxels outside of this surface.\n",
        "\n",
        "In the below example, we set voxels in z=0 plane to be 0s, so we can capture photons escaping from z=1 plane.\n",
        "\n",
        "The saved \"diffuse reflectance\" are stored under the `dref` key of the output dictonary object."
      ],
      "metadata": {
        "id": "a65jBLS9tdXq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "cfg['vol'][:,:,0]=0\n",
        "cfg['issaveref']=1\n",
        "\n",
        "res=pmcx.mcxlab(cfg)\n",
        "\n",
        "res.keys()"
      ],
      "metadata": {
        "id": "-BCp-gJMtpyX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.imshow(np.log10(res['dref'][:, :, 0]))\n",
        "plt.colorbar()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "rnxygE64umkR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# » Photon replay"
      ],
      "metadata": {
        "id": "Vu13dAO4wUK1"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "For details related to the \"photon replay\" algorithm, please check out our below paper:\n",
        "\n",
        "* Ruoyang Yao, Xavier Intes, Qianqian Fang*, [\"A direct approach to compute Jacobians for diffuse optical tomography using perturbation Monte Carlo-based photon 'replay',\"](https://www.osapublishing.org/boe/abstract.cfm?uri=boe-9-10-4588) Biomed. Optics Express 9(10), 4588-4603, (2018)\n"
      ],
      "metadata": {
        "id": "kO1ths515bm7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Using *photon replay* to build Jacobians"
      ],
      "metadata": {
        "id": "rMYcvp2zwbnv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Step 1 - run baseline simulation to get detected photons (and their seeds) at detectors\n",
        "cfg={}\n",
        "cfg['nphoton']=1e7\n",
        "cfg['vol']=np.ones([60, 60, 60], dtype='uint8')\n",
        "cfg['vol'][20:40, 30:40, 20:30]=2\n",
        "cfg['tstart']=0\n",
        "cfg['tend']=5e-9\n",
        "cfg['tstep']=5e-9\n",
        "cfg['srcpos']=[30,30,0]\n",
        "cfg['srcdir']=[0,0,1]\n",
        "cfg['prop']=[[0, 0, 1, 1], [0.005, 1, 0.01, 1.37], [0.1, 10, 0.9, 1]]\n",
        "cfg['detpos']=[[30,20,0,1]]        # to detect photons, one must first define detectors\n",
        "cfg['issavedet']=1                 # cfg.issavedet must be set to 1 or True in order to save detected photons\n",
        "cfg['issrcfrom0']=1                # set this flag to ensure src/det coordinates align with voxel space\n",
        "cfg['issaveseed']=1   ##!important!# set this flag to store detected photon seed data\n",
        "\n",
        "# run baseline simulation, now seed data are stored in res['seeds']]\n",
        "res=pmcx.mcxlab(cfg)\n",
        "res.keys()"
      ],
      "metadata": {
        "id": "YKkvJA6dw25S"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Step 2 - replay detected photon and output Jacobian\n",
        "cfg_replay=cfg;\n",
        "cfg_replay['seed']=res['seeds']        # one must define cfg['seed'] using the returned seeds\n",
        "cfg_replay['detphotons']=res['detp']['data']   # one must define cfg['detphotons'] using the returned detp data\n",
        "cfg_replay['outputtype']='jacobian'    # tell mcx to output absorption (μ_a) Jacobian\n",
        "\n",
        "# run replay, if done properly, the exact the same set of photons will be \"redetected\"\n",
        "res2=pmcx.mcxlab(cfg_replay)\n",
        "res2.keys()\n",
        "print({'redetection':res2['seeds'].shape, 'initial_detection':res['seeds'].shape})"
      ],
      "metadata": {
        "id": "RL30AahxxlnP",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "19f0c007-59fa-4715-9969-9b5f960d1619"
      },
      "execution_count": 27,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "nphoton: 1e+07\n",
            "tstart: 0\n",
            "tstep: 5e-09\n",
            "tend: 5e-09\n",
            "issrcfrom0: 1\n",
            "srcpos: [30, 30, 0, 1]\n",
            "srcdir: [0, 0, 1, 0]\n",
            "issavedet: 1\n",
            "issaveseed: 1\n",
            "dict_keys(['seeds', 'detp', 'flux', 'stat'])\n",
            "{'redetection': (16, 11217), 'initial_detection': (16, 11218)}\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# plot μ_a Jacobian (sensitivity profile, or photon banana)\n",
        "plt.imshow(np.log10(res2['flux'][30,:, :]))\n",
        "plt.colorbar()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "hNDTPCDIydZC",
        "outputId": "d5b4047a-34ab-4191-b2ec-a623ab70d3c2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 468
        }
      },
      "execution_count": 28,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-28-0c9b4b641cc8>:2: RuntimeWarning: divide by zero encountered in log10\n",
            "  plt.imshow(np.log10(res2['flux'][30,:, :]))\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 2 Axes>"
            ],
            "image/png": "[encoded data, 23228 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# » Define domain by JSON based shape descriptors"
      ],
      "metadata": {
        "id": "1zgZKJ1N1Fa3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Define domain by a series of shape constructs\n",
        "\n",
        "---\n"
      ],
      "metadata": {
        "id": "t470_SAB4Nyq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "shape=[]\n",
        "shape.append({'Grid':{'Tag':1, 'Size':[40,60,30]}})\n",
        "shape.append({'Sphere':{'Tag':2, 'O':[20,30,10], 'R':10}})\n",
        "\n",
        "cfg['shapes']=jd.show({'Shapes':shape}, {'string':True})\n",
        "\n",
        "cfg['prop']=[[0, 0, 1, 1], [0.005, 1, 0.01, 1.37], [0.1, 10, 0.9, 1]]"
      ],
      "metadata": {
        "id": "F8Idr-uA1J5L"
      },
      "execution_count": 29,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "cfg['shapes']"
      ],
      "metadata": {
        "id": "j8L27Jgb4Iyy",
        "outputId": "41284b9f-c521-4b9d-9094-042fe768e0ec",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 36
        }
      },
      "execution_count": 30,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'{\"Shapes\": [{\"Grid\": {\"Tag\": 1, \"Size\": [40, 60, 30]}}, {\"Sphere\": {\"Tag\": 2, \"O\": [20, 30, 10], \"R\": 10}}]}'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 30
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "res=pmcx.mcxlab(cfg)"
      ],
      "metadata": {
        "id": "kWTkg2xd2CEI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plt.imshow(np.log10(res['flux'][:, :, 3]))\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "1vR1ST6W2or4"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}