From 8fff1a13ca88a799958aa8a2c47789514d9aba1c Mon Sep 17 00:00:00 2001 From: hayakawa Date: Wed, 1 May 2024 17:33:37 -0700 Subject: [PATCH 1/6] Bug found in dMCdROfRhodMusDetector. Fixed that and cleaned up code in dMuaDetector with comments. Renamed unit tests for DAW and CAW OneLayerDetectors to include "dMC" in name because added derivative detector tests. --- ....cs => pMCdMCCAWOneLayerDetectorsTests.cs} | 934 ++++++++++-------- ....cs => pMCdMCDAWOneLayerDetectorsTests.cs} | 907 +++++++++-------- .../Detectors/dMCdROfRhodMuaDetector.cs | 53 +- .../Detectors/dMCdROfRhodMusDetector.cs | 115 ++- 4 files changed, 1120 insertions(+), 889 deletions(-) rename src/Vts.Test/MonteCarlo/Detectors/{pMCCAWOneLayerDetectorsTests.cs => pMCdMCCAWOneLayerDetectorsTests.cs} (60%) rename src/Vts.Test/MonteCarlo/Detectors/{pMCDAWOneLayerDetectorsTests.cs => pMCdMCDAWOneLayerDetectorsTests.cs} (63%) diff --git a/src/Vts.Test/MonteCarlo/Detectors/pMCCAWOneLayerDetectorsTests.cs b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCCAWOneLayerDetectorsTests.cs similarity index 60% rename from src/Vts.Test/MonteCarlo/Detectors/pMCCAWOneLayerDetectorsTests.cs rename to src/Vts.Test/MonteCarlo/Detectors/pMCdMCCAWOneLayerDetectorsTests.cs index 99ae5ce49..f67fe8600 100644 --- a/src/Vts.Test/MonteCarlo/Detectors/pMCCAWOneLayerDetectorsTests.cs +++ b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCCAWOneLayerDetectorsTests.cs @@ -1,409 +1,525 @@ -using System; -using System.Collections.Generic; -using NUnit.Framework; -using Vts.Common; -using Vts.IO; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Detectors; -using Vts.MonteCarlo.Helpers; -using Vts.MonteCarlo.Sources; -using Vts.MonteCarlo.Tissues; -using Vts.MonteCarlo.PostProcessing; -using Vts.MonteCarlo.PhotonData; - -namespace Vts.Test.MonteCarlo.Detectors -{ - /// - /// These tests execute perturbation Monte Carlo (pMC) on a continuous absorption weighting (CAW) - /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same results, and - /// 2) the tally results match the linux results given the same seed - /// mersenne twister STANDARD_TEST. The linux results assumes photon passes - /// through specular and deweights photon by specular. This test starts photon - /// inside tissue and then multiplies result by specular deweighting to match - /// linux results. The linux results are generated using the post-processing code in - /// the g_post subdirectory. - /// - [TestFixture] - public class pMCCAWOneLayerDetectorsTests - { - private SimulationInput _referenceInputOneLayerTissue; - private SimulationOutput _referenceOutputOneLayerTissue; - private double _factor; - private pMCDatabase _databaseOneLayerTissue; - - private readonly List _listOfTestGeneratedFiles = new() - { - "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed - "DiffuseReflectanceDatabase.txt", - "CollisionInfoDatabase", - "CollisionInfoDatabase.txt", - "file.txt", // file that captures screen output of MC simulation - }; - - [OneTimeTearDown] - public void Clear_folders_and_files() - { - // make sure databases generated from previous tests are deleted - foreach (var file in _listOfTestGeneratedFiles) - { - FileIO.FileDelete(file); - } - } - /// - /// Define SimulationInput to describe homogeneous one layer tissue - /// - /// - /// - [OneTimeSetUp] - public void Execute_Monte_Carlo() - { - // delete previously generated files - Clear_folders_and_files(); - - var simulationOptions = new SimulationOptions( - 0, - RandomNumberGeneratorType.MersenneTwister, - AbsorptionWeightingType.Continuous, - PhaseFunctionType.HenyeyGreenstein, - new List() { DatabaseType.pMCDiffuseReflectance }, - false, // track statistics - 0.0, // RR threshold -> 0 = no RR performed - 0); - var sourceInput = new DirectionalPointSourceInput( - new Position(0.0, 0.0, 0.0), - new Direction(0.0, 0.0, 1.0), - 1); - var detectorInputs = new List() - { - new ROfRhoDetectorInput() { Rho=new DoubleRange(0.0, 10.0, 101)}, - new ROfRhoAndTimeDetectorInput() { Rho=new DoubleRange(0.0, 10.0, 101), Time=new DoubleRange(0.0, 1.0, 101)} , - new ROfFxDetectorInput() { Fx=new DoubleRange(0.0, 0.5, 11)}, - new ROfFxAndTimeDetectorInput() { Fx=new DoubleRange(0.0, 0.5, 11), Time=new DoubleRange(0.0, 1.0, 101)} - }; - _referenceInputOneLayerTissue = new SimulationInput( - 100, - "", // can't create folder in isolated storage - simulationOptions, - sourceInput, - new MultiLayerTissueInput( - new ITissueRegion[] - { - new LayerTissueRegion( - new DoubleRange(double.NegativeInfinity, 0.0), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), - new LayerTissueRegion( - new DoubleRange(0.0, 20.0), - new OpticalProperties(0.01, 1.0, 0.8, 1.4)), - new LayerTissueRegion( - new DoubleRange(20.0, double.PositiveInfinity), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) - } - ), - detectorInputs); - _factor = 1.0 - Optics.Specular( - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); - _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); - - _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); - } - - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference - /// - [Test] - public void Validate_pMC_CAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoAndTimeDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - Time=new DoubleRange(0.0, 1.0, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - - _referenceOutputOneLayerTissue.R_rt[0, 0]), 0.00000000001); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] * _factor - 92.2411018), 0.0000001); - } - [Test] - public void Validate_pMC_CAW_ROfRho_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 0.00000000001); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.922411018), 0.000000001); - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(rho) - /// - [Test] - public void Validate_pMC_CAW_ROfRho_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 1.013156), 0.000001); - } - /// - /// Test to validate that calling dMC results in not a NaN - /// - [Test] - public void Validate_dMC_CAW_dROfRhodMua_dROfRhodMus_produces_not_NaN_results() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new dMCdROfRhodMuaDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 11), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - }, - new dMCdROfRhodMusDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 11), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); - } - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference - /// - [Test] - public void Validate_pMC_CAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxAndTimeDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - Time=new DoubleRange(0.0, 1.0, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 0.00000000001); - } - [Test] - public void Validate_pMC_CAW_ROfFx_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Real - - _referenceOutputOneLayerTissue.R_fx[0].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Imaginary - - _referenceOutputOneLayerTissue.R_fx[0].Imaginary), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.483629), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0 - - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(fx) - /// - [Test] - public void Validate_pMC_CAW_ROfFx_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from prior run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.303789), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.036982), 0.000001); - } - /// - /// Test to validate that calling dMC results in not a NaN - /// - [Test] - public void Validate_dMC_CAW_dROfRhodMua_produces_not_NaN_results() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new dMCdROfRhodMuaDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - }, - new dMCdROfRhodMusDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); - } - } -} - +using System; +using System.Collections.Generic; +using NUnit.Framework; +using Vts.Common; +using Vts.IO; +using Vts.MonteCarlo; +using Vts.MonteCarlo.Detectors; +using Vts.MonteCarlo.Helpers; +using Vts.MonteCarlo.Sources; +using Vts.MonteCarlo.Tissues; +using Vts.MonteCarlo.PostProcessing; +using Vts.MonteCarlo.PhotonData; + +namespace Vts.Test.MonteCarlo.Detectors +{ + /// + /// These tests execute perturbation Monte Carlo (pMC) on a continuous absorption weighting (CAW) + /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same results, and + /// 2) the tally results match the linux results given the same seed + /// mersenne twister STANDARD_TEST. The linux results assumes photon passes + /// through specular and deweights photon by specular. This test starts photon + /// inside tissue and then multiplies result by specular deweighting to match + /// linux results. The linux results are generated using the post-processing code in + /// the g_post subdirectory. + /// + [TestFixture] + public class pMCdMCCAWOneLayerDetectorsTests + { + private SimulationInput _referenceInputOneLayerTissue; + private SimulationOutput _referenceOutputOneLayerTissue; + private double _factor; + private pMCDatabase _databaseOneLayerTissue; + + private readonly List _listOfTestGeneratedFiles = new() + { + "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed + "DiffuseReflectanceDatabase.txt", + "CollisionInfoDatabase", + "CollisionInfoDatabase.txt", + "file.txt", // file that captures screen output of MC simulation + }; + + [OneTimeTearDown] + public void Clear_folders_and_files() + { + // make sure databases generated from previous tests are deleted + foreach (var file in _listOfTestGeneratedFiles) + { + FileIO.FileDelete(file); + } + } + + /// + /// Define SimulationInput to describe homogeneous one layer tissue + /// + /// + /// + [OneTimeSetUp] + public void Execute_Monte_Carlo() + { + // delete previously generated files + Clear_folders_and_files(); + + var simulationOptions = new SimulationOptions( + 0, + RandomNumberGeneratorType.MersenneTwister, + AbsorptionWeightingType.Continuous, + PhaseFunctionType.HenyeyGreenstein, + new List { DatabaseType.pMCDiffuseReflectance }, + false, // track statistics + 0.0, // RR threshold -> 0 = no RR performed + 0); + var sourceInput = new DirectionalPointSourceInput( + new Position(0.0, 0.0, 0.0), + new Direction(0.0, 0.0, 1.0), + 1); + var detectorInputs = new List + { + new ROfRhoDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101)}, + new ROfRhoRecessedDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101),ZPlane=-1.0}, + new ROfRhoAndTimeDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101), Time=new DoubleRange(0.0, 1.0, 101)} , + new ROfFxDetectorInput { Fx=new DoubleRange(0.0, 0.5, 11)}, + new ROfFxAndTimeDetectorInput { Fx=new DoubleRange(0.0, 0.5, 11), Time=new DoubleRange(0.0, 1.0, 101)} + }; + _referenceInputOneLayerTissue = new SimulationInput( + 100, + "", // can't create folder in isolated storage + simulationOptions, + sourceInput, + new MultiLayerTissueInput( + new ITissueRegion[] + { + new LayerTissueRegion( + new DoubleRange(double.NegativeInfinity, 0.0), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), + new LayerTissueRegion( + new DoubleRange(0.0, 20.0), + new OpticalProperties(0.01, 1.0, 0.8, 1.4)), + new LayerTissueRegion( + new DoubleRange(20.0, double.PositiveInfinity), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) + } + ), + detectorInputs); + _factor = 1.0 - Optics.Specular( + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); + _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); + + _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(rho,time) and that + /// R(rho,time) recessed to a height of 0 are equal + /// + [Test] + public void Validate_pMC_CAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoAndTimeDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0.0, 1.0, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 } + }, + new pMCROfRhoAndTimeRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0, 1.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] * _factor - 92.2411018), 0.0000001); + Assert.AreEqual(88, postProcessedOutput.pMC_R_rt_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rtr[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + Assert.AreEqual(88, postProcessedOutput.pMC_R_rtr_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values (0 perturbation) + /// determines results equal to reference for R(rho) and R(rho) recessed + /// when recessed height=0 for pMC and dMC results + /// + [Test] + public void Validate_pMC_dMC_CAW_ROfRho_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new pMCROfRhoRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.922411018), 0.000000001); + // validation value based on previous run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 30.0061), 0.0001); + Assert.AreEqual(88, postProcessedOutput.pMC_R_r_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rr[0] - + _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + Assert.AreEqual(88, postProcessedOutput.pMC_R_rr_TallyCount); + + // validate derivatives with respect to mua and mus with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.612979) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.207432) < 1e-6); + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(rho) + /// + [Test] + public void Validate_pMC_dMC_CAW_ROfRho_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 1.013156), 0.000001); + Assert.IsTrue(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 37.0997) < 1e-4); + Assert.AreEqual(88, postProcessedOutput.pMC_R_r_TallyCount); + // validate derivative values with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.609284) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.192882) < 1e-6); + // and 2nd moments + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r2[0] - 19.5494) < 1e-4); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r2[0] - 5.47322) < 1e-5); + } + + /// + /// Test to validate that calling dMC results in not a NaN + /// + [Test] + public void Validate_dMC_CAW_dROfRhodMua_dROfRhodMus_produces_not_NaN_results() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new dMCdROfRhodMuaDetectorInput + { + Rho = new DoubleRange(0.0, 10, 11), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + }, + new dMCdROfRhodMusDetectorInput + { + Rho = new DoubleRange(0.0, 10, 11), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + // validation value obtained from linux run using above input and seeded the same + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference + /// + [Test] + public void Validate_pMC_CAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxAndTimeDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + Time=new DoubleRange(0.0, 1.0, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, + PerturbedRegionsIndices=new List { 1 } + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 1e-10); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 1e-10); + } + + [Test] + public void Validate_pMC_CAW_ROfFx_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Real - + _referenceOutputOneLayerTissue.R_fx[0].Real), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Imaginary - + _referenceOutputOneLayerTissue.R_fx[0].Imaginary), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.483629), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0= + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(fx) + /// + [Test] + public void Validate_pMC_CAW_ROfFx_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from prior run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.303789), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.036982), 0.000001); + } + + /// + /// Test to validate that calling dMC results in not a NaN + /// + [Test] + public void Validate_dMC_CAW_dROfRhodMua_produces_not_NaN_results() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new dMCdROfRhodMuaDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + }, + new dMCdROfRhodMusDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); + } + } +} + diff --git a/src/Vts.Test/MonteCarlo/Detectors/pMCDAWOneLayerDetectorsTests.cs b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs similarity index 63% rename from src/Vts.Test/MonteCarlo/Detectors/pMCDAWOneLayerDetectorsTests.cs rename to src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs index 1fd46a3b2..f1c1e33eb 100644 --- a/src/Vts.Test/MonteCarlo/Detectors/pMCDAWOneLayerDetectorsTests.cs +++ b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs @@ -1,411 +1,496 @@ -using System; -using System.Collections.Generic; -using NUnit.Framework; -using Vts.Common; -using Vts.IO; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Detectors; -using Vts.MonteCarlo.Helpers; -using Vts.MonteCarlo.Sources; -using Vts.MonteCarlo.Tissues; -using Vts.MonteCarlo.PostProcessing; -using Vts.MonteCarlo.PhotonData; - -namespace Vts.Test.MonteCarlo.Detectors -{ - /// - /// These tests execute perturbation Monte Carlo (pMC) on a discrete absorption weighting (DAW) - /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same - /// (no perturbation) results, and 2) the tally results match the linux results given the same seed - /// mersenne twister STANDARD_TEST. The linux results assumes photon passes - /// through specular and deweights photon by specular. This test starts photon - /// inside tissue and then multiplies result by specular deweighting to match - /// linux results. The linux results are generated using the post-processing code in - /// the g_post subdirectory. - /// - [TestFixture] - public class pMCDAWOneLayerDetectorsTests - { - private SimulationInput _referenceInputOneLayerTissue; - private SimulationOutput _referenceOutputOneLayerTissue; - private double _factor; - private pMCDatabase _databaseOneLayerTissue; - - readonly List listOfTestFiles = new List() - { - "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed - "DiffuseReflectanceDatabase.txt", - "CollisionInfoDatabase", - "CollisionInfoDatabase.txt", - "file.txt", // file that captures screen output of MC simulation - }; - - [OneTimeTearDown] - public void Clear_folders_and_files() - { - // make sure databases generated from previous tests are deleted - foreach (var file in listOfTestFiles) - { - FileIO.FileDelete(file); - } - } - - /// - /// Define SimulationInput to describe homogeneous one layer case and generate reference database - /// - /// - [OneTimeSetUp] - public void Execute_Monte_Carlo() - { - var simulationOptions = new SimulationOptions( - 0, - RandomNumberGeneratorType.MersenneTwister, - AbsorptionWeightingType.Discrete, - PhaseFunctionType.HenyeyGreenstein, - new List() {DatabaseType.pMCDiffuseReflectance}, - false, // track statistics - 0.0, // RR threshold -> 0 = no RR performed - 0); - var sourceInput = new DirectionalPointSourceInput( - new Position(0.0, 0.0, 0.0), - new Direction(0.0, 0.0, 1.0), - 1); - var detectorInputs = new List() - { - new ROfRhoDetectorInput() {Rho=new DoubleRange(0.0, 10.0, 101)}, - new ROfRhoRecessedDetectorInput() { Rho=new DoubleRange(0.0, 10.0, 101),ZPlane=-1.0}, - new ROfRhoAndTimeDetectorInput() { Rho = new DoubleRange(0.0, 10.0, 101),Time = new DoubleRange(0.0, 1.0, 101)}, - new ROfFxDetectorInput() {Fx=new DoubleRange(0.0, 0.5, 11)}, - new ROfFxAndTimeDetectorInput() { Fx = new DoubleRange(0.0, 0.5, 11),Time = new DoubleRange(0.0, 1.0, 101)} - - }; - _referenceInputOneLayerTissue = new SimulationInput( - 100, - "", // can't create folder in isolated storage - simulationOptions, - sourceInput, - new MultiLayerTissueInput( - new ITissueRegion[] - { - new LayerTissueRegion( - new DoubleRange(double.NegativeInfinity, 0.0), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), - new LayerTissueRegion( - new DoubleRange(0.0, 20.0), - new OpticalProperties(0.01, 1.0, 0.8, 1.4)), - new LayerTissueRegion( - new DoubleRange(20.0, double.PositiveInfinity), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) - } - ), - detectorInputs); - _factor = 1.0 - Optics.Specular( - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); - _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); - - _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); - } - - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(rho,time) and that - /// R(rho,time) recessed to a height of 0 are equal - /// - [Test] - public void Validate_pMC_DAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoAndTimeDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - Time=new DoubleRange(0.0, 1.0, 101), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - }, - new pMCROfRhoAndTimeRecessedDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - Time=new DoubleRange(0, 1.0, 101), - ZPlane=0.0, - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - }, - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - - _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0]*_factor - 61.5238307), 0.0000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_rt_TallyCount); - - // validation value obtained from non-pMC non-recessed run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rtr[0, 0] - - _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); - Assert.AreEqual(89, postProcessedOutput.pMC_R_rtr_TallyCount); - } - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(rho) and R(rho) recessed - /// when Height=0, and R(rho,maxdepth) recessed when Height=0 - /// - [Test] - public void Validate_pMC_DAW_ROfRho_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - }, - new pMCROfRhoRecessedDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - ZPlane=0.0, - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - }, - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0]*_factor - 0.615238307), 0.000000001); - // validation value based on previous run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 20.022918), 0.000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); - - // validation value obtained from non-pMC non-recessed run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rr[0] - - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); - Assert.AreEqual(89, postProcessedOutput.pMC_R_rr_TallyCount); - - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(rho) - /// - [Test] - public void Validate_pMC_DAW_ROfRho_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.7226588), 0.0000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); - } - - /// - /// Test to validate that calling dMC results in not a NaN - /// - [Test] - public void Validate_dMC_DAW_dROfRhodMua_produces_not_NaN_results() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new dMCdROfRhodMuaDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - }, - new dMCdROfRhodMusDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); - Assert.AreEqual(89, postProcessedOutput.dMCdMua_R_r_TallyCount); - Assert.AreEqual(89, postProcessedOutput.dMCdMus_R_r_TallyCount); - } - - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(fx,time) - /// - [Test] - public void Validate_pMC_DAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxAndTimeDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - Time=new DoubleRange(0.0, 1.0, 101), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 0.00000000001); - // validation value obtained from prior run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - 6.858014), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - 0.339772), 0.000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_fxt_TallyCount); - } - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(fx) - /// - [Test] - public void Validate_pMC_DAW_ROfFx_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - _referenceOutputOneLayerTissue.R_fx[1].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - _referenceOutputOneLayerTissue.R_fx[1].Imaginary), 0.00000000001); - // validation value based on previous run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.328865), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.083909), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.467357), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0 - Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(fx) - /// - [Test] - public void Validate_pMC_DAW_ROfFx_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from prior run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.304018), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.029895), 0.000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); - } - - } -} - +using System; +using System.Collections.Generic; +using NUnit.Framework; +using Vts.Common; +using Vts.IO; +using Vts.MonteCarlo; +using Vts.MonteCarlo.Detectors; +using Vts.MonteCarlo.Helpers; +using Vts.MonteCarlo.Sources; +using Vts.MonteCarlo.Tissues; +using Vts.MonteCarlo.PostProcessing; +using Vts.MonteCarlo.PhotonData; + +namespace Vts.Test.MonteCarlo.Detectors +{ + /// + /// These tests execute perturbation Monte Carlo (pMC) on a discrete absorption weighting (DAW) + /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same + /// (no perturbation) results, and 2) the tally results match the linux results given the same seed + /// mersenne twister STANDARD_TEST. The linux results assumes photon passes + /// through specular and deweights photon by specular. This test starts photon + /// inside tissue and then multiplies result by specular deweighting to match + /// linux results. The linux results are generated using the post-processing code in + /// the g_post subdirectory. + /// + [TestFixture] + public class pMCdMCDAWOneLayerDetectorsTests + { + private SimulationInput _referenceInputOneLayerTissue; + private SimulationOutput _referenceOutputOneLayerTissue; + private double _factor; + private pMCDatabase _databaseOneLayerTissue; + + private readonly List _listOfTestFiles = new() + { + "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed + "DiffuseReflectanceDatabase.txt", + "CollisionInfoDatabase", + "CollisionInfoDatabase.txt", + "file.txt", // file that captures screen output of MC simulation + }; + + [OneTimeTearDown] + public void Clear_folders_and_files() + { + // make sure databases generated from previous tests are deleted + foreach (var file in _listOfTestFiles) + { + FileIO.FileDelete(file); + } + } + + /// + /// Define SimulationInput to describe homogeneous one layer case and generate reference database + /// + /// + [OneTimeSetUp] + public void Execute_Monte_Carlo() + { + // delete previously generated files + Clear_folders_and_files(); + + var simulationOptions = new SimulationOptions( + 0, + RandomNumberGeneratorType.MersenneTwister, + AbsorptionWeightingType.Discrete, + PhaseFunctionType.HenyeyGreenstein, + new List {DatabaseType.pMCDiffuseReflectance}, + false, // track statistics + 0.0, // RR threshold -> 0 = no RR performed + 0); + var sourceInput = new DirectionalPointSourceInput( + new Position(0.0, 0.0, 0.0), + new Direction(0.0, 0.0, 1.0), + 1); + var detectorInputs = new List + { + new ROfRhoDetectorInput {Rho=new DoubleRange(0.0, 10.0, 101), TallySecondMoment = true}, + new ROfRhoRecessedDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101),ZPlane=-1.0}, + new ROfRhoAndTimeDetectorInput { Rho = new DoubleRange(0.0, 10.0, 101),Time = new DoubleRange(0.0, 1.0, 101)}, + new ROfFxDetectorInput {Fx=new DoubleRange(0.0, 0.5, 11)}, + new ROfFxAndTimeDetectorInput { Fx = new DoubleRange(0.0, 0.5, 11),Time = new DoubleRange(0.0, 1.0, 101)} + + }; + _referenceInputOneLayerTissue = new SimulationInput( + 100, + "", // can't create folder in isolated storage + simulationOptions, + sourceInput, + new MultiLayerTissueInput( + new ITissueRegion[] + { + new LayerTissueRegion( + new DoubleRange(double.NegativeInfinity, 0.0), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), + new LayerTissueRegion( + new DoubleRange(0.0, 20.0), + new OpticalProperties(0.01, 1.0, 0.8, 1.4)), + new LayerTissueRegion( + new DoubleRange(20.0, double.PositiveInfinity), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) + } + ), + detectorInputs); + _factor = 1.0 - Optics.Specular( + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); + _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); + + _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(rho,time) and that + /// R(rho,time) recessed to a height of 0 are equal + /// + [Test] + public void Validate_pMC_DAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoAndTimeDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0.0, 1.0, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 } + }, + new pMCROfRhoAndTimeRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0, 1.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0]*_factor - 61.5238307), 0.0000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_rt_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rtr[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + Assert.AreEqual(89, postProcessedOutput.pMC_R_rtr_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values (0 perturbation) + /// determines results equal to reference for R(rho) and R(rho) recessed + /// when recessed height=0 for pMC and dMC results + /// + [Test] + public void Validate_pMC_dMC_DAW_ROfRho_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new pMCROfRhoRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0]*_factor - 0.615238307), 0.000000001); + // validation value based on previous run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 20.022918), 0.000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rr[0] - + _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + Assert.AreEqual(89, postProcessedOutput.pMC_R_rr_TallyCount); + + // validate derivatives with respect to mua and mus with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.166157) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.213279) < 1e-6); + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(rho) + /// + [Test] + public void Validate_pMC_dMC_DAW_ROfRho_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.7226588), 0.0000001); + Assert.IsTrue(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 28.10877) < 1e-4); + Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); + // validate derivative values with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.187384) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.235936) < 1e-6); + // and 2nd moments + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r2[0] - 1.80735) < 1e-5); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r2[0] - 5.22420) < 1e-5); + } + + /// + /// Test to validate that calling dMC results in not a NaN + /// + [Test] + public void Validate_dMC_DAW_dROfRhodMua_produces_not_NaN_results() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new dMCdROfRhodMuaDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + }, + new dMCdROfRhodMusDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); + Assert.AreEqual(89, postProcessedOutput.dMCdMua_R_r_TallyCount); + Assert.AreEqual(89, postProcessedOutput.dMCdMus_R_r_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(fx,time) + /// + [Test] + public void Validate_pMC_DAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxAndTimeDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + Time=new DoubleRange(0.0, 1.0, 101), + PerturbedOps=new List + { // perturbed ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 } + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 0.00000000001); + // validation value obtained from prior run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - 6.858014), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - 0.339772), 0.000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_fxt_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(fx) + /// + [Test] + public void Validate_pMC_DAW_ROfFx_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + PerturbedOps=new List + { // perturbed ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - _referenceOutputOneLayerTissue.R_fx[1].Real), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - _referenceOutputOneLayerTissue.R_fx[1].Imaginary), 0.00000000001); + // validation value based on previous run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.328865), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.083909), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.467357), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0 + Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(fx) + /// + [Test] + public void Validate_pMC_DAW_ROfFx_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + // validation value obtained from prior run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.304018), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.029895), 0.000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); + } + + } +} + diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs index 73d04588a..41dda30f9 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs @@ -1,4 +1,4 @@ -using System; + using System; using System.Collections.Generic; using System.Linq; using System.Runtime.Serialization; @@ -138,8 +138,8 @@ public void Initialize(ITissue tissue, Random rng) TallyCount = 0; // if the data arrays are null, create them (only create second moment if TallySecondMoment is true) - Mean = Mean ?? new double[Rho.Count - 1]; - SecondMoment = SecondMoment ?? (TallySecondMoment ? new double[Rho.Count - 1] : null); + Mean ??= new double[Rho.Count - 1]; + SecondMoment ??= (TallySecondMoment ? new double[Rho.Count - 1] : null); // initialize any other necessary class fields here _perturbedOps = PerturbedOps; @@ -155,6 +155,7 @@ public void Initialize(ITissue tissue, Random rng) /// absorption weighting type protected void SetAbsorbAction(AbsorptionWeightingType awt) { + // AbsorptionWeightingType.Analog cannot have derivatives so not a case switch (awt) { // note: pMC is not applied to analog processing, @@ -194,24 +195,18 @@ private double AbsorbContinuous(IList numberOfCollisions, IList pa { var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region + // NOTE: following code only works for single perturbed region because derivative of + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that foreach (var i in _perturbedRegionsIndices) { + // rearranged to be more numerically stable weightFactor *= - -pathLength[i] * // dMua* factor - (Math.Exp(-perturbedOps[i].Mua * pathLength[i]) / Math.Exp(-_referenceOps[i].Mua * pathLength[i])); // mua pert - if (numberOfCollisions[i] > 0) - { - // the following is more numerically stable + -pathLength[i] * // dMua* factor Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * pathLength[i] / numberOfCollisions[i]), numberOfCollisions[i]); - } - else - { - weightFactor *= Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * pathLength[i]); - } } return weightFactor; } @@ -220,26 +215,18 @@ private double AbsorbDiscrete(IList numberOfCollisions, IList path { var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region + // NOTE: following code only works for single perturbed region because derivative of + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that foreach (var i in _perturbedRegionsIndices) { - if (numberOfCollisions[i] > 0) - { - weightFactor *= - -pathLength[i] * // dMua* factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - else // numberOfCollisions[i] in pert region is 0 - { - weightFactor *= - -pathLength[i] * // dMua* factor - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i]); - } + // rearranged to be more numerically stable + weightFactor *= + -pathLength[i] * // dMua* factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); } return weightFactor; } diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs index 65db23655..0e8497ddd 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs @@ -1,4 +1,4 @@ -using System; + using System; using System.Collections.Generic; using System.Linq; using System.Runtime.Serialization; @@ -137,10 +137,10 @@ public void Initialize(ITissue tissue, Random rng) TallyCount = 0; // if the data arrays are null, create them (only create second moment if TallySecondMoment is true) - Mean = Mean ?? new double[Rho.Count - 1]; - SecondMoment = SecondMoment ?? (TallySecondMoment ? new double[Rho.Count - 1] : null); + Mean ??= new double[Rho.Count - 1]; + SecondMoment ??= (TallySecondMoment ? new double[Rho.Count - 1] : null); - // intialize any other necessary class fields here + // initialize any other necessary class fields here _perturbedOps = PerturbedOps; _perturbedRegionsIndices = PerturbedRegionsIndices; _referenceOps = tissue.Regions.Select(r => r.RegionOP).ToList(); @@ -154,6 +154,7 @@ public void Initialize(ITissue tissue, Random rng) /// absorption weighting type protected void SetAbsorbAction(AbsorptionWeightingType awt) { + // AbsorptionWeightingType.Analog cannot have derivatives so not a case switch (awt) { // note: pMC is not applied to analog processing, @@ -193,30 +194,57 @@ private double AbsorbContinuous(IList numberOfCollisions, IList pa { var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region + // NOTE: following code only works for single perturbed region because derivative of + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that foreach (var i in _perturbedRegionsIndices) { - weightFactor *= - Math.Exp(-perturbedOps[i].Mus * pathLength[i]) / Math.Exp(-_referenceOps[i].Mus * pathLength[i]); // Mus pert if (numberOfCollisions[i] - 1 > 0) { + var dum = numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * + pathLength[i] / (numberOfCollisions[i] - 1)), + numberOfCollisions[i] - 1); + var dum2 = pathLength[i] * // dMus* 2nd factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp( + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); weightFactor *= - numberOfCollisions[i] * (1.0 / _referenceOps[i].Mus) * // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * - pathLength[i] / (numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1) // minus one here - - pathLength[i] * // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); // one here + numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * + pathLength[i] / (numberOfCollisions[i] - 1)), + numberOfCollisions[i] - 1) // (j-1) here + - pathLength[i] * // dMus* 2nd factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua ) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); // j here } - else + else // number of collisions in pert region is 1 { - weightFactor *= Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * pathLength[i]); + weightFactor *= + 1 / _referenceOps[i].Mus * // dMus* 1st factor + Math.Exp( + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mus) * + pathLength[i]) + - pathLength[i] * // dMus* 2nd factor + (_perturbedOps[i].Mus / _referenceOps[i].Mus) * + Math.Exp( + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mus) * + pathLength[i]); } } return weightFactor; @@ -226,41 +254,56 @@ private double AbsorbDiscrete(IList numberOfCollisions, IList path { var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region + // NOTE: following code only works for single perturbed region because derivative of + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that foreach (var i in _perturbedRegionsIndices) { if (numberOfCollisions[i] - 1 > 0) { + var dum = numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * + pathLength[i] / (numberOfCollisions[i] - 1)), + numberOfCollisions[i] - 1); + var dum2= pathLength[i] * // dMus* 2nd factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp( + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); weightFactor *= - numberOfCollisions[i]*(1.0/_referenceOps[i].Mus)* // dMus* 1st factor + numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor Math.Pow( - _perturbedOps[i].Mus/_referenceOps[i].Mus* - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - - _referenceOps[i].Mus)* + _perturbedOps[i].Mus/_referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * pathLength[i]/(numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1) // minus one here - - pathLength[i]* // dMus* 2nd factor + numberOfCollisions[i] - 1) // (j-1) here + - pathLength[i] * // dMus* 2nd factor Math.Pow( - _perturbedOps[i].Mus/_referenceOps[i].Mus* + _perturbedOps[i].Mus/_referenceOps[i].Mus * Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - - _referenceOps[i].Mus)* + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - + _referenceOps[i].Mua) * pathLength[i]/numberOfCollisions[i]), - numberOfCollisions[i]); // one here + numberOfCollisions[i]); // j here } else // number of collisions in pert region is 1 { weightFactor *= - numberOfCollisions[i] * (1.0 / _referenceOps[i].Mus) * // dMus* 1st factor + 1 / _referenceOps[i].Mus * // dMus* 1st factor Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mus) * pathLength[i]) - pathLength[i] * // dMus* 2nd factor (_perturbedOps[i].Mus / _referenceOps[i].Mus) * Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - + -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mus) * pathLength[i]); } From 661d23776ad0f028b637370cacf616d46ea8aae8 Mon Sep 17 00:00:00 2001 From: Janaka Ranasinghesagara Date: Thu, 2 May 2024 10:24:41 -0700 Subject: [PATCH 2/6] Applied simplified equations. --- .../Detectors/dMCdROfRhodMusDetector.cs | 111 +++--------------- 1 file changed, 14 insertions(+), 97 deletions(-) diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs index 0e8497ddd..a13d6ea41 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs @@ -198,54 +198,13 @@ private double AbsorbContinuous(IList numberOfCollisions, IList pa // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that foreach (var i in _perturbedRegionsIndices) { - if (numberOfCollisions[i] - 1 > 0) - { - var dum = numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i] / (numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1); - var dum2 = pathLength[i] * // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - weightFactor *= - numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i] / (numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1) // (j-1) here - - pathLength[i] * // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua ) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); // j here - } - else // number of collisions in pert region is 1 - { - weightFactor *= - 1 / _referenceOps[i].Mus * // dMus* 1st factor - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mus) * - pathLength[i]) - - pathLength[i] * // dMus* 2nd factor - (_perturbedOps[i].Mus / _referenceOps[i].Mus) * - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mus) * - pathLength[i]); - } + // rearranged to be more numerically stable + weightFactor *= + ((numberOfCollisions[i] / _perturbedOps[i].Mus) - pathLength[i]) * // dMus* factor + Math.Pow(_perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); } return weightFactor; } @@ -258,55 +217,13 @@ private double AbsorbDiscrete(IList numberOfCollisions, IList path // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that foreach (var i in _perturbedRegionsIndices) { - if (numberOfCollisions[i] - 1 > 0) - { - var dum = numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i] / (numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1); - var dum2= pathLength[i] * // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - weightFactor *= - numberOfCollisions[i] / _referenceOps[i].Mus * // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus/_referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i]/(numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1) // (j-1) here - - pathLength[i] * // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus/_referenceOps[i].Mus * - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mua) * - pathLength[i]/numberOfCollisions[i]), - numberOfCollisions[i]); // j here - } - else // number of collisions in pert region is 1 - { - weightFactor *= - 1 / _referenceOps[i].Mus * // dMus* 1st factor - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mus) * - pathLength[i]) - - pathLength[i] * // dMus* 2nd factor - (_perturbedOps[i].Mus / _referenceOps[i].Mus) * - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - - _referenceOps[i].Mus) * - pathLength[i]); - } + // rearranged to be more numerically stable + weightFactor *= + ((numberOfCollisions[i]/ _perturbedOps[i].Mus) -pathLength[i] )* // dMus* factor + Math.Pow(_perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); } return weightFactor; } From 6529314f02cf3f42b485a8a43bfb85a066654b37 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 2 May 2024 11:11:38 -0700 Subject: [PATCH 3/6] Merged my edits with prior push. --- .../Detectors/dMCdROfRhodMuaDetector.cs | 48 ++++++++----------- .../Detectors/dMCdROfRhodMusDetector.cs | 33 +++++-------- 2 files changed, 31 insertions(+), 50 deletions(-) diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs index 41dda30f9..e512a282e 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs @@ -193,42 +193,32 @@ public void Tally(Photon photon) private double AbsorbContinuous(IList numberOfCollisions, IList pathLength, IList perturbedOps) { - var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region because derivative of // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that - foreach (var i in _perturbedRegionsIndices) - { - // rearranged to be more numerically stable - weightFactor *= - -pathLength[i] * // dMua* factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - return weightFactor; + // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation + var i = _perturbedRegionsIndices[0]; + // rearranged to be more numerically stable + return -pathLength[i] * // dMua* factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); } private double AbsorbDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) { - var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region because derivative of - // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that - foreach (var i in _perturbedRegionsIndices) - { - // rearranged to be more numerically stable - weightFactor *= - -pathLength[i] * // dMua* factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - return weightFactor; + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that. + // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation + var i = _perturbedRegionsIndices[0]; + // rearranged to be more numerically stable + return -pathLength[i] * // dMua* factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); } /// diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs index a13d6ea41..ab819d41f 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs @@ -192,43 +192,34 @@ public void Tally(Photon photon) private double AbsorbContinuous(IList numberOfCollisions, IList pathLength, IList perturbedOps) { - var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region because derivative of // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that - foreach (var i in _perturbedRegionsIndices) - { - // rearranged to be more numerically stable - weightFactor *= - ((numberOfCollisions[i] / _perturbedOps[i].Mus) - pathLength[i]) * // dMus* factor + // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation + var i = _perturbedRegionsIndices[0]; + // rearranged to be more numerically stable + return (numberOfCollisions[i] / _perturbedOps[i].Mus - pathLength[i]) * Math.Pow(_perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - + _referenceOps[i].Mus - _referenceOps[i].Mua) * pathLength[i] / numberOfCollisions[i]), numberOfCollisions[i]); - } - return weightFactor; } private double AbsorbDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) { - var weightFactor = 1.0; - // NOTE: following code only works for single perturbed region because derivative of // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that - foreach (var i in _perturbedRegionsIndices) - { - // rearranged to be more numerically stable - weightFactor *= - ((numberOfCollisions[i]/ _perturbedOps[i].Mus) -pathLength[i] )* // dMus* factor + // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation + var i = _perturbedRegionsIndices[0]; + // rearranged to be more numerically stable + return (numberOfCollisions[i]/ _perturbedOps[i].Mus -pathLength[i] )* // dMus* factor Math.Pow(_perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - + _referenceOps[i].Mus - _referenceOps[i].Mua) * pathLength[i] / numberOfCollisions[i]), numberOfCollisions[i]); - } - return weightFactor; } - /// /// method to normalize detector results after numPhotons launched /// From 906b80b6dc71eb2df0f9578bdf4d6043827b23e7 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 2 May 2024 12:17:35 -0700 Subject: [PATCH 4/6] Cleaned up code per two new sonarcloud issues, i.e. two methods had identical code. That is because derivative Radon-Nikodym factor is same for DAW and CAW. --- .../Detectors/dMCdROfRhodMuaDetector.cs | 32 +++++++------------ .../Detectors/dMCdROfRhodMusDetector.cs | 32 +++++++------------ 2 files changed, 22 insertions(+), 42 deletions(-) diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs index e512a282e..807170383 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs @@ -158,17 +158,15 @@ protected void SetAbsorbAction(AbsorptionWeightingType awt) // AbsorptionWeightingType.Analog cannot have derivatives so not a case switch (awt) { - // note: pMC is not applied to analog processing, - // only DAW and CAW + // note: pMC is not applied to analog processing, only DAW and CAW case AbsorptionWeightingType.Continuous: - _absorbAction = AbsorbContinuous; - break; case AbsorptionWeightingType.Discrete: default: - _absorbAction = AbsorbDiscrete; + _absorbAction = AbsorbContinuousOrDiscrete; break; } } + /// /// method to tally to detector /// @@ -191,7 +189,14 @@ public void Tally(Photon photon) SecondMoment[ir] += photon.DP.Weight * weightFactor * photon.DP.Weight * weightFactor; } - private double AbsorbContinuous(IList numberOfCollisions, IList pathLength, IList perturbedOps) + /// + /// The following method works for both discrete or continuous absorption weighting + /// + /// photon number of collisions in perturbed region + /// photon path length in perturbed region + /// perturbed optical properties of perturbed region + /// derivative with respect to mus + private double AbsorbContinuousOrDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) { // NOTE: following code only works for single perturbed region because derivative of // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that @@ -206,21 +211,6 @@ private double AbsorbContinuous(IList numberOfCollisions, IList pa numberOfCollisions[i]); } - private double AbsorbDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) - { - // NOTE: following code only works for single perturbed region because derivative of - // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that. - // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation - var i = _perturbedRegionsIndices[0]; - // rearranged to be more numerically stable - return -pathLength[i] * // dMua* factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - /// /// method to normalize detector results after numPhotons launched /// diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs index ab819d41f..05e2404e9 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs @@ -157,17 +157,15 @@ protected void SetAbsorbAction(AbsorptionWeightingType awt) // AbsorptionWeightingType.Analog cannot have derivatives so not a case switch (awt) { - // note: pMC is not applied to analog processing, - // only DAW and CAW + // note: pMC is not applied to analog processing, only DAW and CAW case AbsorptionWeightingType.Continuous: - _absorbAction = AbsorbContinuous; - break; case AbsorptionWeightingType.Discrete: default: - _absorbAction = AbsorbDiscrete; + _absorbAction = AbsorbContinuousOrDiscrete; break; } } + /// /// method to tally to detector /// @@ -190,7 +188,14 @@ public void Tally(Photon photon) SecondMoment[ir] += photon.DP.Weight * weightFactor * photon.DP.Weight * weightFactor; } - private double AbsorbContinuous(IList numberOfCollisions, IList pathLength, IList perturbedOps) + /// + /// The following method works for both discrete or continuous absorption weighting + /// + /// photon number of collisions in perturbed region + /// photon path length in perturbed region + /// perturbed optical properties of perturbed region + /// derivative with respect to mua + private double AbsorbContinuousOrDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) { // NOTE: following code only works for single perturbed region because derivative of // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that @@ -205,21 +210,6 @@ private double AbsorbContinuous(IList numberOfCollisions, IList pa numberOfCollisions[i]); } - private double AbsorbDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) - { - // NOTE: following code only works for single perturbed region because derivative of - // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that - // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation - var i = _perturbedRegionsIndices[0]; - // rearranged to be more numerically stable - return (numberOfCollisions[i]/ _perturbedOps[i].Mus -pathLength[i] )* // dMus* factor - Math.Pow(_perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - /// /// method to normalize detector results after numPhotons launched /// From 042d94472d817acd9f9964db1fe13b941c3d0e50 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 2 May 2024 15:58:43 -0700 Subject: [PATCH 5/6] Cleaned up issue with 2 issues with switch statement on the absorption weighting type: 1) added Analog option so switch encompasses all enum types, and 2) removal of the empty case clause. --- .../Detectors/dMCdROfRhodMuaDetector.cs | 15 +++++++-------- .../Detectors/dMCdROfRhodMusDetector.cs | 15 +++++++-------- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs index 807170383..04a99d148 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs @@ -155,16 +155,15 @@ public void Initialize(ITissue tissue, Random rng) /// absorption weighting type protected void SetAbsorbAction(AbsorptionWeightingType awt) { - // AbsorptionWeightingType.Analog cannot have derivatives so not a case - switch (awt) + // AbsorptionWeightingType.Analog cannot have derivatives so exception + _absorbAction = awt switch { // note: pMC is not applied to analog processing, only DAW and CAW - case AbsorptionWeightingType.Continuous: - case AbsorptionWeightingType.Discrete: - default: - _absorbAction = AbsorbContinuousOrDiscrete; - break; - } + AbsorptionWeightingType.Continuous => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Discrete => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Analog => throw new ArgumentException(@"Analog is not allowed with this detector", nameof(awt)), + _ => throw new ArgumentOutOfRangeException(typeof(AbsorptionWeightingType).ToString()) + }; } /// diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs index 05e2404e9..8f14c7cd0 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs @@ -154,16 +154,15 @@ public void Initialize(ITissue tissue, Random rng) /// absorption weighting type protected void SetAbsorbAction(AbsorptionWeightingType awt) { - // AbsorptionWeightingType.Analog cannot have derivatives so not a case - switch (awt) + // AbsorptionWeightingType.Analog cannot have derivatives so exception + _absorbAction = awt switch { // note: pMC is not applied to analog processing, only DAW and CAW - case AbsorptionWeightingType.Continuous: - case AbsorptionWeightingType.Discrete: - default: - _absorbAction = AbsorbContinuousOrDiscrete; - break; - } + AbsorptionWeightingType.Continuous => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Discrete => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Analog => throw new ArgumentException(@"Analog is not allowed with this detector", nameof(awt)), + _ => throw new ArgumentOutOfRangeException(typeof(AbsorptionWeightingType).ToString()) + }; } /// From fdd6f3a58bbc719530f75f9ecfd78f6369d1cc85 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 2 May 2024 16:35:10 -0700 Subject: [PATCH 6/6] Added unit tests to test specification of AbsorptionWeightingType.Analog and the resulting exception. --- .../pMCdMCDAWOneLayerDetectorsTests.cs | 53 +++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs index f1c1e33eb..5b1ccf058 100644 --- a/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs +++ b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs @@ -491,6 +491,59 @@ public void Validate_pMC_DAW_ROfFx_nonzero_perturbation_one_layer_tissue() Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); } + [Test] + public void Test_Analog_absorption_weighting_type_throws_argument_exception() + { + var test1 = new DMuaDetectorTest(); + Assert.Throws(() => + { + try + { + test1.SetAbsorbAction(AbsorptionWeightingType.Analog); + } + catch (Exception e) + { + Assert.AreEqual("Analog is not allowed with this detector (Parameter 'awt')", e.Message); + throw; + } + }); + var test2 = new DMusDetectorTest(); + Assert.Throws(() => + { + try + { + test2.SetAbsorbAction(AbsorptionWeightingType.Analog); + } + catch (Exception e) + { + Assert.AreEqual("Analog is not allowed with this detector (Parameter 'awt')", e.Message); + throw; + } + }); + } + + /// + /// Expose protected method in a new class that inherits the class under test + /// + public class DMuaDetectorTest : dMCdROfRhodMuaDetector + { + public new void SetAbsorbAction(AbsorptionWeightingType awt) + { + base.SetAbsorbAction(awt); + } + } + + /// + /// Expose protected method in a new class that inherits the class under test + /// + public class DMusDetectorTest : dMCdROfRhodMusDetector + { + public new void SetAbsorbAction(AbsorptionWeightingType awt) + { + base.SetAbsorbAction(awt); + } + } + } }