pmcx fails to read cfg['angleinvcdf'] and cfg['invcdf'] array data

[fangq]

this was reported by @vzoutenbier in the following post in #194

#194 (comment)

the issue seems that the supplied data are only read partially.

Hi Dr. Fang,

Thanks for your quick response. I am using v 0.3.2 so I expect it to have these updates :)

I wrote a script, pasted below, which is able to reproduce my issue across multiple machines. It allows for relatively easy tests. I show some tests below which indicate the issue is within pmcx, or "between the keyboard and the chair".

it seems that only the first entry of a list/numpy array is considered when defining cfg['angleinvcdf']. Below I put in an np.array, list of two entries, and a list of one entry and get the same outputs for all three.

cfg['angleinvcdf'] = np.linspace(1/6,1/5,5) -> (Output A, below) thin cone, not expected
cfg['angleinvcdf'] = [1/6,1/5]                      -> (Output A,) thin cone, not expected
cfg['angleinvcdf'] = **[**1/6**]**                            -> (Output A) thin cone, expected

Output A

This shows that no matter what the length of the input array, it always evaluates the first entry as the only entry.

Changing the first entry of the list does make a difference, we can use it to adjust the emission conical angle:

cfg['angleinvcdf'] = np.linspace(1/6,1/5,5) -> (Output A, above)
cfg['angleinvcdf']= np.linspace(1/12,1/5,5) -> (Output C, below) Changes angle of thin cone

Output C

So, this code only considers the first entry for the output emission.

I also found that the angleinvcdf must be defined in an array, or it evaluates incorrectly:

cfg['angleinvcdf'] = np.array([1/6])              -> (Output A) 
cfg['angleinvcdf'] = **[**1/6**]**                            -> (Output A) thin cone, expected
cfg['angleinvcdf'] = 1/6                               -> (Output B) light emitted exactly along srcdir

Output B

So, it doesnt matter if it's a pythonic list, or a numpy array.

All these results are true, independent of the 4th entry of the source dir.
cfg['srcdir']=[0,0,1,0]
or
cfg['srcdir']=[0,0,1,1]
which should be the difference between interpolation between angles and discrete points- but I don't see any changes.

Something that I notice is that when running the simulation, the terminal output says:

nphoton: 1e+07
tstart: 0
tstep: 5e-09
tend: 5e-09
isreflect: 0
issrcfrom0: 1
srcpos: [30, 30, 15, 1]
srcdir: [0, 0, 1, 0]

which shows that pmcx has appended a fourth entry to the srcpos, this is typically the radius of the source. I can assign this to be anything non zero and the output is the same. If I assign it to 0 the simulation freezes.

Any leads or help here greatly appreciated.

Best regards,
Vincent Zoutenbier

This is basically a verbatim translation from your demo mcxlab example code demo_mcxlab_launchangle.m with pythonic plotting along 3 planes, translated to python. Testing within your group and/or any feedback greatly appreciated.

(replace [hashtag] with a #, these are comments that show as headers in markdown on this website.)

import pmcx
import numpy as np
import matplotlib.pyplot as plt
from Software import tmcx
from copy import deepcopy

gpus=pmcx.gpuinfo()
gpus[0]

# Method 1, setup simulation
cfg={}
cfg['nphoton']=1e7
cfg['vol']=np.ones([61,61,61], dtype='uint8')
cfg['srcpos']=[30,30,15]
cfg['srcdir']=[0,0,1,0]
cfg['prop']=[[0,0,1,1],[0.001,0,1,1]]
cfg['tstart']=0
cfg['tend']=5e-9
cfg['tstep']=5e-9
cfg['isreflect']=0
cfg['issrcfrom0']=1

# define angleinvcdf in launch angle using cfg['angleinvcdf']
cfg['angleinvcdf']=np.linspace(1/6,1/5,5)  # launch angle is uniformly distributed between [pi/6 and pi/5] with interpolation (odd-number length)

# Run the simulation
res=pmcx.run(cfg)

# plot every point with negligable offset
res['flux']=res['flux']+1e-40

# plot cross sections
fig = plt.figure()
gs = fig.add_gridspec(2, 2)
(ax1, ax2), (ax3, ax4) = gs.subplots()
vmin = np.log10(res['flux'].max())-3
vmax = np.log10(res['flux'].max())
ax1.imshow(np.log10(res['flux'][:,:,45]),interpolation='nearest',vmin=vmin, vmax=vmax)
ax1.set_title('Z=45')
ax2.set_axis_off()
ax3.imshow(np.log10(np.squeeze(res['flux'][:,30,:])).T,interpolation='nearest',vmin=vmin, vmax=vmax)
ax3.set_title('Y=30')
ax4.imshow(np.log10(np.squeeze(res['flux'][30,:,:])).T,interpolation='nearest',vmin=vmin, vmax=vmax)
ax4.set_title('X=30')

plt.show()

[fangq]

Something that I notice is that when running the simulation, the terminal output says:

nphoton: 1e+07
tstart: 0
tstep: 5e-09
tend: 5e-09
isreflect: 0
issrcfrom0: 1
srcpos: [30, 30, 15, 1]
srcdir: [0, 0, 1, 0]

which shows that pmcx has appended a fourth entry to the srcpos, this is typically the radius of the source. I can assign this to be anything non zero and the output is the same. If I assign it to 0 the simulation freezes.

Any leads or help here greatly appreciated.

for your above comment - the 4th element of srcpos is the initial weight of each packet. by default, you should set it to 1. if you change it to 0, the launch will stall because it refuses to launch 0-weight packets.

[vzoutenbier]

Great! Thanks for the quick fix and responses!

It looks like it is spreading the photons angles as intended :)

Have a good Wednesday
VZ

[vzoutenbier]

Hi Dr. Fang,

I have a question, that may just be a geometry question rather than a package issue. I am trying to simulate an LED source, I have the invCDF angles calculated, that is, angles which carry equal emission probability between them.

However, including zero, or small angles in general, seems to create an artifact I don't understand. If I run a simulation (same code as above) with the invCDF angles of [0,0.5] I expect an isotropic source that is only launching photons along the source Dir with direction components Vz >=0, and is azimuthally symmertric. That is, it should have an equal probability of emissions between theta 0 and pi/2. However- I see a strong forward propagation component, which appears to give more phtons along Dir. See plot below.

Do you know why this could be?

I thought it may have to do with the way that photons can be generated, if you have an equal probability between 0 and pi/2 AND azimuthal symmetry then there can be a buildup of photons around theta = 0 because all azimuthal angles contribute to the flux at small angles, which is not true of photons with larger azimuthal angles which explore a region of the volume that cannot be approximately accessed by other combinations of theta/phi. if this is true, maybe we need to change the probability function for theta given that more Phi angles can contribute? Im not sure. Any help here is appreciated.

Best regards
Vincent

PS See you at Photonics West?

[fangq]

@vzoutenbier, is the invcdf you computed based on angle or based on u=cos(angle)? here is the help info for invcdf, as you can see it uses u as the variable, not angle

%      cfg.invcdf:     user-specified scattering phase function. To use this, one must define
%                      a vector with monotonically increasing value between -1 and 1
%                      defining the discretized inverse function of the cummulative density function (CDF)
%                      of u=cos(theta), i.e. inv(CDF(u))=inv(CDF(cos(theta))), where theta [0-pi] is the
%                      zenith angle of the scattering event, and u=cos(theta) is between -1 and 1.
%                      Please note that the defined inv(CDF) is relative to u, not to theta. This is because
%                      it is relatively easy to compute as P(u) is the primary form of phase function
%                      see <demo_mcxlab_phasefun.m>

PS See you at Photonics West?

I haven't decided yet - my student, Shijie Yan will give a talk on our most recent works on mmc.

[vzoutenbier]

Totally mised that documentation blurb. Thanks for pointing that out! I wasn't doing any cos() corrections.

I'm a bit confused still.

If I probe what angles are emitted with various angles,

at theta = 0, u = cos(0) = 1
cfg['angleinvcdf']=np.cos(0)

at theta = pi/4, u = cos(pi/4) = 0.707...
cfg['angleinvcdf']=np.cos(np.pi/4)

at theta = pi/3, u = cos(pi/3) = 0.5
cfg['angleinvcdf']=np.cos(np.pi/3)

at theta = 75deg, u = cos(75/180np.pi) = 0.259
cfg['angleinvcdf']=np.cos(75/180np.pi)

at theta = pi/2, u = cos(pi/2) = 0
cfg['angleinvcdf']=np.cos(np.pi/2)

plotting them all together;
cfg['srcdir']=[0,0,1,1] # <- discrete output angles
cfg['angleinvcdf']=[0,0.259,0.5,0.707,1]

this shows the whole space is accessible with theta [0,pi/2] : u=[1,0].

If I input anything that makes u negative, I get the error:
ValueError: cfg.angleinvcdf contains invalid data; it must be a monotonically increasing vector with all values between 0 and 1

which is clear enough. However: i'm not sure that the photons are evenly being distributed through the different angles.

To simulate an LED, I'm interested only in defining forward propagating photons (along the Dir vector, u = 0 through 0.5). I can sample this evenly and see that the angles are indeed evenly spaced out

cfg['srcdir']=[0,0,1,1]
cfg['angleinvcdf']=[0,0.1,0.2,0.3,0.4,0.5]

Why is it that the for 0 direction, forward propagating, the beam appears brighter than the other points, which should have equal weighting?

[vzoutenbier]

Hi Fang, any chance you can comment on why photons are evenly being distributed through the different angles in my last image above? :)

[fangq]

@vzoutenbier, I am so sorry for my late reply. my earlier comment regarding cfg.invcdf was irrelevant, because the field you have question is cfg.angleinvcdf.

There are two modes to control launch angle generation using cfg.angleinvcdf. One mode is discrete, one mode is continuous, which is controlled by the 4th element of cfg.srcdir, which in normal circumstances defines the focal length of the beam.

• when cfg.srcdir(4) = 0 (default), the angles are generated by interpolation - which I assume that's what you wanted.
• when cfg.srcdir(4) = 1, the numbers in cfg.angleinvcdf are assumed to be probabilities at each discrete angle with an interval decided by the length of the vector.

here is the relevant documentation
https://github.com/fangq/mcx/blob/v2024.2/mcxlab/mcxlab.m#L129-L147

here is the demo
https://github.com/fangq/mcx/blob/v2024.2/mcxlab/examples/demo_mcxlab_launchangle.m

in order to create a continuous distribution, you should set the 4th element of srcdir to 0.

[vzoutenbier]

Thanks for these resources!

Using your demo_.._launchangle.m file I was able to make a function to create a function which generates an even distribution photons from the source through the angle of 0 to pi/2. Because of the nonlinearity of cosine, the uniformity of the source output is HIGHLY dependent on the number of points used in the angleinvcdf for small angles. Using a low number of points creates the very strong directly forward scattering artifact, which is what I was chasing in my previous comments. I show this below by showing the output at 5, 50, 500, and 5000 sampling points for angleinvcdf.

Below is my code for anyone else that may have a similar issue, which takes the number of input points for angleinvcdf and creates sampling points for a uniform distribution of source photons between theta_min and theta_max. I know there are other ways to achieve this same effect in MCX, but this is a good basis to use to create arbitrary angle outputs...

function normalized_theta_values = calculateNormalizedTheta(num_points)
% This function calculates normalized theta values between 0 and pi/2
% that are equally spaced in probability described by a cosine function.
% The values are normalized by pi.

% Define the range of theta
theta_min = 0;
theta_max = pi / 2;

% Generate equally spaced probabilities
probabilities = linspace(cos(theta_min), cos(theta_max), num_points);

% Find the corresponding theta values for these probabilities
theta_values = acos(probabilities);

% Normalize by pi
normalized_theta_values = theta_values / pi;

end

[fangq]

@vzoutenbier, thanks for the follow up. are you trying to create a Lambertian source? if so, can you try setting cfg.srcdir(4) to -inf? mcx has built-in support to Lambertian source, see source code

https://github.com/fangq/mcx/blob/v2024.2/src/mcx_core.cu#L1538-L1544
and documentation:
https://github.com/fangq/mcx/blob/v2024.2/mcxlab/mcxlab.m#L78-L79