Skip to content

Commit 61e7157

Browse files
hayakawa16hayakawa16
committed
Got fluence_vs_number_of_photons.ipynb working. But needed to start a new branch because somehow (ugh!) I cleaned Kernel but still results in file I pushed. So needed to add other files I have on this branch.
1 parent 0349a81 commit 61e7157

File tree

3 files changed

+676
-0
lines changed

3 files changed

+676
-0
lines changed

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

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

‎monte-carlo/analog_vs_caw_efficiency.ipynb‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
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": "code",
100+
"execution_count": null,
101+
"id": "3d0619f6-acab-4da4-aa45-d48d990a595a",
102+
"metadata": {},
103+
"outputs": [],
104+
"source": [
105+
"Setup the values for the simulations and plot the results using Plotly\n",
106+
"\n",
107+
"Analog vs CAW"
108+
]
109+
},
110+
{
111+
"cell_type": "code",
112+
"execution_count": null,
113+
"id": "b16d74a4",
114+
"metadata": {},
115+
"outputs": [],
116+
"source": [
117+
"# Setup the detector input for the simulation\n",
118+
"detectorRange = DoubleRange(start=0, stop=10, number=101)\n",
119+
"detectorInput = ROfRhoDetectorInput()\n",
120+
"detectorInput.Rho = detectorRange\n",
121+
"detectorInput.TallySecondMoment = True\n",
122+
"detectorInput.Name = \"ROfRho\"\n",
123+
"detectors = Array.CreateInstance(IDetectorInput,1)\n",
124+
"detectors[0] = detectorInput\n",
125+
"\n",
126+
"# Setup the tissue input for the simulation\n",
127+
"regions = Array.CreateInstance(ITissueRegion, 3)\n",
128+
"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",
129+
"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",
130+
"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",
131+
"\n",
132+
"simulationOptions1 = SimulationOptions()\n",
133+
"simulationOptions1.AbsorptionWeightingType = AbsorptionWeightingType.Analog\n",
134+
"# create a SimulationInput object to define the simulation\n",
135+
"simulationInput1 = SimulationInput()\n",
136+
"simulationInput1.N=10000\n",
137+
"simulationInput1.OutputName = \"MonteCarloROfRho-Analog\"\n",
138+
"simulationInput1.DetectorInputs= detectors\n",
139+
"simulationInput1.Options = simulationOptions1\n",
140+
"simulationInput1.Tissue = MultiLayerTissueInput(regions)\n",
141+
"\n",
142+
"simulationOptions2 = SimulationOptions()\n",
143+
"simulationOptions2.AbsorptionWeightingType = AbsorptionWeightingType.Continuous\n",
144+
"# create a SimulationInput object to define the simulation\n",
145+
"simulationInput2 = SimulationInput()\n",
146+
"simulationInput2.N=10000\n",
147+
"simulationInput2.OutputName = \"MonteCarloROfRho-CAW\"\n",
148+
"simulationInput2.DetectorInputs = detectors\n",
149+
"simulationInput2.Options = simulationOptions2\n",
150+
"simulationInput2.Tissue = MultiLayerTissueInput(regions)\n",
151+
"\n",
152+
"# create the simulations\n",
153+
"simulation1 = MonteCarloSimulation(simulationInput1)\n",
154+
"simulation2 = MonteCarloSimulation(simulationInput2)\n",
155+
"\n",
156+
"# run the simulations\n",
157+
"start_time = time.time()\n",
158+
"simulationOutput1 = simulation1.Run()\n",
159+
"end_time = time.time()\n",
160+
"elapsed_time = end_time - start_time\n",
161+
"print(f\"Elapsed time Analog: {elapsed_time:.6f} seconds\")\n",
162+
"start_time = time.time()\n",
163+
"simulationOutput2 = simulation2.Run()\n",
164+
"end_time = time.time()\n",
165+
"elapsed_time = end_time - start_time\n",
166+
"print(f\"Elapsed time CAW: {elapsed_time:.6f} seconds\")\n",
167+
"\n",
168+
"# determine standard deviation and plot the results using Plotly\n",
169+
"detectorResults1 = Array.CreateInstance(ROfRhoDetector,1)\n",
170+
"detectorResults1[0] = simulationOutput1.ResultsDictionary[\"ROfRho\"]\n",
171+
"Reflectance1 = [r for r in detectorResults1[0].Mean]\n",
172+
"SecondMoment1 = [s for s in detectorResults1[0].SecondMoment]\n",
173+
"StandardDeviation1 = np.sqrt((SecondMoment1 - np.multiply(Reflectance1, Reflectance1)) / simulationInput1.N)\n",
174+
"RelativeError1 = np.divide(StandardDeviation1, Reflectance1)\n",
175+
"detectorMidpoints1 = [mp for mp in detectorRange]\n",
176+
"\n",
177+
"detectorResults2 = Array.CreateInstance(ROfRhoDetector,1)\n",
178+
"detectorResults2[0] = simulationOutput2.ResultsDictionary[\"ROfRho\"]\n",
179+
"Reflectance2 = [r for r in detectorResults2[0].Mean]\n",
180+
"SecondMoment2 = [s for s in detectorResults2[0].SecondMoment]\n",
181+
"StandardDeviation2 = np.sqrt((SecondMoment2 - np.multiply(Reflectance2, Reflectance2)) / simulationInput2.N)\n",
182+
"RelativeError2 = np.divide(StandardDeviation2, Reflectance2)\n",
183+
"detectorMidpoints2 = [mp for mp in detectorRange]\n",
184+
"\n",
185+
"# plot reflectance with 1-sigma error bars and relative error difference\n",
186+
"chart = make_subplots(rows=2, cols=1)\n",
187+
"xLabel = \"ρ [mm]\"\n",
188+
"yLabel = \"log(R(ρ)) [mm-2]\"\n",
189+
"# reflectance with 1-sigma error bars: omit last data point because includes reflectance beyond last rho bin\n",
190+
"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",
191+
"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",
192+
"chart.update_traces(error_y_thickness=1)\n",
193+
"chart.update_layout(xaxis_title=xLabel, yaxis_title=yLabel, xaxis_range=[0,10]) \n",
194+
"chart.update_yaxes(type=\"log\", row=1, col=1)\n",
195+
"# relative error difference\n",
196+
"relativeErrorDifference = RelativeError1 - RelativeError2\n",
197+
"chart.add_trace(go.Scatter(x=detectorMidpoints1[:-2], y=relativeErrorDifference[:-1], mode='lines', showlegend=False), row=2, col=1)\n",
198+
"chart.add_hline(y=0.0, line_dash=\"dash\", line_color=\"black\", row=2, col=1)\n",
199+
"chart['layout']['yaxis2']['title']='Analog RE - CAW RE' \n",
200+
"chart['layout']['xaxis2']['title']=xLabel \n",
201+
"chart['layout']['xaxis2']['range']=[0,10]\n",
202+
"\n",
203+
"chart.show()"
204+
]
205+
}
206+
],
207+
"metadata": {
208+
"kernelspec": {
209+
"display_name": "Python 3 (ipykernel)",
210+
"language": "python",
211+
"name": "python3"
212+
},
213+
"language_info": {
214+
"codemirror_mode": {
215+
"name": "ipython",
216+
"version": 3
217+
},
218+
"file_extension": ".py",
219+
"mimetype": "text/x-python",
220+
"name": "python",
221+
"nbconvert_exporter": "python",
222+
"pygments_lexer": "ipython3",
223+
"version": "3.13.7"
224+
}
225+
},
226+
"nbformat": 4,
227+
"nbformat_minor": 5
228+
}

0 commit comments

Comments
 (0)