Scaling Jacobians with voxel side length
Hirvi Pauliina
Dear Dr. Fang and all MCX-users,
First of all, thank you for this interesting group and very valuable discussions!
I would like to ask a question related to the Jacobians with respect to mua. As we know, MCX outputs these in the replay mode if we set 'jacobian' as the output type. I use MCX in MATLAB environment, thus MCXLAB.
I have voxel side length other than 1 mm, thus I set cfg.unitinmm = voxel side length. My question is: After obtaining the Jacobians, shouldn't I multiply the voxel-wise values with cfg.unitinmm to get correct units ([mm])? So, jacobian_final = cfg.unitinmm.*jacobian?
This is my interpretation after reading the source codes (mcx_core.cu and mcx_utils.c). MCX takes cfg.unitinmm into account when calculating the total detected weights (in replay.weight, mcx_utils.c line 1097), but when the voxel-wise path lengths are scaled with the weight, the path length unit seems to be in grid unit (Lmove in mcx_core.cu line 778).
In this part of the code, mua and mus are in grid unit, thus path lengths in grid unit are otherwise appropriate. For fluence and flux, the units are fixed during normalization (mcx_core.cu line 1606), but Jacobian is just scaled with unitless detected weight. In my opinion, this gives the sensitivity with respect to mua in grid unit [1/grid unit], thus we need to multiply with voxel side length to get mm as Jacobian unit. Am I correct or what have I missed?
Also, does setting cfg.seed = -1 (negative) in MCXLAB mean that the seeds are (repeatedly) randomly generated?
I would be extremely grateful for comments!
Best regards,
Pauliina Hirvi
Qianqian Fang
Dear Dr. Fang and all MCX-users,
First of all, thank you for this interesting group and very valuable discussions!
I would like to ask a question related to the Jacobians with respect to mua. As we know, MCX outputs these in the replay mode if we set 'jacobian' as the output type. I use MCX in MATLAB environment, thus MCXLAB.
I have voxel side length other than 1 mm, thus I set cfg.unitinmm = voxel side length. My question is: After obtaining the Jacobians, shouldn't I multiply the voxel-wise values with cfg.unitinmm to get correct units ([mm])? So, jacobian_final = cfg.unitinmm.*jacobian?
This is my interpretation after reading the source codes (mcx_core.cu and mcx_utils.c). MCX takes cfg.unitinmm into account when calculating the total detected weights (in replay.weight, mcx_utils.c line 1097), but when the voxel-wise path lengths are scaled with the weight, the path length unit seems to be in grid unit (Lmove in mcx_core.cu line 778).
In this part of the code, mua and mus are in grid unit, thus path lengths in grid unit are otherwise appropriate. For fluence and flux, the units are fixed during normalization (mcx_core.cu line 1606), but Jacobian is just scaled with unitless detected weight. In my opinion, this gives the sensitivity with respect to mua in grid unit [1/grid unit], thus we need to multiply with voxel side length to get mm as Jacobian unit. Am I correct or what have I missed?
good catch! the replay was implemented without considering the voxel size.
I just committed this change, hopefully will fix the issue
https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0
please let me know if you are able to get consistent results when change unitinmm values.
Also, does setting cfg.seed = -1 (negative) in MCXLAB mean that the seeds are (repeatedly) randomly generated?
clock, so the result is not repeatable. If you want to repeatable simulations, you
should manually set it to a large number. By default, it is set to 0x623F9A9E
https://github.com/fangq/mcx/blob/master/src/mcx_utils.c#L123
hope this helps.
Qianqian
I would be extremely grateful for comments!
Best regards,
Pauliina Hirvi
--
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 post to this group, send email to mcx-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mcx-users.
For more options, visit https://groups.google.com/d/optout.
Hirvi Pauliina
Dear Dr. Fang,
Thank you for Your reply, and sorry for my late response.
Correct me if I have misunderstood, but I think that now in the fixed version the jacobian is divided by cfg.unitinmm, instead of multiplying by it. My interpretation was that you should multiply (jacobian_final = cfg.unitinmm.*jacobian), since when
the voxel-wise path lengths are scaled with the weight, the path lengths are in grid unit.
Now in mcx_core.cu we first multiply the scale with cfg.unitinmm, but then perform scale=1.f/scale, and in the function mcx_normalize (mcx_utils.c), the values are multiplied with the scale, thus divided with cfg.unitinmm. So I suggest performing
scale=1.f/scale before multiplication by cfg.unitinmm.
Furthermore, if the user doesn't wish to normalize the values (for example I don't need values divided by total detected weight at this point), then the obtained jacobians are in grid units, not mm. Maybe this could be commented, or then multiplication
by cfg.unitinmm could be part of the general calculation?
Relating to the question I asked about cfg.seed = -1 (negative): With " repeatedly" I referred to the fact that by default MCX reseeds the RNG after 1e7 steps in each thread, thus the values are "repeatedly" randomly generated. I just wasn't sure whether cfg.seed=-1 works in the matlab version. Thanks for confirming this!
Thank you for the help and valuable discussion!
Best regards,
Pauliina Hirvi
Lähetetty: 29. marraskuuta 2017 22:26
Vastaanottaja: mcx-...@googlegroups.com; Hirvi Pauliina
Aihe: Re: [mcx-users] Scaling Jacobians with voxel side length
Qianqian Fang
Dear Dr. Fang,
Thank you for Your reply, and sorry for my late response.
Correct me if I have misunderstood, but I think that now in the fixed version the jacobian is divided by cfg.unitinmm, instead of multiplying by it. My interpretation was that you should multiply (jacobian_final = cfg.unitinmm.*jacobian), since when the voxel-wise path lengths are scaled with the weight, the path lengths are in grid unit.
I am not sure which line you are looking at, but I did multiply cfg.unitinmm to
normalize the Jacobian, see line#1694 of this patch
https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694
Now in mcx_core.cu we first multiply the scale with cfg.unitinmm, but then perform scale=1.f/scale, and in the function mcx_normalize (mcx_utils.c), the values are multiplied with the scale, thus divided with cfg.unitinmm. So I suggest performing scale=1.f/scale before multiplication by cfg.unitinmm.
Furthermore, if the user doesn't wish to normalize the values (for example I don't need values divided by total detected weight at this point), then the obtained jacobians are in grid units, not mm. Maybe this could be commented, or then multiplication by cfg.unitinmm could be part of the general calculation?
you can see all added comments at
https://github.com/fangq/mcx/commit/7e6e081ea4a8ef3d4deca8168a2963dfb2be0878
as you can see from this comment:
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L921-L925
pretty much everything calculated directly from the GPU kernel has a length unit of a
grid size. Because normalization is done in the host, if one skips normalization, I would
anticipate all length units to be in the grid-unit.
Relating to the question I asked about cfg.seed = -1 (negative): With " repeatedly" I referred to the fact that by default MCX reseeds the RNG after 1e7 steps in each thread, thus the values are "repeatedly" randomly generated. I just wasn't sure whether cfg.seed=-1 works in the matlab version. Thanks for confirming this!
https://github.com/fangq/mcx/blob/master/src/mcx_core.h#L129
the default RNG, i.e. xorshift128+, has quite good quality and period length, so
I don't see rejuvenation is needed any more. so, you can see the gpu_rng_reseed()
function for this RNG is empty.
https://github.com/fangq/mcx/blob/master/src/xorshift128p_rand.cu#L110-L114
reseed may be helpful for the old default RNG (Logistic-lattic) to prevent
collapsing (i.e. all RNG states become 0), but not a must because the LL5 RNG
was stable in most cases.
Also, the RNG "rejuvenation" is not exactly using a new seed (indecently from host), but
dynamically calculated (hashed) from the current RNG states to increase entropy.
hope this helps.
Qianqian
Hirvi Pauliina
Dear Dr. Fang,
Thank you for your quick reply!
Yes, the scale is correctly multiplied by cfg.unitinmm on line #1694 of the patch (https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694).
However, right after this on line #1696 the scale is inverted (scale = 1/scale). Thus,
the function mcx_normalize (line #446 here
https://github.com/fangq/mcx/blob/master/src/mcx_utils.c) divides the Jacobian by cfg.unitinmm, though it should multiply. Maybe I misunderstand the line #1696 here? Sorry for this trouble.
Best wishes,
Pauliina Hirvi
Lähetetty: 3. tammikuuta 2018 19:32
Vastaanottaja: Hirvi Pauliina; mcx-...@googlegroups.com
Aihe: Re: VS: [mcx-users] Scaling Jacobians with voxel side length
Dear Dr. Fang,
Thank you for Your reply, and sorry for my late response.
Correct me if I have misunderstood, but I think that now in the fixed version the jacobian is divided by cfg.unitinmm, instead of multiplying by it. My interpretation was that you should multiply (jacobian_final = cfg.unitinmm.*jacobian), since when the voxel-wise path lengths are scaled with the weight, the path lengths are in grid unit.
hi Pauliina
I am not sure which line you are looking at, but I did multiply cfg.unitinmm to
normalize the Jacobian, see line#1694 of this patch
https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694
Now in mcx_core.cu we first multiply the scale with cfg.unitinmm, but then perform scale=1.f/scale, and in the function mcx_normalize (mcx_utils.c), the values are multiplied with the scale, thus divided with cfg.unitinmm. So I suggest performing scale=1.f/scale before multiplication by cfg.unitinmm.
Furthermore, if the user doesn't wish to normalize the values (for example I don't need values divided by total detected weight at this point), then the obtained jacobians are in grid units, not mm. Maybe this could be commented, or then multiplication by cfg.unitinmm could be part of the general calculation?
during the holiday break, I added some detailed comments to the mcx source code,
you can see all added comments at
https://github.com/fangq/mcx/commit/7e6e081ea4a8ef3d4deca8168a2963dfb2be0878
as you can see from this comment:
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L921-L925
pretty much everything calculated directly from the GPU kernel has a length unit of a
grid size. Because normalization is done in the host, if one skips normalization, I would
anticipate all length units to be in the grid-unit.
Relating to the question I asked about cfg.seed = -1 (negative): With " repeatedly" I referred to the fact that by default MCX reseeds the RNG after 1e7 steps in each thread, thus the values are "repeatedly" randomly generated. I just wasn't sure whether cfg.seed=-1 works in the matlab version. Thanks for confirming this!
I did implement the reseed option after 10^7 steps, but currently, this is obsolete.
https://github.com/fangq/mcx/blob/master/src/mcx_core.h#L129
the default RNG, i.e. xorshift128+, has quite good quality and period length, so
I don't see rejuvenation is needed any more. so, you can see the gpu_rng_reseed()
function for this RNG is empty.
https://github.com/fangq/mcx/blob/master/src/xorshift128p_rand.cu#L110-L114
Hirvi Pauliina
Dear Dr.Fang,
Unfortunately I am still wondering about the scaling on lines #1694-1696 on the patch (https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694). I think that on line #1696 the scale is inverted (scale = 1/scale), thus the code divides jacobian by cfg.unitinmm, when it should multiply. I would be gratefult if you could explain!
In addition, I was wondering whether you are planning to enable anisotropic voxels (voxels with different x,y and z side length) in the
near future?
Best regards,
Pauliina Hirvi
Lähetetty: 3. tammikuuta 2018 21:35
Vastaanottaja: Qianqian Fang; mcx-...@googlegroups.com
Aihe: VS: VS: [mcx-users] Scaling Jacobians with voxel side length
Dear Dr. Fang,
Thank you for your quick reply!
Yes, the scale is correctly multiplied by cfg.unitinmm on line #1694 of the patch (https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694).
However, right after this on line #1696 the scale is inverted (scale = 1/scale). Thus,
the function mcx_normalize (line #446 here
https://github.com/fangq/mcx/blob/master/src/mcx_utils.c) divides the Jacobian by cfg.unitinmm, though it should multiply. Maybe I misunderstand the line #1696 here? Sorry for this trouble.
Best wishes,
Pauliina Hirvi
Lähetetty: 3. tammikuuta 2018 19:32
Vastaanottaja: Hirvi Pauliina; mcx-...@googlegroups.com
Aihe: Re: VS: [mcx-users] Scaling Jacobians with voxel side length
Dear Dr. Fang,
Thank you for Your reply, and sorry for my late response.
Correct me if I have misunderstood, but I think that now in the fixed version the jacobian is divided by cfg.unitinmm, instead of multiplying by it. My interpretation was that you should multiply (jacobian_final = cfg.unitinmm.*jacobian), since when the voxel-wise path lengths are scaled with the weight, the path lengths are in grid unit.
hi Pauliina
I am not sure which line you are looking at, but I did multiply cfg.unitinmm to
normalize the Jacobian, see line#1694 of this patch
https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694
Now in mcx_core.cu we first multiply the scale with cfg.unitinmm, but then perform scale=1.f/scale, and in the function mcx_normalize (mcx_utils.c), the values are multiplied with the scale, thus divided with cfg.unitinmm. So I suggest performing scale=1.f/scale before multiplication by cfg.unitinmm.
Furthermore, if the user doesn't wish to normalize the values (for example I don't need values divided by total detected weight at this point), then the obtained jacobians are in grid units, not mm. Maybe this could be commented, or then multiplication by cfg.unitinmm could be part of the general calculation?
during the holiday break, I added some detailed comments to the mcx source code,
you can see all added comments at
https://github.com/fangq/mcx/commit/7e6e081ea4a8ef3d4deca8168a2963dfb2be0878
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L921-L925
pretty much everything calculated directly from the GPU kernel has a length unit of a
grid size. Because normalization is done in the host, if one skips normalization, I would
anticipate all length units to be in the grid-unit.
Relating to the question I asked about cfg.seed = -1 (negative): With " repeatedly" I referred to the fact that by default MCX reseeds the RNG after 1e7 steps in each thread, thus the values are "repeatedly" randomly generated. I just wasn't sure whether cfg.seed=-1 works in the matlab version. Thanks for confirming this!
I did implement the reseed option after 10^7 steps, but currently, this is obsolete.
https://github.com/fangq/mcx/blob/master/src/mcx_core.h#L129
I don't see rejuvenation is needed any more. so, you can see the gpu_rng_reseed()
function for this RNG is empty.
https://github.com/fangq/mcx/blob/master/src/xorshift128p_rand.cu#L110-L114
Qianqian Fang
Dear Dr.Fang,
Unfortunately I am still wondering about the scaling on lines #1694-1696 on the patch (https://github.com/fangq/mcx/commit/fff057cc550606eed5fedcdf720d1196b26e06a0#diff-0083a506345d0d19caffd23f50b59bcdR1694). I think that on line #1696 the scale is inverted (scale = 1/scale), thus the code divides jacobian by cfg.unitinmm, when it should multiply. I would be gratefult if you could explain!
sorry, it took me a while to get back to you (and understand the problem).
you are right. the normalization factor should multiply unitinmm instead.
previously, my derivation was based on the unit of the Jacobian, which is
dPhi/dmua with 1/mm unit. That's why I divided unitinmm in my first attempt.
now I realized that matching the unit may work for fluence normalization
(because the quantity accumulated is J or count), but it does not work for
Jacobian, because the quantity accumulated is Lmove (grid unit), so I have
to convert that to mm first.
attached is a script that I verified the updated formula.
when moving from 1x1x1mm voxel to 0.5^3 voxels, the Jacobian should drop
by a factor of 8 and the sum should maintain the same.
( what we are solving, i.e. the Jacobian, is not dPhi/dmua itself, but
dPhi/dmua integrated inside each voxel). Given a mua perturbation,
the change in the detector should be constant independent of the
voxel size.
I also updated the code to incorporate this change
https://github.com/fangq/mcx/commit/6b55b4b2a3ed113a8b81953465068c378c253e0a
thanks again for pointing out this inconsistency.
Note: running 5e8 photons in the attached script took 82 seconds for the
0.5^3 voxel domain using Titan V, if it takes too long on your GPU, you
can drop it to 1e7.
In addition, I was wondering whether you are planning to enable anisotropic voxels (voxels with different x,y and z side length) in the near future?
MCX 2018 release (v1.0! finally) and then add big feature changes such
as this one.
if you have been working on this and are willing to share your code, feel
free to send me a pull request from your git repo.
Qianqian