Mutiple voxels as illumination sources

[KestutisMa]

Mutiple voxels as illumination sources

[xopto]

This feature is not yet included in the repository (it will be likely included with the next major commit). In the meantime you can try to use the code from the attached files:

1. replace xopto/mcvox/mcsource/voxel.py with the attached voxel.py (implements IsotropicVoxels)
2. replace xopto/mcvox/mcsource/__init__.py with the attached __init__.py (imports IsotropicVoxels)

The attached example isotropicvoxels.py includes a multi-voxel source (IsotropicVoxels) where light is emitted from centrally positioned line of isotropic voxels along the z axis. All voxels have equal intensity in this example. Use the weights property of the source to set the intensities of voxels to custom values (from 0.0 to 1.0). Please note that this is an experimental implementation.

isotropicvoxels.zip
voxel.zip
init.zip

[xopto]

IsotropicVoxels class must be used instead of IsotropicVoxel (this gives ValueError: could not broadcast input array from shape (100,3) into shape (3,)).

A couple of things to check.

1. Ignore the numba warnings (NumbaDeprecationWarning ...)

2. Check that xopto/mcvox/mcsource/voxel.py contains IsotropicVoxels class. The replaced file should define IsotropicVoxel and IsotropicVoxels class. The origin/old voxel.py file only implements/defines IsotropicVoxel.

3. Check the xopto/mcvox/mcsource/__init__.py file which should have a line "from .voxel import IsotropicVoxel, IsotropicVoxels". The original/old file will only imports IsotropicVoxel class: "from .voxel import IsotropicVoxel".
   (The first 2 points should resolve "AttributeError: module 'xopto.mcvox.mcsource' has no attribute 'IsotropicVoxels'. Did you mean: 'IsotropicVoxel'?")

4. The numpy error (AttributeError: module 'numpy' has no attribute 'int'.) can be resolved by modifying line 119 of xopto/mcvox/source/voxel.py from "self._voxel = np.zeros((3,), dtype=np.int)" to "self._voxel = np.zeros((3,), dtype=np.int32)". The problem lies with "np.int" that has been deprecated by numpy. Using np.int32 instead of np.int should solve the issue.

5. Restart python if running the code interactively. This will reload the modified modules.

6. Make sure to use IsotropicVoxels class instead of IsotropicVoxel. The first argument of IsotropicVoxels will take a 2D array of voxel indices (one voxel per row ... [[ind_x1, ind_y1, ind_z1], [ind_x2, ind_y2, ind_z2], ...]). Optionally, use the second argument of IsotropicVoxels to define the intensities of voxels as a vector of floating-point values from 0.0 to 1.0.

Let us know if this solves the issues.

[xopto]

The fluence seems stretched because the axes units in the posted plot are indices instead of mm. Note that the z-axis ends at 25 (index of the last voxel along the z-axis). A plot displayed in mm units would cover a span of 10 mm in all 3 axes and there would be no stretch (x and y axis from -5 to 5 mm, z-axis from 0 to 10 mm). If you are displaying images with matplotlib.pyplot.imshow, use the extent parameter to specify the physical dimensions of the slice/image along the individual axes. Without this information, imshow assumes that pixel/voxel size is equal in all axes and that the size is 1.

The unknown error usually means that the OpenCL engine ran out of resources (e.g. too many Mc instances are created) but there can be many other reasons for this error. Unfortunately, it is not possible to reproduce the error with the provided code. The shape of the index array ((1, 3) or (3,)) is not triggering the error (at least not directly).

[KestutisMa]

Thank you for your prompt reply.
I tried with files you provided, but running isotropicvoxels.py resulted in some errors:

    ~/Doc/J/Mo/pyxopto/examples/mcvox    master !2 ?4  python isotropicvoxels.py                                                                                 ✔  MonteCarloPhoton  
/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/cl/clrng.py:184: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _rng_init(x: np.uint64, a: np.uint32, fora: np.uint32, n_rng: int, xinit: np.uint64) -> int:
Traceback (most recent call last):
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/pyxopto/examples/mcvox/isotropicvoxels.py", line 56, in <module>
    source = mc.mcsource.IsotropicVoxels(src_voxels)
AttributeError: module 'xopto.mcvox.mcsource' has no attribute 'IsotropicVoxels'. Did you mean: 'IsotropicVoxel'?

    ~/Doc/J/Mo/pyxopto/examples/mcvox    master !2 ?4  python isotropicvoxels.py                                                                               1 ✘  MonteCarloPhoton  
/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/cl/clrng.py:184: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _rng_init(x: np.uint64, a: np.uint32, fora: np.uint32, n_rng: int, xinit: np.uint64) -> int:
Traceback (most recent call last):
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/pyxopto/examples/mcvox/isotropicvoxels.py", line 56, in <module>
    source = mc.mcsource.IsotropicVoxel(src_voxels)
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/mcvox/mcsource/voxel.py", line 119, in __init__
    self._voxel = np.zeros((3,), dtype=np.int)
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/numpy/__init__.py", line 319, in __getattr__
    raise AttributeError(__former_attrs__[attr])
AttributeError: module 'numpy' has no attribute 'int'.
`np.int` was a deprecated alias for the builtin `int`. To avoid this error in existing code, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:
    https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'inf'?

    ~/Doc/J/Mo/pyxopto/examples/mcvox    master !2 ?4  python isotropicvoxels.py                                                                               1 ✘  MonteCarloPhoton  
/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/cl/clrng.py:184: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _rng_init(x: np.uint64, a: np.uint32, fora: np.uint32, n_rng: int, xinit: np.uint64) -> int:
Traceback (most recent call last):
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/pyxopto/examples/mcvox/isotropicvoxels.py", line 56, in <module>
    source = mc.mcsource.IsotropicVoxel(src_voxels)
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/mcvox/mcsource/voxel.py", line 120, in __init__
    self._set_voxel(voxel)
  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/mcvox/mcsource/voxel.py", line 125, in _set_voxel
    self._voxel[:] = voxel
ValueError: could not broadcast input array from shape (100,3) into shape (3,)

I tried to fix first two errors, as it was self explainable in code, but stuck at last error:

  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/venv/lib/python3.10/site-packages/pyxopto-0.2.2-py3.10.egg/xopto/mcvox/mcsource/voxel.py", line 125, in _set_voxel
    self._voxel[:] = voxel
ValueError: could not broadcast input array from shape (100,3) into shape (3,)

Kestutis

[KestutisMa]

It worked, many thanks. Problem was because earlier I installed pyxopto with pip, while later cloned with git. So some files were conflicting and I was editing wrong files.

[KestutisMa]

Dear pyxopto Team,
I just observed that if simulation voxels dimensions are not equal, then I got resulting fluence "streched" over some axis, for example if just one voxel (10,10,10) is used as source and domain geometry dimensions are (100,50,25):

        self.voxels = mc.mcgeometry.Voxels(
            mc.mcgeometry.Axis(-5.0e-3, 5.0e-3, 100),
            mc.mcgeometry.Axis(-5.0e-3, 5.0e-3, 50),
            mc.mcgeometry.Axis( 0.0e-3, 10.0e-3, 25)
        )
        src_voxels = np.array([[10,10,10]])
        self.source = mc.mcsource.IsotropicVoxels(voxels=src_voxels)

        self.fluence = mc.mcfluence.Fluence(
            mc.mcfluence.Axis(-5.0e-3, 5.0e-3, 100),
            mc.mcfluence.Axis(-5.0e-3, 5.0e-3, 50),
            mc.mcfluence.Axis( 0.0e-3, 10.0e-3, 25)
        )

this gives fluence look like this, volume plot and volume slice plot of fluence.data array:

Also, trying to simulate when src_voxels = np.array([[10,10,90]]) gives error:

  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/pyxopto/xopto/mcvox/mc.py", line 1002, in run
    ).wait()
pyopencl._cl.LogicError: clWaitForEvents failed: <unknown error -9999>

[KestutisMa]

Thank you very much for quick response, I see now.

The code generating -9999 error:

        voxels = mc.mcgeometry.Voxels(
            mc.mcgeometry.Axis(-5.0e-3, 5.0e-3, 100),
            mc.mcgeometry.Axis(-5.0e-3, 5.0e-3, 50),
            mc.mcgeometry.Axis( 0.0e-3, 10.0e-3, 25)
        )

        materials = mc.mcmaterial.Materials([
            mc.mcmaterial.Material(n=1.0, mua=0.0, mus=0.0, pf=pf),
            mc.mcmaterial.Material(n=1.33, mua=0.01e2, mus=10.0e2, pf=pf),
            mc.mcmaterial.Material(n=1.33, mua=10.0e2, mus=100.0e2, pf=pf)
        ])

        src_voxels = np.array([[10,10,90]])
        source = mc.mcsource.IsotropicVoxels(voxels=src_voxels)

        fluence = mc.mcfluence.Fluence(
            mc.mcfluence.Axis(-5.0e-3, 5.0e-3, 100),
            mc.mcfluence.Axis(-5.0e-3, 5.0e-3, 50),
            mc.mcfluence.Axis( 0.0e-3, 10.0e-3, 25)
        )

        detectors = mc.mcdetector.Detectors(
            top=mc.mcdetector.Cartesian(
                mc.mcdetector.Axis(-5.0e-3, 5e-3, 100)
            ),
            specular=mc.mcdetector.Total()
        )

        mc_obj = mc.Mc(voxels, materials, source, detectors, fluence=fluence,
                cl_devices=cl_device,
                options=[mc.mcoptions.McDebugMode(debug)])

        mc_obj.rmax = 25e-3

        Z, Y, X = mc_obj.voxels.grid
        mc_obj.voxels.material[:] = 2

        nphotons = 1e6 if not debug else 1
        results = None

I just modified line:
src_voxels = np.array([[10,10,10]]) => src_voxels = np.array([[10,10,90]])

[KestutisMa]

I set correct aspect radio of dimensions in mm, and fluence now looks like expected. But I still got -9999 error if I set source voxel to (10,10,90), also got error if (90,10,10): The voxel index (90.0, 10.0, 10.0) is not valid!.

[KestutisMa]

If I set domain dimensions to (10,20,30) then I do not get -9999 error.
But if I set source voxel to (19,2,2) - it results to zero fluence at all voxels, and if I set source voxel to (2,2,19) error appears: ValueError: The voxel index (2.0, 2.0, 19.0) is not valid!

[KestutisMa]

I just modified line 389 in voxel.py:

        return (self.shape[0] > index[2] >= 0) and \
               (self.shape[1] > index[1] >= 0) and \
               (self.shape[2] > index[0] >= 0)

and this error is gone, and fluence looks good:

  File "/home/lab/Documents/Jurevic/MonteCarloPhoton/pyxopto/xopto/mcvox/mcsource/voxel.py", line 389, in cl_pack
    raise ValueError('The voxel index ({}, {}, {}) is not '
ValueError: The voxel index (2.0, 2.0, 19.0) is not valid!

this error was because arrays are addressed zyx in opencl, while python xyz,

Thank you,
Kestutis

[xopto]

There was/is some inconsistency in the implementation of arguments that define voxel indices. Some functions were assuming (ind_x, ind_y, ind_z) while others (ind_z, ind_y, ind_x). This resulted in improper validation of arguments leading to unexpected errors. The API is now pushed towards (ind_x, ind_y, ind_z) except when directly indexing/slicing the numpy arrays. The two attached files should fix the issue. Note that with these modifications the indices of source voxels must be given as (x, y, z).

mcgeometry_voxel.zip
File in mcgeometry_voxel.zip replaces mcvox/mcgeometry/voxel.py.

mcsource_voxel.zip
File in mcsource_voxel.zip replaces mcvox/mcsource/voxel.py.

We are currently preparing a major update of the master branch so updating your local repository with these changes is the simplest solution.

[KestutisMa]

Hi @xopto I still often get error if count of iterations is more than 10-20:

pyopencl._cl.LogicError: clWaitForEvents failed: <unknown error -9999>

I am on Linux, NVIDIA Quadro P2200 GPU, which mem usage is ~1.5GB out of 5 GB, when calculating.
Simulation domain grid: 30x30x10.
Error is occurring at random intervals, that is changing every time I run script and it doesn't depend on GPU mem usage, I closed all other applications.

OpenCL source files fused in 5.410708 ms.
Executing OpenCL code on: [<pyopencl.Device 'Quadro P2200' on 'NVIDIA CUDA' at 0x49fa560>]
OpenCL build options: []
DEBUG:pyopencl.cache:build program: binary cache hit (key: 87c198e3b7cc6ed6231c65f0e7ca068b)
Source code built in 1.098 ms.
OpenCL engine created in 32.833 ms.
DEBUG:pytools.persistent_dict:pyopencl-invoker-cache-v41: disk cache hit [key=6df3abc0553494dec9e292a29b880a21db34e2cf4a9aaa1c1adff46813d93bdc]
DEBUG:pytools.persistent_dict:pyopencl-invoker-cache-v41: disk cache hit [key=fa89a01ab1b76d41f8bfa666b77eb2dd1c9900bd25ff7d205b1158db76f0be01]
DEBUG:pytools.persistent_dict:pyopencl-invoker-cache-v41: disk cache hit [key=ca8250de1f59d1d7b3ddf91bf59e8e7fae4b20e674cbf44396bd757d067211f5]
McKernel processed 1,000,000 packets in 10,240 threads:
    uploaded/built in:     34.567 ms
    executed in      :     34.715 ms
    downloaded in    :      0.608 ms

Where could be the cause of such behaviour?
At the quick glance, isn't 10,240 threads too much?

[KestutisMa]

Interesting note, dimensions 4x4x4 works till end (400 iterations), but 8x8x8 gives same -9999 error after few dozens iterations (sometimes more. sometime less).

[xopto]

Can you prepare and post a minimal example that reproduces this error?
The number of threads seems to be within what one would expect for that GPU.

[KestutisMa]

I did simple test script and it works with dummy input data.
Maybe there is issue in my code, since I use it inside a class.
I will do more extensive testing.

[KestutisMa]

@xopto, After more testing I reproduced error, minimal example here
In three test runs, error appeared after 3, 5, 20 iterations respectively.
Maybe this is because some resources are not released after each subsequent iteration?

[KestutisMa]

Could you please confirm that error is reproducible using minimal example?
I placed all code in single for loop, so Mc object is constructed and destructed every cycle, like new run.
But I still get same error after few dozen cycles. I am wondering if that could be GPU drivers specific problem.
Or maybe OpenCL is not releasing all resources between cycles.

[KestutisMa]

Hey, I've found solution. I just removed voxel sources from last elements in z direction (z = z -1). Now I got no error.
Interesting note, that doesn't work if I do same in x or y direction, despite simulation domain being symmetric in test case.