Skip to content

Commit 6bbbd42

Browse files
hayakawa16hayakawa16
authored
Merge pull request #26 from VirtualPhotonics/feature/25c-add-more-monte-carlo-scripts-for-possible-short-course-use
Feature/25c add more monte carlo scripts for possible short course use
2 parents 0349a81 + 48d3007 commit 6bbbd42

File tree

4 files changed

+868
-0
lines changed

4 files changed

+868
-0
lines changed

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

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
1+
# Spatially-Resolved Reflectance Predictions for Analog versus CAW Simulations
2+
#
3+
# Goal: This exercise compares error estimates of spatially-resolved
4+
# reflectance using Analog versus Continuous Absorption Weighting (CAW) i
5+
# simulations.
6+
#
7+
# Import the Operating System so we can access the files for the VTS library
8+
from pythonnet import load
9+
load('coreclr')
10+
import clr
11+
import os
12+
import time
13+
file = '../libraries/Vts.dll'
14+
clr.AddReference(os.path.abspath(file))
15+
import numpy as np
16+
import plotly.graph_objects as go
17+
from plotly.subplots import make_subplots
18+
from Vts import *
19+
from Vts.Common import *
20+
from Vts.Extensions import *
21+
from Vts.Modeling.Optimizers import *
22+
from Vts.Modeling.ForwardSolvers import *
23+
from Vts.SpectralMapping import *
24+
from Vts.Factories import *
25+
from Vts.MonteCarlo import *
26+
from Vts.MonteCarlo.Sources import *
27+
from Vts.MonteCarlo.Tissues import *
28+
from Vts.MonteCarlo.Detectors import *
29+
from Vts.MonteCarlo.Factories import *
30+
from Vts.MonteCarlo.PhotonData import *
31+
from Vts.MonteCarlo.PostProcessing import *
32+
from System import Array, Object, Double
33+
# Setup the values for the Analog and CAW simulations and plot the results
34+
# Setup the detector input for the simulation
35+
detectorRange = DoubleRange(start=0, stop=10, number=101)
36+
detectorInput = ROfRhoDetectorInput()
37+
detectorInput.Rho = detectorRange
38+
detectorInput.TallySecondMoment = True
39+
detectorInput.Name = "ROfRho"
40+
detectors = Array.CreateInstance(IDetectorInput,1)
41+
detectors[0] = detectorInput
42+
43+
# Setup the tissue input for the simulation
44+
regions = Array.CreateInstance(ITissueRegion, 3)
45+
regions[0] = LayerTissueRegion(zRange=DoubleRange(Double.NegativeInfinity, 0.0), op=OpticalProperties(mua=0.0, musp=1E-10, g=1.0, n=1.0)) # air
46+
regions[1] = LayerTissueRegion(zRange=DoubleRange(0.0, 100.0), op=OpticalProperties(mua=0.01, musp=1.0, g=0.8, n=1.4)) # tissue
47+
regions[2] = LayerTissueRegion(zRange=DoubleRange(100.0, Double.PositiveInfinity), op=OpticalProperties(mua=0.0, musp=1E-10, g=1.0, n=1.0)) # air
48+
49+
simulationOptions1 = SimulationOptions()
50+
simulationOptions1.AbsorptionWeightingType = AbsorptionWeightingType.Analog
51+
# create a SimulationInput object to define the simulation
52+
simulationInput1 = SimulationInput()
53+
simulationInput1.N=10000
54+
simulationInput1.OutputName = "MonteCarloROfRho-Analog"
55+
simulationInput1.DetectorInputs= detectors
56+
simulationInput1.Options = simulationOptions1
57+
simulationInput1.Tissue = MultiLayerTissueInput(regions)
58+
59+
simulationOptions2 = SimulationOptions()
60+
simulationOptions2.AbsorptionWeightingType = AbsorptionWeightingType.Continuous
61+
# create a SimulationInput object to define the simulation
62+
simulationInput2 = SimulationInput()
63+
simulationInput2.N=10000
64+
simulationInput2.OutputName = "MonteCarloROfRho-CAW"
65+
simulationInput2.DetectorInputs = detectors
66+
simulationInput2.Options = simulationOptions2
67+
simulationInput2.Tissue = MultiLayerTissueInput(regions)
68+
69+
# create the simulations
70+
simulation1 = MonteCarloSimulation(simulationInput1)
71+
simulation2 = MonteCarloSimulation(simulationInput2)
72+
73+
# run the simulations
74+
start_time = time.time()
75+
simulationOutput1 = simulation1.Run()
76+
end_time = time.time()
77+
elapsed_time = end_time - start_time
78+
print(f"Elapsed time Analog: {elapsed_time:.6f} seconds")
79+
start_time = time.time()
80+
simulationOutput2 = simulation2.Run()
81+
end_time = time.time()
82+
elapsed_time = end_time - start_time
83+
print(f"Elapsed time CAW: {elapsed_time:.6f} seconds")
84+
85+
# determine standard deviation and plot the results using Plotly
86+
detectorResults1 = Array.CreateInstance(ROfRhoDetector,1)
87+
detectorResults1[0] = simulationOutput1.ResultsDictionary["ROfRho"]
88+
Reflectance1 = [r for r in detectorResults1[0].Mean]
89+
SecondMoment1 = [s for s in detectorResults1[0].SecondMoment]
90+
StandardDeviation1 = np.sqrt((SecondMoment1 - np.multiply(Reflectance1, Reflectance1)) / simulationInput1.N)
91+
RelativeError1 = np.divide(StandardDeviation1, Reflectance1)
92+
detectorMidpoints1 = [mp for mp in detectorRange]
93+
94+
detectorResults2 = Array.CreateInstance(ROfRhoDetector,1)
95+
detectorResults2[0] = simulationOutput2.ResultsDictionary["ROfRho"]
96+
Reflectance2 = [r for r in detectorResults2[0].Mean]
97+
SecondMoment2 = [s for s in detectorResults2[0].SecondMoment]
98+
StandardDeviation2 = np.sqrt((SecondMoment2 - np.multiply(Reflectance2, Reflectance2)) / simulationInput2.N)
99+
RelativeError2 = np.divide(StandardDeviation2, Reflectance2)
100+
detectorMidpoints2 = [mp for mp in detectorRange]
101+
102+
# plot reflectance with 1-sigma error bars and relative error difference
103+
chart = make_subplots(rows=2, cols=1)
104+
xLabel = "ρ [mm]"
105+
yLabel = "log(R(ρ)) [mm-2]"
106+
# reflectance with 1-sigma error bars: omit last data point because includes reflectance beyond last rho bin
107+
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)
108+
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)
109+
chart.update_traces(error_y_thickness=1)
110+
chart.update_layout(xaxis_title=xLabel, yaxis_title=yLabel, xaxis_range=[0,10])
111+
chart.update_yaxes(type="log", row=1, col=1)
112+
# relative error difference
113+
relativeErrorDifference = RelativeError1 - RelativeError2
114+
chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=relativeErrorDifference[:-1], mode='lines', showlegend=False), row=2, col=1)
115+
chart.add_hline(y=0.0, line_dash="dash", line_color="black", row=2, col=1)
116+
chart['layout']['yaxis2']['title']='Analog RE - CAW RE'
117+
chart['layout']['xaxis2']['title']=xLabel
118+
chart['layout']['xaxis2']['range']=[0,10]
119+
120+
chart.show(renderer="browser")

‎monte-carlo/analog_vs_caw_efficiency.ipynb‎

Lines changed: 226 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,226 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "38575b9a-5c9a-4dab-8b74-08c28120efcd",
6+
"metadata": {},
7+
"source": [
8+
"# Spatially-Resolved Reflectance Predictions for Analog versus CAW Simulations\n",
9+
"**Carole Hayakawa**\n",
10+
"\n",
11+
"**August 2025**\n",
12+
"\n",
13+
"Goal: This exercise compares error estimates of spatially-resolved reflectance using Analog versus Continuous Absorption Weighting (CAW) simulations. \n",
14+
"It is assumed that\n",
15+
"\n",
16+
"* [.NET 8](https://dotnet.microsoft.com/en-us/download/dotnet/8.0) has been installed\n",
17+
"\n",
18+
"* 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"
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": null,
24+
"id": "1c12174d",
25+
"metadata": {},
26+
"outputs": [],
27+
"source": [
28+
"#Import the Operating System so we can access the files for the VTS library\n",
29+
"import os\n",
30+
"current_directory = os.getcwd()\n",
31+
"library_directory = current_directory.replace(\"monte-carlo\", \"libraries\")\n",
32+
"vts_path = os.path.join(library_directory, \"Vts.dll\")\n",
33+
"#Import Math\n",
34+
"import math"
35+
]
36+
},
37+
{
38+
"cell_type": "markdown",
39+
"id": "7f248374",
40+
"metadata": {},
41+
"source": [
42+
"Use pip to install PythonNet Plotly and Numpy"
43+
]
44+
},
45+
{
46+
"cell_type": "code",
47+
"execution_count": null,
48+
"id": "b08ccbd2",
49+
"metadata": {},
50+
"outputs": [],
51+
"source": [
52+
"pip install pythonnet plotly numpy"
53+
]
54+
},
55+
{
56+
"cell_type": "markdown",
57+
"id": "82d0ef92",
58+
"metadata": {},
59+
"source": [
60+
"Import the Core CLR runtime from PythonNet and add the reference for the VTS library and its dependencies\n",
61+
"\n",
62+
"Import the namespaces from the Python libraries and the VTS library"
63+
]
64+
},
65+
{
66+
"cell_type": "code",
67+
"execution_count": null,
68+
"id": "38947713",
69+
"metadata": {},
70+
"outputs": [],
71+
"source": [
72+
"from pythonnet import set_runtime\n",
73+
"set_runtime(\"coreclr\")\n",
74+
"import clr\n",
75+
"clr.AddReference(vts_path)\n",
76+
"import numpy as np\n",
77+
"from plotly.subplots import make_subplots\n",
78+
"import plotly.graph_objects as go\n",
79+
"from System import Array, Double\n",
80+
"import time\n",
81+
"from Vts import *\n",
82+
"from Vts.Common import *\n",
83+
"from Vts.Extensions import *\n",
84+
"from Vts.Modeling.Optimizers import *\n",
85+
"from Vts.Modeling.ForwardSolvers import *\n",
86+
"from Vts.SpectralMapping import *\n",
87+
"from Vts.Factories import *\n",
88+
"from Vts.MonteCarlo import *\n",
89+
"from Vts.MonteCarlo.Sources import *\n",
90+
"from Vts.MonteCarlo.Tissues import *\n",
91+
"from Vts.MonteCarlo.Detectors import *\n",
92+
"from Vts.MonteCarlo.Factories import *\n",
93+
"from Vts.MonteCarlo.PhotonData import *\n",
94+
"from Vts.MonteCarlo.PostProcessing import *\n",
95+
"from System import Array"
96+
]
97+
},
98+
{
99+
"cell_type": "markdown",
100+
"id": "2385db51-7a57-4753-af03-0ac9d62ecd40",
101+
"metadata": {},
102+
"source": [
103+
"Setup the values for the simulations and plot the results using Plotly\n",
104+
"\n",
105+
"Analog vs CAW"
106+
]
107+
},
108+
{
109+
"cell_type": "code",
110+
"execution_count": null,
111+
"id": "b16d74a4",
112+
"metadata": {},
113+
"outputs": [],
114+
"source": [
115+
"# Setup the detector input for the simulation\n",
116+
"detectorRange = DoubleRange(start=0, stop=10, number=101)\n",
117+
"detectorInput = ROfRhoDetectorInput()\n",
118+
"detectorInput.Rho = detectorRange\n",
119+
"detectorInput.TallySecondMoment = True\n",
120+
"detectorInput.Name = \"ROfRho\"\n",
121+
"detectors = Array.CreateInstance(IDetectorInput,1)\n",
122+
"detectors[0] = detectorInput\n",
123+
"\n",
124+
"# Setup the tissue input for the simulation\n",
125+
"regions = Array.CreateInstance(ITissueRegion, 3)\n",
126+
"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",
127+
"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",
128+
"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",
129+
"\n",
130+
"simulationOptions1 = SimulationOptions()\n",
131+
"simulationOptions1.AbsorptionWeightingType = AbsorptionWeightingType.Analog\n",
132+
"# create a SimulationInput object to define the simulation\n",
133+
"simulationInput1 = SimulationInput()\n",
134+
"simulationInput1.N=10000\n",
135+
"simulationInput1.OutputName = \"MonteCarloROfRho-Analog\"\n",
136+
"simulationInput1.DetectorInputs= detectors\n",
137+
"simulationInput1.Options = simulationOptions1\n",
138+
"simulationInput1.Tissue = MultiLayerTissueInput(regions)\n",
139+
"\n",
140+
"simulationOptions2 = SimulationOptions()\n",
141+
"simulationOptions2.AbsorptionWeightingType = AbsorptionWeightingType.Continuous\n",
142+
"# create a SimulationInput object to define the simulation\n",
143+
"simulationInput2 = SimulationInput()\n",
144+
"simulationInput2.N=10000\n",
145+
"simulationInput2.OutputName = \"MonteCarloROfRho-CAW\"\n",
146+
"simulationInput2.DetectorInputs = detectors\n",
147+
"simulationInput2.Options = simulationOptions2\n",
148+
"simulationInput2.Tissue = MultiLayerTissueInput(regions)\n",
149+
"\n",
150+
"# create the simulations\n",
151+
"simulation1 = MonteCarloSimulation(simulationInput1)\n",
152+
"simulation2 = MonteCarloSimulation(simulationInput2)\n",
153+
"\n",
154+
"# run the simulations\n",
155+
"start_time = time.time()\n",
156+
"simulationOutput1 = simulation1.Run()\n",
157+
"end_time = time.time()\n",
158+
"elapsed_time = end_time - start_time\n",
159+
"print(f\"Elapsed time Analog: {elapsed_time:.6f} seconds\")\n",
160+
"start_time = time.time()\n",
161+
"simulationOutput2 = simulation2.Run()\n",
162+
"end_time = time.time()\n",
163+
"elapsed_time = end_time - start_time\n",
164+
"print(f\"Elapsed time CAW: {elapsed_time:.6f} seconds\")\n",
165+
"\n",
166+
"# determine standard deviation and plot the results using Plotly\n",
167+
"detectorResults1 = Array.CreateInstance(ROfRhoDetector,1)\n",
168+
"detectorResults1[0] = simulationOutput1.ResultsDictionary[\"ROfRho\"]\n",
169+
"Reflectance1 = [r for r in detectorResults1[0].Mean]\n",
170+
"SecondMoment1 = [s for s in detectorResults1[0].SecondMoment]\n",
171+
"StandardDeviation1 = np.sqrt((SecondMoment1 - np.multiply(Reflectance1, Reflectance1)) / simulationInput1.N)\n",
172+
"RelativeError1 = np.divide(StandardDeviation1, Reflectance1)\n",
173+
"detectorMidpoints1 = [mp for mp in detectorRange]\n",
174+
"\n",
175+
"detectorResults2 = Array.CreateInstance(ROfRhoDetector,1)\n",
176+
"detectorResults2[0] = simulationOutput2.ResultsDictionary[\"ROfRho\"]\n",
177+
"Reflectance2 = [r for r in detectorResults2[0].Mean]\n",
178+
"SecondMoment2 = [s for s in detectorResults2[0].SecondMoment]\n",
179+
"StandardDeviation2 = np.sqrt((SecondMoment2 - np.multiply(Reflectance2, Reflectance2)) / simulationInput2.N)\n",
180+
"RelativeError2 = np.divide(StandardDeviation2, Reflectance2)\n",
181+
"detectorMidpoints2 = [mp for mp in detectorRange]\n",
182+
"\n",
183+
"# plot reflectance with 1-sigma error bars and relative error difference\n",
184+
"chart = make_subplots(rows=2, cols=1)\n",
185+
"xLabel = \"ρ [mm]\"\n",
186+
"yLabel = \"log(R(ρ)) [mm-2]\"\n",
187+
"# reflectance with 1-sigma error bars: omit last data point because includes reflectance beyond last rho bin\n",
188+
"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",
189+
"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",
190+
"chart.update_traces(error_y_thickness=1)\n",
191+
"chart.update_layout(xaxis_title=xLabel, yaxis_title=yLabel, xaxis_range=[0,10]) \n",
192+
"chart.update_yaxes(type=\"log\", row=1, col=1)\n",
193+
"# relative error difference\n",
194+
"relativeErrorDifference = RelativeError1 - RelativeError2\n",
195+
"chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=relativeErrorDifference[:-1], mode='lines', showlegend=False), row=2, col=1)\n",
196+
"chart.add_hline(y=0.0, line_dash=\"dash\", line_color=\"black\", row=2, col=1)\n",
197+
"chart['layout']['yaxis2']['title']='Analog RE - CAW RE' \n",
198+
"chart['layout']['xaxis2']['title']=xLabel \n",
199+
"chart['layout']['xaxis2']['range']=[0,10]\n",
200+
"\n",
201+
"chart.show()"
202+
]
203+
}
204+
],
205+
"metadata": {
206+
"kernelspec": {
207+
"display_name": "Python 3 (ipykernel)",
208+
"language": "python",
209+
"name": "python3"
210+
},
211+
"language_info": {
212+
"codemirror_mode": {
213+
"name": "ipython",
214+
"version": 3
215+
},
216+
"file_extension": ".py",
217+
"mimetype": "text/x-python",
218+
"name": "python",
219+
"nbconvert_exporter": "python",
220+
"pygments_lexer": "ipython3",
221+
"version": "3.13.7"
222+
}
223+
},
224+
"nbformat": 4,
225+
"nbformat_minor": 5
226+
}

0 commit comments

Comments
 (0)