From d4db79ce474ceae55664b020f09c01a6cfaddd89 Mon Sep 17 00:00:00 2001 From: Pauliina Hirvi Date: Thu, 9 Apr 2026 17:09:18 +0300 Subject: [PATCH] =?UTF-8?q?Update=20the=20code=20to=20also=20scale=20"b"/"?= =?UTF-8?q?wptof"=20(otWPTOF)=20replay=20output=20with=20the=20current=20s?= =?UTF-8?q?cattering=20coefficient=20in=20SVMC=20mode.=20Fix=20the=20scatt?= =?UTF-8?q?ering=20coefficient=20to=20be=20the=20true=20physical=20scatter?= =?UTF-8?q?ing=20coefficient,=20if=20cfg.unitinmm=20is=20not=201=20and=20m?= =?UTF-8?q?us=20has=20been=20scaled=20(cfg->prop[i].mus=20*=3D=20cfg->unit?= =?UTF-8?q?inmm).=20NOTE:=20Currently=20compilation=20fails=20since=20cfg-?= =?UTF-8?q?>unitinmm=20is=20undefined=20on=20line=202255!=20Where=20is=20i?= =?UTF-8?q?t=20stored=20=E2=80=93=20gcfg=20does=20not=20have=20the=20field?= =?UTF-8?q?=3F=20Otherwise=20compilation=20is=20successful.=20In=20additio?= =?UTF-8?q?n,=20fix=20the=20inner=20derivative=20sign=20when=20converting?= =?UTF-8?q?=20derivatives=20with=20respect=20to=20the=20diffusion=20coeffi?= =?UTF-8?q?cient=20(D)=20to=20derivatives=20with=20respect=20to=20the=20sc?= =?UTF-8?q?attering=20coefficient=20(mus=20or=20musp)=20using=20the=20chai?= =?UTF-8?q?n=20rule=20in=20the=20adjoint=20mode.=20In=20my=20simulations,?= =?UTF-8?q?=20this=20fix=20matches=20the=20signs=20of=20the=20adjoint=20an?= =?UTF-8?q?d=20replay=20scattering=20Jacobians.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/mcx_core.cu | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 3992af05..bb561e79 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -2250,13 +2250,16 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], theta = replayweight[photonidx]; - if (gcfg->outputtype == otRFmus || gcfg->outputtype == otWP) { + if (gcfg->outputtype == otRFmus || gcfg->outputtype == otWP || gcfg->outputtype == otWPTOF) { if (issvmc) { - theta = (prop.mus == 0.f) ? 1.f : __fdividef(theta, prop.mus); + theta = (prop.mus == 0.f) ? 0.f : __fdividef(cfg->unitinmm * theta, prop.mus); } if (gcfg->outputtype == otWP) { tmp0 = theta; + } else if (gcfg->outputtype == otWPTOF) { + tmp0 = photontof[photonidx]; + tmp0 *= theta; } else { ctheta = ppath[gcfg->w0offset + gcfg->srcnum]; stheta = ppath[gcfg->w0offset + gcfg->srcnum + 1]; @@ -2266,7 +2269,6 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], } } else { tmp0 = (gcfg->outputtype == otDCS) ? (1.f - ctheta) : 1.f; - tmp0 = (gcfg->outputtype == otWPTOF) ? photontof[photonidx] : tmp0; tmp0 *= theta; } @@ -4382,8 +4384,8 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { float onemg = 1.f - cfg->prop[medid].g; if (mus > 0.f && onemg > 0.f) { - /** adjoint_musp: J_D * 3*D^2 = J_D / (3*(1-g)^2*mus^2) */ - opscale = 1.f / (3.f * onemg * onemg * mus * mus); + /** adjoint_musp: J_D * (-3*D^2) = -J_D / (3*(1-g)^2*mus^2) */ + opscale = -1.f / (3.f * onemg * onemg * mus * mus); } } @@ -4448,11 +4450,11 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) { float onemg = 1.f - cfg->prop[medid].g; if (mus > 0.f && onemg > 0.f) { - /** D = 1/(3*(1-g)*mus). adjoint_mus: J * 3*D^2*(1-g) = J / (3*(1-g)*mus^2) - * adjoint_musp: J * 3*D^2 = J / (3*(1-g)^2*mus^2) */ + /** D = 1/(3*(1-g)*mus). adjoint_mus: J * (-3*D^2*(1-g)) = - J / (3*(1-g)*mus^2) + * adjoint_musp: J * (-3*D^2) = - J / (3*(1-g)^2*mus^2) */ opscale = (cfg->outputtype == otAdjointMus) - ? 1.f / (3.f * onemg * mus * mus) - : 1.f / (3.f * onemg * onemg * mus * mus); + ? -1.f / (3.f * onemg * mus * mus) + : -1.f / (3.f * onemg * onemg * mus * mus); } }