templates/preproc_tally.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import kdsource as kds\n",
    "import mcpl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Open tally"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input cell:\n",
    "It is recommiended to complete template in a duplicated cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outputfile = \"filename\" # Simulation output file\n",
    "tallyname = \"tally\"     # Tally name\n",
    "spectrum = \"spectrum\"   # Energy decay spectrum (only for activation tallies)\n",
    "\n",
    "tally = kds.T4Tally(outputfile, tallyname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = None # Index of cell of non-plotted variable. None to average all cells.\n",
    "\n",
    "[fig, [scores,errs]] = tally.plot2D(['z', 'x'], cell=cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save as particle list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcplname = None # Photon list file name (None: create default name)\n",
    "\n",
    "mcplname = tally.save_tracks(tracksfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create KDSource"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open particle list\n",
    "\n",
    "mcplfile = mcplfile\n",
    "readformat = \"mcpl\"\n",
    "pt = \"p\"\n",
    "trasl = None\n",
    "rot = None\n",
    "x2z = False\n",
    "set_params = True\n",
    "plist = kds.PList(mcplfile, readformat, pt=pt, trasl=trasl, rot=rot, switch_x2z=x2z, set_params=set_params)\n",
    "\n",
    "# Create Geometry\n",
    "\n",
    "trasl = None\n",
    "rot = None\n",
    "volparams = [tally.grids[0][0], tally.grids[0][-1],\n",
    "             tally.grids[1][0], tally.grids[1][-1],\n",
    "             tally.grids[2][0], tally.grids[2][-1]]\n",
    "metric = kds.Geometry([kds.Lethargy(10), kds.Vol(*volparams), kds.Isotrop()], trasl=trasl, rot=rot)\n",
    "\n",
    "# Create KDSource\n",
    "\n",
    "J = 1 # Activation source intensity [1/s]\n",
    "s = kds.KDSource(plist, metric, bw=\"silv\", J=J)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluate statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# With MCPL\n",
    "\n",
    "file = s.plist.filename\n",
    "mcpl.dump_stats(mcpl.collect_stats(file))\n",
    "mcpl.plot_stats(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# With kds.Stats\n",
    "\n",
    "# Select variables ranges (box)\n",
    "maskfun = kds.Box(None, None)\n",
    "\n",
    "# Weighting function\n",
    "fact_dosim = kds.H10(pt=plist.pt)\n",
    "weightfun = lambda parts: fact_dosim(parts[:,0])\n",
    "\n",
    "parts,ws = plist.get()\n",
    "stats = kds.Stats(parts, ws, weightfun=weightfun, maskfun=maskfun)\n",
    "\n",
    "N,I,err = stats.mean_weight(steps=100)\n",
    "plt.show()\n",
    "N,mn,err = stats.mean(var=1, steps=100)\n",
    "plt.show()\n",
    "N,mn,err = stats.std(var=1, steps=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimize BW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set importance for each variable\n",
    "\n",
    "var_imp = [1E-3, 1,1,1, 1E-3,1E-3,1E-3] # Fit only position\n",
    "parts,ws = s.plist.get(1E4)\n",
    "scaling = s.geom.std(parts=parts, weights=ws)\n",
    "scaling /= var_imp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 354,
     "status": "ok",
     "timestamp": 1634594740542,
     "user": {
      "displayName": "Inti Abbate",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "15235311645725350304"
     },
     "user_tz": 180
    },
    "id": "GKc55zPoouGB"
   },
   "outputs": [],
   "source": [
    "# Number of particles to use for optimization.\n",
    "# Use -1 to use all particles in list (recommended)\n",
    "N = -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "m8yv8nu8oLYH"
   },
   "source": [
    "Choose one of the available bandwidth optimization methods. Recommended method is Method 3 (adaptive MLCV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 523,
     "status": "ok",
     "timestamp": 1634594742089,
     "user": {
      "displayName": "Inti Abbate",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "15235311645725350304"
     },
     "user_tz": 180
    },
    "id": "gFTxFighhq2_",
    "outputId": "afa22c99-0c96-4992-e857-c790d36ca2b5",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Method 1: Silverman's Rule: Simple and fast method.\n",
    "# BW is chosen based on only on the number of particles, and dimension of\n",
    "# geometry.\n",
    "\n",
    "s.bw_method = \"silv\"\n",
    "s.fit(N, scaling=scaling)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 450
    },
    "executionInfo": {
     "elapsed": 48768,
     "status": "ok",
     "timestamp": 1634598202519,
     "user": {
      "displayName": "Inti Abbate",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "15235311645725350304"
     },
     "user_tz": 180
    },
    "id": "Owolp6SzqcHc",
    "outputId": "9d26d939-b4d8-4a45-fb50-b98984530cb9",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Method 2: Non-adaptive Maximum Likelihood Cross-Validation:\n",
    "# Creates a grid of non-adaptive bandwidths and evaluates the\n",
    "# cross-validation scores on each one, which is an indicator of the\n",
    "# quality of the estimation. Selects the bandwidth that optimizes\n",
    "# CV score.\n",
    "\n",
    "# MLCV optimization with reduced N\n",
    "s.bw_method = \"mlcv\"\n",
    "grid = np.logspace(-1,1,20)\n",
    "N_cv = int(1E4) # Use a smaller N to reduce computation times\n",
    "s.fit(N_cv, scaling=scaling, grid=grid)\n",
    "\n",
    "# Adapt optimized bandwidth\n",
    "bw = s.kde.bw\n",
    "dim = s.geom.dim\n",
    "bw *= kds.bw_silv(dim,N)/kds.bw_silv(dim,N_cv) # Apply Silverman factor\n",
    "s = kds.KDSource(plist, geom, bw=bw) # Create new KDSource with adapted BW\n",
    "s.fit(N=N, scaling=scaling)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 943
    },
    "executionInfo": {
     "elapsed": 133707,
     "status": "ok",
     "timestamp": 1634598497325,
     "user": {
      "displayName": "Inti Abbate",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "15235311645725350304"
     },
     "user_tz": 180
    },
    "id": "e_UaDW8Ipdtr",
    "outputId": "dddb7a90-758e-4c10-b990-6bd7abadba5f",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Method 3: Adaptive Maximum Likelihood Cross-Validation:\n",
    "# Creates a grid of adaptive bandwidths and evaluates the\n",
    "# cross-validation scores on each one, which is an indicator of the\n",
    "# quality of the estimation. Selects the bandwidth that optimizes\n",
    "# CV score.\n",
    "# kNN is used to generate the seed adaptive bandwidth.\n",
    "\n",
    "# kNN bandwidth\n",
    "s.bw_method = \"knn\"\n",
    "batch_size = 10000 # Batch size for KNN search\n",
    "k = 10             # Numer of neighbors per batch\n",
    "s.fit(N, scaling=scaling, batch_size=batch_size, k=k)\n",
    "bw_knn = s.kde.bw\n",
    "\n",
    "# MLCV optimization of previously calculated kNN bandwidth\n",
    "s.bw_method = \"mlcv\"\n",
    "N_cv = int(1E4)      # Use a smaller N to reduce computation times\n",
    "seed = bw_knn[:N_cv] # Use kNN BW as seed (first N elements)\n",
    "grid = np.logspace(-1,1,20)\n",
    "s.fit(N_cv, scaling=scaling, seed=seed, grid=grid)\n",
    "bw_cv = s.kde.bw\n",
    "\n",
    "# Extend MLCV optimization to full KNN BW\n",
    "bw_knn_cv = bw_knn * bw_cv[0]/bw_knn[0]   # Apply MLCV factor\n",
    "dim = s.geom.dim\n",
    "bw_knn_cv *= kds.bw_silv(dim,len(bw_knn))/kds.bw_silv(dim,len(bw_cv)) # Apply Silverman factor\n",
    "s = kds.KDSource(plist, geom, bw=bw_knn_cv) # Create new KDSource with full BW\n",
    "s.fit(N=N, scaling=scaling)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save KDSource"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.scaling[0] = 0 # Discrete energies\n",
    "# Isotropic emission\n",
    "s.scaling[4] = np.inf\n",
    "s.scaling[5] = np.inf\n",
    "s.scaling[6] = np.inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmlfilename = None # XML source parameters file name\n",
    "bwfilename = None  # BW file name. Only for adaptive BW\n",
    "\n",
    "xmlfilename = s.save(xmlfilename, bwfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare tracks and resampled"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample particles from KDE source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "resampledfile = tracksfile.split('.')[0]+\"_resampled\"\n",
    "N_resampled = N\n",
    "!kdtool resample $sourcefilename -n $N_resampled -o $resampledfile\n",
    "resampledfile += \".mcpl.gz\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate histograms of the tracks list and resampled list\n",
    "Modify templates to build the desired plots to compare the tracks and KDE distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fast 1D histograms using mcpl library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var = \"x\"\n",
    "\n",
    "# Tracks list histogram\n",
    "bins_tracks,hist = stats_tracks[var][\"hist_bins\"],stats_tracks[var][\"hist\"]\n",
    "hist_tracks /= np.sum(hist_tracks)\n",
    "plt.stairs(hist, edges=bins_tracks, baseline=None, label=\"Original\")\n",
    "\n",
    "# Resampled list histogram\n",
    "bins_kde,hist = stats_resampled[var][\"hist_bins\"],stats_resampled[var][\"hist\"]\n",
    "hist_kde /= np.sum(hist_kde)\n",
    "plt.stairs(hist, edges=bins_kde, baseline=None, label=\"Resampled\")\n",
    "\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.xlabel(var)\n",
    "plt.ylabel(\"J\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1D histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbins = 100\n",
    "bins = np.logspace(-10,10,nbins+1)\n",
    "grid = (bins[1:] + bins[:-1]) / 2\n",
    "\n",
    "# Tracks list histogram\n",
    "mcplfile = mcpl.MCPLFile(plfile) # plfile should be an MCPL file for this to work\n",
    "hist_tracks = np.zeros(nbins)\n",
    "errs_tracks = np.zeros(nbins)\n",
    "I = 0\n",
    "for pb in mcplfile.particle_blocks:\n",
    "    data = pb.ekin\n",
    "    ws = pb.weight\n",
    "    wssum = ws.sum()\n",
    "    I += wssum\n",
    "    if wssum > 0:\n",
    "        hist_tracks += np.histogram(data, bins=bins, weights=ws, density=True)[0] * wssum\n",
    "        errs_tracks += np.histogram(data, bins=bins, weights=ws**2)[0]\n",
    "hist_tracks *= J / I\n",
    "errs_tracks = np.sqrt(errs_tracks) / (bins[1:]-bins[:-1]) * J / I * grid\n",
    "sp=plt.stairs(hist_tracks, edges=bins, baseline=None, label=\"Tracks\")\n",
    "c=sp.get_edgecolor()\n",
    "plt.errorbar(grid, hist_tracks, yerr=errs_tracks, capsize=1, linewidth=1, fmt='none', ecolor=c)\n",
    "\n",
    "# Resampled list histogram\n",
    "mcplfile = mcpl.MCPLFile(resampledfile)\n",
    "hist_kde = np.zeros(nbins)\n",
    "errs_kde = np.zeros(nbins)\n",
    "I = 0\n",
    "for pb in mcplfile.particle_blocks:\n",
    "    data = pb.ekin\n",
    "    ws = pb.weight\n",
    "    wssum = ws.sum()\n",
    "    I += wssum\n",
    "    if wssum > 0:\n",
    "        hist_kde += np.histogram(data, bins=bins, weights=ws, density=True)[0] * wssum\n",
    "        errs_kde += np.histogram(data, bins=bins, weights=ws**2)[0]\n",
    "hist_kde *= J / I\n",
    "errs_kde = np.sqrt(errs_kde) / (bins[1:]-bins[:-1]) * J / I * grid\n",
    "sp=plt.stairs(hist_kde, edges=bins, baseline=None, label=\"Resampled\")\n",
    "c=sp.get_edgecolor()\n",
    "plt.errorbar(grid, hist_kde, yerr=errs_kde, capsize=1, linewidth=1, fmt='none', ecolor=c)\n",
    "\n",
    "plt.grid()\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.legend()\n",
    "plt.xlabel(r\"Energy [MeV]\")\n",
    "plt.ylabel(r\"$J\\ \\left[\\frac{n}{MeV\\ s}\\right]$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2D histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbins1 = 100\n",
    "bins1 = np.linspace(-10,10,nbins1+1)\n",
    "nbins2 = 100\n",
    "bins2 = np.linspace(-10,10,nbins2+1)\n",
    "\n",
    "# Tracks list histogram\n",
    "mcplfile = mcpl.MCPLFile(plfile) # plfile should be an MCPL file for this to work\n",
    "hist_tracks = np.zeros((nbins2,nbins1))\n",
    "I = 0\n",
    "for pb in mcplfile.particle_blocks:\n",
    "    data1 = pb.x\n",
    "    data2 = pb.y\n",
    "    ws = pb.weight\n",
    "    wssum = ws.sum()\n",
    "    I += wssum\n",
    "    if wssum > 0:\n",
    "        h,_,_ = np.histogram2d(data2, data1, bins=[bins2,bins1], weights=ws, density=True)\n",
    "        hist_tracks += h * wssum\n",
    "hist_tracks *= J / I\n",
    "plt.pcolormesh(bins1, bins2, hist_tracks, cmap='jet', rasterized=True)\n",
    "plt.title(\"Tracks\")\n",
    "plt.xlabel(\"x [cm]\")\n",
    "plt.ylabel(\"y [cm]\")\n",
    "cbar = plt.colorbar(label=r\"$\\frac{n}{cm^2 s}$\")\n",
    "plt.show()\n",
    "\n",
    "# Resampled list histogram\n",
    "mcplfile = mcpl.MCPLFile(resampledfile)\n",
    "hist_kde = np.zeros((nbins2,nbins1))\n",
    "I = 0\n",
    "for pb in mcplfile.particle_blocks:\n",
    "    data1 = pb.x\n",
    "    data2 = pb.y\n",
    "    ws = pb.weight\n",
    "    I += ws.sum()\n",
    "    h,_,_ = np.histogram2d(data2, data1, bins=[bins2,bins1], weights=ws, density=True)\n",
    "    hist_kde += h * ws.sum()\n",
    "hist_kde *= J / I\n",
    "plt.pcolormesh(bins1, bins2, hist_kde, cmap='jet', rasterized=True)\n",
    "plt.title(\"KDE\")\n",
    "plt.xscale(\"log\")\n",
    "plt.xlabel(\"x [cm]\")\n",
    "plt.ylabel(\"y [cm]\")\n",
    "cbar = plt.colorbar(label=r\"$\\frac{n}{cm^2 s}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}