Replicating Photon Trajectories

[adonald...@gmail.com]

Hello Dr. Fang and all,

I have a simulation volume with a region of interest whose optical properties are different for different simulations. So far, my strategy has been to simply use different simulation volumes (one with the highly-absorbing ROI and one without) for each simulation, but I also realize that I could use the same volume and alter the optical properties of the appropriate materials as needed.

For the sake of troubleshooting, I'd like to run the simulations with each individual detected photon taking the same path through the medium each time. I thought I could do this similar to the approach in demo_mcxlab_replay.m, but that method is still returning partial pathlengths (PPLs) that are in different orders for each simulation. That is, the same PPL values are present for both simulation outputs (both sum to the same value), but these values are in a different order. Example output:

    PPLs (sim 1)  |  PPLs (sim 2)

--------------------------|--------------------------

    1.0e+03 *        |
                             |
    1.0868             |  1.0868
    1.1157             |  1.0205
    1.0205             |  1.1157
    1.1311             |  1.1311
    1.0386             |  1.0386
    1.1042             |  1.1042
    1.0531             |  1.0531
    1.1317             |  1.1317
    1.1444             |  1.1444
    1.1287             |  1.1287
    1.1243             |  1.1243

Any advice on how to go about this? Of course, let me know if you need more info from me; I'll include my code below. Thank you for putting together this great software package!

- Anthony

Code:

% Setup

clear cfg seeds detp detp1 detp2

cfg.nphoton=1e7; % define how many photons to simulate (normally use 1e7)

cfg.vol=uint8(ones(100,60,60)); % generate simulation volume (add ROI later)

cfg.unitinmm = 1;

% Assign Sources/Detectors

rho = 30; % S-D separation (in mm)

if mod(rho,2)==0

sx = 50-rho/2; dx = 50+rho/2;% x-coordinates for src and det

else

sx = 50-rho/2+0.5; dx = 50+rho/2+0.5;% x-coordinates for src and det

end

yplane = 30;% y-dim along which all sources & detectors will be aligned

cfg.srcpos = [sx yplane 1];

cfg.srcdir=[0 0 1]; % source direction

cfg.detpos = [dx yplane 1 1];

% Other

cfg.savedetflag = 'pw'; % want to save PPLs & initial photon packet weights

cfg.gpuid=1; % use the first GPU (default), can use strings of '01' masks to use multiple GPU, such as '0110'

cfg.autopilot=1;

cfg.tstart=0; % define the start time of the simulation (in s)

cfg.tend=5e-9; % define the maximum duration of the photon packet (in s)

cfg.tstep=5e-9;

% Set Optical Properties

g = 0.88;% anisotropy value

cfg.prop=[0 0 1 1 % medium 0: the environment

0.0046 0.83/(1-g) g 1.3 % medium 1: bulk breast tissue

0.0062 0.83/(1-g) g 1.3];% medium 2: ROI (only used by cfgA)

% Run Simulations  

[~,detp,~,seeds] = mcxlab(cfg);

cfg.seed = seeds.data;

cfg.detphotons = detp.data;

[~,detp1] = mcxlab(cfg);

[~,detp2] = mcxlab(cfg);

% Check PPL Alignment

ppl0 = detp.ppath; ppl1 = detp1.ppath; ppl2 = detp2.ppath;

fprintf('\n\tPPL 1\t PPL 2\n')

disp([ppl1(end-10:end,1) ppl2(end-10:end,1)])% each row should be the same

disp(['Total PPLs (sim 1, sim 2): ' num2str(sum(ppl1(:,1))) ' ' num2str(sum(ppl2(:,1)))])% should be equal

disp(['Number of mismatched PPLs: ' num2str(sum(~(round(ppl1,4)==round(ppl2,4))))])% barring some rounding error, this should be 0 for both media

[Qianqian Fang]

hi Anthony,

this is a known behavior. see our comment in mmc's replay example (mcx-replay is similar)

https://github.com/fangq/mmc/blob/v2024.2/mmclab/example/demo_example_replay.m#L61-L68

when running simulations in a multi-threaded environment, the output order could be different - this is because the execution order of the threads many not be exactly reproduced between multiple runs. This makes the output data shuffle the orders.

however, if you call ismember, you could see that the two outputs are actually identical.

all(ismember(ppl1, ppl2, 'rows'))

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/fed1b359-6c07-4830-97da-838aa4125a91n%40googlegroups.com.

[Hirvi Pauliina]

Dear Dr Fang and Anthony,

Just a minor (documentation related) follow-up question on replicating photon trajectories:

I observed that by explicitly setting cfg.seed=0 or cfg.seed=-1, I get different trajectory data over repeated simulations. This is in line with the comment "if this integer is 0 or negative, MCX uses the current system clock ... as the master seed" [https://groups.google.com/g/mcx-users/c/E9VgPb4q62M/m/d6G_9Dr7AgAJ]. If I don't specify the seed explicitly, the trajectories are replicated over iterated simulations.

Based on the documentation, I originally thought that the default value for -E/cfg.seed would be 0 [Github, README.txt] or -1 [https://mcx.space/wiki/index.cgi?Doc/mcx_help], so the simulations would not be replicable by default. The output files also read mcx.seed=0, if seed is not given. But since the results seem to be replicable, I assume the following is the current status: "if user seed is not given, mcx uses 1648335518 as the default seed" [https://groups.google.com/g/mcx-users/c/E9VgPb4q62M/m/d6G_9Dr7AgAJ], making simulations replicable by default?

I didn't yet find where the input seed value is handled in the source code, so I attached a short script for testing the replicability. With Traj_repeat_PH(5,0) and Traj_repeat_PH(5,-1) I get repeat_level=0%, whereas Traj_repeat_PH(5,[]) returns 100%.

Best wishes,

Pauliina

---

From: 'Qianqian Fang' via mcx-users <mcx-...@googlegroups.com>
Sent: Tuesday, August 6, 2024 7:22 AM
To: mcx-...@googlegroups.com <mcx-...@googlegroups.com>; adonald...@gmail.com <adonald...@gmail.com>
Subject: Re: [mcx-users] Replicating Photon Trajectories

 

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/6d4ae236-fb89-41ce-921b-487e676da9c3%40neu.edu.

[Qianqian Fang]

hi Paullina,

glad to hear from you! thanks for pointing out this discrepancy.

I confirm that the default RNG seed is 1648335518 (if user does not do set -E/--seed or cfg.seed).

if user sets -E/--seed to 0 or any negative number (except -999 - used internally for replay .mch input), then mcx uses the current system clock as seed (thus result is not reproducible)

https://github.com/fangq/mcx/blob/v2024.2/src/mcx_core.cu#L3115-L3119

I have updated the documentation to reflect this default behavior with the following two commits

https://github.com/fangq/mcx/commit/fd2db485d9ad29ba4c999b1af1f5d18eaf46a3be

https://github.com/fangq/mcx/commit/3c480b6f69840a3e49046fda943f6bf3e6e762a6

let me know if I miss any other documentation on this.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/AS8P192MB1320CE426879556B8916985EE9BF2%40AS8P192MB1320.EURP192.PROD.OUTLOOK.COM.

[adonald...@gmail.com]

Dr. Fang and Pauliina,

Thank you for the explanation; that makes sense. I'll just run a single simulation with the ROI and change the absorption coefficient as appropriate after the fact. Since the scattering doesn't change in the ROI relative to the background, this should work fine.

And thank you, Pauliina, for that follow-up question! I'd wondered about that but forgot to ask, so I'm glad you brought that up.

Best,

Anthony