Feature/25c add more monte carlo scripts for possible short course use

monte-carlo/analog-vs-caw-efficiency.py

@@ -0,0 +1,120 @@
+# Spatially-Resolved Reflectance Predictions for Analog versus CAW Simulations 
+#
+# Goal: This exercise compares error estimates of spatially-resolved 
+# reflectance using Analog versus Continuous Absorption Weighting (CAW) i
+# simulations. 
+#
+# Import the Operating System so we can access the files for the VTS library
+from pythonnet import load
+load('coreclr')
+import clr
+import os
+import time
+file = '../libraries/Vts.dll'
+clr.AddReference(os.path.abspath(file))
+import numpy as np
+import plotly.graph_objects as go
+from plotly.subplots import make_subplots
+from Vts import *
+from Vts.Common import *
+from Vts.Extensions import *
+from Vts.Modeling.Optimizers import *
+from Vts.Modeling.ForwardSolvers import *
+from Vts.SpectralMapping import *
+from Vts.Factories import *
+from Vts.MonteCarlo import *
+from Vts.MonteCarlo.Sources import *
+from Vts.MonteCarlo.Tissues import *
+from Vts.MonteCarlo.Detectors import *
+from Vts.MonteCarlo.Factories import *
+from Vts.MonteCarlo.PhotonData import *
+from Vts.MonteCarlo.PostProcessing import *
+from System import Array, Object, Double
+# Setup the values for the Analog and CAW simulations and plot the results
+# Setup the detector input for the simulation
+detectorRange = DoubleRange(start=0, stop=10, number=101)
+detectorInput = ROfRhoDetectorInput()
+detectorInput.Rho = detectorRange
+detectorInput.TallySecondMoment = True
+detectorInput.Name = "ROfRho"
+detectors = Array.CreateInstance(IDetectorInput,1)
+detectors[0] = detectorInput
+
+# Setup the tissue input for the simulation
+regions = Array.CreateInstance(ITissueRegion, 3)
+regions[0] = LayerTissueRegion(zRange=DoubleRange(Double.NegativeInfinity, 0.0), op=OpticalProperties(mua=0.0, musp=1E-10, g=1.0, n=1.0)) # air
+regions[1] = LayerTissueRegion(zRange=DoubleRange(0.0, 100.0), op=OpticalProperties(mua=0.01, musp=1.0, g=0.8, n=1.4)) # tissue
+regions[2] = LayerTissueRegion(zRange=DoubleRange(100.0, Double.PositiveInfinity), op=OpticalProperties(mua=0.0, musp=1E-10, g=1.0, n=1.0)) # air
+
+simulationOptions1 = SimulationOptions()
+simulationOptions1.AbsorptionWeightingType = AbsorptionWeightingType.Analog
+# create a SimulationInput object to define the simulation
+simulationInput1 = SimulationInput()
+simulationInput1.N=10000
+simulationInput1.OutputName = "MonteCarloROfRho-Analog"
+simulationInput1.DetectorInputs= detectors
+simulationInput1.Options = simulationOptions1
+simulationInput1.Tissue = MultiLayerTissueInput(regions)
+
+simulationOptions2 = SimulationOptions()
+simulationOptions2.AbsorptionWeightingType = AbsorptionWeightingType.Continuous
+# create a SimulationInput object to define the simulation
+simulationInput2 = SimulationInput()
+simulationInput2.N=10000
+simulationInput2.OutputName = "MonteCarloROfRho-CAW"
+simulationInput2.DetectorInputs = detectors
+simulationInput2.Options = simulationOptions2
+simulationInput2.Tissue = MultiLayerTissueInput(regions)
+
+# create the simulations
+simulation1 = MonteCarloSimulation(simulationInput1)
+simulation2 = MonteCarloSimulation(simulationInput2)
+
+# run the simulations
+start_time = time.time()
+simulationOutput1 = simulation1.Run()
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"Elapsed time Analog: {elapsed_time:.6f} seconds")
+start_time = time.time()
+simulationOutput2 = simulation2.Run()
+end_time = time.time()
+elapsed_time = end_time - start_time
+print(f"Elapsed time CAW: {elapsed_time:.6f} seconds")
+
+# determine standard deviation and plot the results using Plotly
+detectorResults1 = Array.CreateInstance(ROfRhoDetector,1)
+detectorResults1[0] = simulationOutput1.ResultsDictionary["ROfRho"]
+Reflectance1 = [r for r in detectorResults1[0].Mean]
+SecondMoment1 = [s for s in detectorResults1[0].SecondMoment]
+StandardDeviation1 = np.sqrt((SecondMoment1 - np.multiply(Reflectance1, Reflectance1)) / simulationInput1.N)
+RelativeError1 = np.divide(StandardDeviation1, Reflectance1)
+detectorMidpoints1 = [mp for mp in detectorRange]
+
+detectorResults2 = Array.CreateInstance(ROfRhoDetector,1)
+detectorResults2[0] = simulationOutput2.ResultsDictionary["ROfRho"]
+Reflectance2 = [r for r in detectorResults2[0].Mean]
+SecondMoment2 = [s for s in detectorResults2[0].SecondMoment]
+StandardDeviation2 = np.sqrt((SecondMoment2 - np.multiply(Reflectance2, Reflectance2)) / simulationInput2.N)
+RelativeError2 = np.divide(StandardDeviation2, Reflectance2)
+detectorMidpoints2 = [mp for mp in detectorRange]
+
+# plot reflectance with 1-sigma error bars and relative error difference
+chart = make_subplots(rows=2, cols=1)
+xLabel = "ρ [mm]"
+yLabel = "log(R(ρ)) [mm-2]"
+# reflectance with 1-sigma error bars: omit last data point because includes reflectance beyond last rho bin
+chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=Reflectance1[:-1], error_y=dict(type='data',array=StandardDeviation1[:-1],visible=True), mode='markers', name='Analog'), row=1, col=1)
+chart.add_trace(go.Scatter(x=detectorMidpoints2[:-2], y=Reflectance2[:-1], error_y=dict(type='data',array=StandardDeviation2[:-1],visible=True), mode='markers', name='CAW'), row=1, col=1)
+chart.update_traces(error_y_thickness=1)
+chart.update_layout(xaxis_title=xLabel, yaxis_title=yLabel, xaxis_range=[0,10])
+chart.update_yaxes(type="log", row=1, col=1)
+# relative error difference
+relativeErrorDifference = RelativeError1 - RelativeError2
+chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=relativeErrorDifference[:-1], mode='lines', showlegend=False), row=2, col=1)
+chart.add_hline(y=0.0, line_dash="dash", line_color="black", row=2, col=1)
+chart['layout']['yaxis2']['title']='Analog RE - CAW RE'
+chart['layout']['xaxis2']['title']=xLabel
+chart['layout']['xaxis2']['range']=[0,10]
+
+chart.show(renderer="browser")

monte-carlo/analog_vs_caw_efficiency.ipynb

@@ -0,0 +1,228 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "38575b9a-5c9a-4dab-8b74-08c28120efcd",
+   "metadata": {},
+   "source": [
+    "# Spatially-Resolved Reflectance Predictions for Analog versus CAW Simulations\n",
+    "**Carole Hayakawa**\n",
+    "\n",
+    "**August 2025**\n",
+    "\n",
+    "Goal: This exercise compares error estimates of spatially-resolved reflectance using Analog versus Continuous Absorption Weighting (CAW) simulations. \n",
+    "It is assumed that\n",
+    "\n",
+    "* [.NET 8](https://dotnet.microsoft.com/en-us/download/dotnet/8.0) has been installed\n",
+    "\n",
+    "* The latest [VTS libraries](https://github.com/VirtualPhotonics/Vts.Scripting.Python/releases) have been downloaded from the zip file in releases and extracted to the libraries folder"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1c12174d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Import the Operating System so we can access the files for the VTS library\n",
+    "import os\n",
+    "current_directory = os.getcwd()\n",
+    "library_directory = current_directory.replace(\"monte-carlo\", \"libraries\")\n",
+    "vts_path = os.path.join(library_directory, \"Vts.dll\")\n",
+    "#Import Math\n",
+    "import math"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f248374",
+   "metadata": {},
+   "source": [
+    "Use pip to install PythonNet Plotly and Numpy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b08ccbd2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pip install pythonnet plotly numpy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "82d0ef92",
+   "metadata": {},
+   "source": [
+    "Import the Core CLR runtime from PythonNet and add the reference for the VTS library and its dependencies\n",
+    "\n",
+    "Import the namespaces from the Python libraries and the VTS library"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "38947713",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pythonnet import set_runtime\n",
+    "set_runtime(\"coreclr\")\n",
+    "import clr\n",
+    "clr.AddReference(vts_path)\n",
+    "import numpy as np\n",
+    "from plotly.subplots import make_subplots\n",
+    "import plotly.graph_objects as go\n",
+    "from System import Array, Double\n",
+    "import time\n",
+    "from Vts import *\n",
+    "from Vts.Common import *\n",
+    "from Vts.Extensions import *\n",
+    "from Vts.Modeling.Optimizers import *\n",
+    "from Vts.Modeling.ForwardSolvers import *\n",
+    "from Vts.SpectralMapping import *\n",
+    "from Vts.Factories import *\n",
+    "from Vts.MonteCarlo import *\n",
+    "from Vts.MonteCarlo.Sources import *\n",
+    "from Vts.MonteCarlo.Tissues import *\n",
+    "from Vts.MonteCarlo.Detectors import *\n",
+    "from Vts.MonteCarlo.Factories import *\n",
+    "from Vts.MonteCarlo.PhotonData import *\n",
+    "from Vts.MonteCarlo.PostProcessing import *\n",
+    "from System import Array"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3d0619f6-acab-4da4-aa45-d48d990a595a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Setup the values for the simulations and plot the results using Plotly\n",
+    "\n",
+    "Analog vs CAW"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b16d74a4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Setup the detector input for the simulation\n",
+    "detectorRange = DoubleRange(start=0, stop=10, number=101)\n",
+    "detectorInput = ROfRhoDetectorInput()\n",
+    "detectorInput.Rho = detectorRange\n",
+    "detectorInput.TallySecondMoment = True\n",
+    "detectorInput.Name = \"ROfRho\"\n",
+    "detectors = Array.CreateInstance(IDetectorInput,1)\n",
+    "detectors[0] = detectorInput\n",
+    "\n",
+    "# Setup the tissue input for the simulation\n",
+    "regions = Array.CreateInstance(ITissueRegion, 3)\n",
+    "regions[0] = LayerTissueRegion(zRange=DoubleRange(Double.NegativeInfinity, 0.0), op=OpticalProperties(mua=0.0, musp=1E-10, g=1.0, n=1.0)) # air\n",
+    "regions[1] = LayerTissueRegion(zRange=DoubleRange(0.0, 100.0), op=OpticalProperties(mua=0.01, musp=1.0, g=0.8, n=1.4)) # tissue\n",
+    "regions[2] = LayerTissueRegion(zRange=DoubleRange(100.0, Double.PositiveInfinity), op=OpticalProperties(mua=0.0, musp=1E-10, g=1.0, n=1.0)) # air\n",
+    "\n",
+    "simulationOptions1 = SimulationOptions()\n",
+    "simulationOptions1.AbsorptionWeightingType = AbsorptionWeightingType.Analog\n",
+    "# create a SimulationInput object to define the simulation\n",
+    "simulationInput1 = SimulationInput()\n",
+    "simulationInput1.N=10000\n",
+    "simulationInput1.OutputName = \"MonteCarloROfRho-Analog\"\n",
+    "simulationInput1.DetectorInputs= detectors\n",
+    "simulationInput1.Options = simulationOptions1\n",
+    "simulationInput1.Tissue = MultiLayerTissueInput(regions)\n",
+    "\n",
+    "simulationOptions2 = SimulationOptions()\n",
+    "simulationOptions2.AbsorptionWeightingType = AbsorptionWeightingType.Continuous\n",
+    "# create a SimulationInput object to define the simulation\n",
+    "simulationInput2 = SimulationInput()\n",
+    "simulationInput2.N=10000\n",
+    "simulationInput2.OutputName = \"MonteCarloROfRho-CAW\"\n",
+    "simulationInput2.DetectorInputs = detectors\n",
+    "simulationInput2.Options = simulationOptions2\n",
+    "simulationInput2.Tissue = MultiLayerTissueInput(regions)\n",
+    "\n",
+    "# create the simulations\n",
+    "simulation1 = MonteCarloSimulation(simulationInput1)\n",
+    "simulation2 = MonteCarloSimulation(simulationInput2)\n",
+    "\n",
+    "# run the simulations\n",
+    "start_time = time.time()\n",
+    "simulationOutput1 = simulation1.Run()\n",
+    "end_time = time.time()\n",
+    "elapsed_time = end_time - start_time\n",
+    "print(f\"Elapsed time Analog: {elapsed_time:.6f} seconds\")\n",
+    "start_time = time.time()\n",
+    "simulationOutput2 = simulation2.Run()\n",
+    "end_time = time.time()\n",
+    "elapsed_time = end_time - start_time\n",
+    "print(f\"Elapsed time CAW: {elapsed_time:.6f} seconds\")\n",
+    "\n",
+    "# determine standard deviation and plot the results using Plotly\n",
+    "detectorResults1 = Array.CreateInstance(ROfRhoDetector,1)\n",
+    "detectorResults1[0] = simulationOutput1.ResultsDictionary[\"ROfRho\"]\n",
+    "Reflectance1 = [r for r in detectorResults1[0].Mean]\n",
+    "SecondMoment1 = [s for s in detectorResults1[0].SecondMoment]\n",
+    "StandardDeviation1 = np.sqrt((SecondMoment1 - np.multiply(Reflectance1, Reflectance1)) / simulationInput1.N)\n",
+    "RelativeError1 = np.divide(StandardDeviation1, Reflectance1)\n",
+    "detectorMidpoints1 = [mp for mp in detectorRange]\n",
+    "\n",
+    "detectorResults2 = Array.CreateInstance(ROfRhoDetector,1)\n",
+    "detectorResults2[0] = simulationOutput2.ResultsDictionary[\"ROfRho\"]\n",
+    "Reflectance2 = [r for r in detectorResults2[0].Mean]\n",
+    "SecondMoment2 = [s for s in detectorResults2[0].SecondMoment]\n",
+    "StandardDeviation2 = np.sqrt((SecondMoment2 - np.multiply(Reflectance2, Reflectance2)) / simulationInput2.N)\n",
+    "RelativeError2 = np.divide(StandardDeviation2, Reflectance2)\n",
+    "detectorMidpoints2 = [mp for mp in detectorRange]\n",
+    "\n",
+    "# plot reflectance with 1-sigma error bars and relative error difference\n",
+    "chart = make_subplots(rows=2, cols=1)\n",
+    "xLabel = \"ρ [mm]\"\n",
+    "yLabel = \"log(R(ρ)) [mm-2]\"\n",
+    "# reflectance with 1-sigma error bars: omit last data point because includes reflectance beyond last rho bin\n",
+    "chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=Reflectance1[:-1], error_y=dict(type='data',array=StandardDeviation1[:-1],visible=True), mode='markers', name='Analog'), row=1, col=1)\n",
+    "chart.add_trace(go.Scatter(x=detectorMidpoints2[:-2], y=Reflectance2[:-1], error_y=dict(type='data',array=StandardDeviation2[:-1],visible=True), mode='markers', name='CAW'), row=1, col=1)\n",
+    "chart.update_traces(error_y_thickness=1)\n",
+    "chart.update_layout(xaxis_title=xLabel, yaxis_title=yLabel, xaxis_range=[0,10])    \n",
+    "chart.update_yaxes(type=\"log\", row=1, col=1)\n",
+    "# relative error difference\n",
+    "relativeErrorDifference = RelativeError1 - RelativeError2\n",
+    "chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=relativeErrorDifference[:-1], mode='lines', showlegend=False), row=2, col=1)\n",
+    "chart.add_hline(y=0.0, line_dash=\"dash\", line_color=\"black\", row=2, col=1)\n",
+    "chart['layout']['yaxis2']['title']='Analog RE - CAW RE' \n",
+    "chart['layout']['xaxis2']['title']=xLabel \n",
+    "chart['layout']['xaxis2']['range']=[0,10]\n",
+    "\n",
+    "chart.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.13.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}