Commit 0ccc8ee8 authored by Danny Hindson's avatar Danny Hindson
Browse files

Incorporate review comments

Main change is change IObject::generatePointInObject to return boost::optional<V3D>
instead of using an out parameter combined with a boolean return value
parent e2af417f
......@@ -38,15 +38,16 @@ public:
size_t maxScatterPtAttempts);
std::tuple<double, double> calculate(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &finalPos,
double lambdaBefore, double lambdaAfter,
std::string &debugString) const;
double lambdaBefore, double lambdaAfter);
std::string getDebugString() const { return debugString; };
private:
const IBeamProfile &m_beamProfile;
const MCInteractionVolume m_scatterVol;
MCInteractionVolume m_scatterVol;
const size_t m_nevents;
const size_t m_maxScatterAttempts;
const double m_error;
std::string debugString;
};
} // namespace Algorithms
......
......@@ -9,6 +9,7 @@
#include "MantidAlgorithms/DllConfig.h"
#include "MantidGeometry/Objects/BoundingBox.h"
#include <boost/optional.hpp>
namespace Mantid {
namespace API {
......@@ -45,13 +46,18 @@ public:
double calculateAbsorption(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &startPos,
const Kernel::V3D &endPos, double lambdaBefore,
double lambdaAfter) const;
double lambdaAfter);
std::string generateScatterPointStats() const;
Kernel::V3D generatePoint(Kernel::PseudoRandomNumberGenerator &rng) const;
Kernel::V3D generatePoint(Kernel::PseudoRandomNumberGenerator &rng);
private:
mutable int m_sampleScatterPoints = 0;
mutable std::vector<int> m_envScatterPoints;
int getComponentIndex(Kernel::PseudoRandomNumberGenerator &rng);
boost::optional<Kernel::V3D>
generatePointInObjectByIndex(int componentIndex,
Kernel::PseudoRandomNumberGenerator &rng);
void UpdateScatterPointCounts(int componentIndex);
int m_sampleScatterPoints = 0;
std::vector<int> m_envScatterPoints;
const boost::shared_ptr<Geometry::IObject> m_sample;
const Geometry::SampleEnvironment *m_env;
const Geometry::BoundingBox m_activeRegion;
......
......@@ -276,10 +276,10 @@ MatrixWorkspace_uptr MonteCarloAbsorption::doSimulation(
} else {
// elastic case already initialized
}
std::string debugString;
std::tie(outY[j], std::ignore) =
strategy.calculate(rng, detPos, lambdaIn, lambdaOut, debugString);
g_log.debug(debugString);
strategy.calculate(rng, detPos, lambdaIn, lambdaOut);
g_log.debug(strategy.getDebugString());
// Ensure we have the last point for the interpolation
if (lambdaStepSize > 1 && j + lambdaStepSize >= nbins && j + 1 != nbins) {
......
......@@ -43,12 +43,12 @@ MCAbsorptionStrategy::MCAbsorptionStrategy(const IBeamProfile &beamProfile,
* where it is detected
* @param lambdaBefore Wavelength, in \f$\\A^-1\f$, before scattering
* @param lambdaAfter Wavelength, in \f$\\A^-1\f$, after scattering
* @param debugString string to return any debug messages
* @return A tuple of the <correction factor, associated error>.
*/
std::tuple<double, double> MCAbsorptionStrategy::calculate(
Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &finalPos,
double lambdaBefore, double lambdaAfter, std::string &debugString) const {
std::tuple<double, double>
MCAbsorptionStrategy::calculate(Kernel::PseudoRandomNumberGenerator &rng,
const Kernel::V3D &finalPos,
double lambdaBefore, double lambdaAfter) {
const auto scatterBounds = m_scatterVol.getBoundingBox();
double factor(0.0);
for (size_t i = 0; i < m_nevents; ++i) {
......
......@@ -40,6 +40,7 @@ MCInteractionVolume::MCInteractionVolume(
}
try {
m_env = &sample.getEnvironment();
assert(m_env);
if (m_env->nelements() == 0) {
throw std::invalid_argument(
"MCInteractionVolume() - Sample enviroment has zero components.");
......@@ -61,6 +62,58 @@ const Geometry::BoundingBox &MCInteractionVolume::getBoundingBox() const {
return m_sample->getBoundingBox();
}
/**
* Randomly select a component across the sample/environment
* @param rng A reference to a PseudoRandomNumberGenerator where
* nextValue should return a flat random number between 0.0 & 1.0
* @return The randomly selected component index
*/
int MCInteractionVolume::getComponentIndex(
Kernel::PseudoRandomNumberGenerator &rng) {
int componentIndex;
// the sample has componentIndex -1, env components are number 0 upwards
if (m_env) {
componentIndex = rng.nextInt(0, static_cast<int>(m_env->nelements())) - 1;
} else {
componentIndex = -1;
}
return componentIndex;
}
/**
* Generate a point in an object identified by an index
* @param componentIndex Index of the sample/environment component where
* the sample is -1
* @param rng A reference to a PseudoRandomNumberGenerator where
* nextValue should return a flat random number between 0.0 & 1.0
* @return The generated point
*/
boost::optional<Kernel::V3D> MCInteractionVolume::generatePointInObjectByIndex(
int componentIndex, Kernel::PseudoRandomNumberGenerator &rng) {
boost::optional<Kernel::V3D> pointGenerated;
if (componentIndex == -1) {
pointGenerated = m_sample->generatePointInObject(rng, m_activeRegion, 1);
} else {
pointGenerated = m_env->getComponent(componentIndex)
.generatePointInObject(rng, m_activeRegion, 1);
}
return pointGenerated;
}
/**
* Update the scatter point counts
* @param componentIndex Index of the sample/environment component where
* the sample is -1
*/
void MCInteractionVolume::UpdateScatterPointCounts(int componentIndex) {
if (componentIndex == -1) {
m_sampleScatterPoints++;
} else {
m_envScatterPoints[componentIndex]++;
}
}
/**
* Generate point randomly across one of the components of the environment
* including the sample itself in the selection. The method first selects
......@@ -70,33 +123,15 @@ const Geometry::BoundingBox &MCInteractionVolume::getBoundingBox() const {
* nextValue should return a flat random number between 0.0 & 1.0
* @return The generated point
*/
Kernel::V3D MCInteractionVolume::generatePoint(
Kernel::PseudoRandomNumberGenerator &rng) const {
Kernel::V3D
MCInteractionVolume::generatePoint(Kernel::PseudoRandomNumberGenerator &rng) {
for (size_t i = 0; i < m_maxScatterAttempts; i++) {
Kernel::V3D point;
int componentIndex;
if (m_env) {
componentIndex = rng.nextInt(0, static_cast<int>(m_env->nelements())) - 1;
} else {
componentIndex = -1;
}
bool pointGenerated;
if (componentIndex == -1) {
pointGenerated =
m_sample->generatePointInObject(rng, m_activeRegion, 1, point);
} else {
pointGenerated =
m_env->getComponent(componentIndex)
.generatePointInObject(rng, m_activeRegion, 1, point);
}
int componentIndex = getComponentIndex(rng);
boost::optional<Kernel::V3D> pointGenerated =
generatePointInObjectByIndex(componentIndex, rng);
if (pointGenerated) {
if (componentIndex == -1) {
m_sampleScatterPoints++;
} else {
m_envScatterPoints[componentIndex]++;
}
return point;
UpdateScatterPointCounts(componentIndex);
return *pointGenerated;
}
}
throw std::runtime_error("MCInteractionVolume::generatePoint() - Unable to "
......@@ -119,7 +154,7 @@ Kernel::V3D MCInteractionVolume::generatePoint(
*/
double MCInteractionVolume::calculateAbsorption(
Kernel::PseudoRandomNumberGenerator &rng, const Kernel::V3D &startPos,
const Kernel::V3D &endPos, double lambdaBefore, double lambdaAfter) const {
const Kernel::V3D &endPos, double lambdaBefore, double lambdaAfter) {
// Generate scatter point. If there is an environment present then
// first select whether the scattering occurs on the sample or the
// environment. The attenuation for the path leading to the scatter point
......@@ -174,10 +209,10 @@ std::string MCInteractionVolume::generateScatterPointStats() const {
scatterPointSummary << "Scatter point counts:" << std::endl;
int totalScatterPoints = m_sampleScatterPoints;
for (std::vector<int>::size_type i = 0; i < m_envScatterPoints.size(); i++) {
totalScatterPoints += m_envScatterPoints[i];
}
int totalScatterPoints =
std::accumulate(m_envScatterPoints.begin(), m_envScatterPoints.end(),
m_sampleScatterPoints);
scatterPointSummary << "Total scatter points: " << totalScatterPoints
<< std::endl;
......
......@@ -57,9 +57,8 @@ public:
const double lambdaBefore(2.5), lambdaAfter(3.5);
double factor(0.0), error(0.0);
std::string debugString;
std::tie(factor, error) =
mcabsorb.calculate(rng, endPos, lambdaBefore, lambdaAfter, debugString);
mcabsorb.calculate(rng, endPos, lambdaBefore, lambdaAfter);
TS_ASSERT_DELTA(0.0043828472, factor, 1e-08);
TS_ASSERT_DELTA(1.0 / std::sqrt(nevents), error, 1e-08);
}
......@@ -86,10 +85,8 @@ public:
EXPECT_CALL(rng, nextValue()).WillRepeatedly(Return(0.5));
const double lambdaBefore(2.5), lambdaAfter(3.5);
const V3D endPos(0.7, 0.7, 1.4);
std::string debugString;
TS_ASSERT_THROWS(
mcabs.calculate(rng, endPos, lambdaBefore, lambdaAfter, debugString),
const std::runtime_error &)
TS_ASSERT_THROWS(mcabs.calculate(rng, endPos, lambdaBefore, lambdaAfter),
const std::runtime_error &)
}
private:
......
......@@ -16,6 +16,9 @@
#include <gmock/gmock.h>
using Mantid::Algorithms::MCInteractionVolume;
using namespace ::testing;
using namespace Mantid::Kernel;
using namespace MonteCarloTesting;
class MCInteractionVolumeTest : public CxxTest::TestSuite {
public:
......@@ -30,7 +33,6 @@ public:
// Success cases
//----------------------------------------------------------------------------
void test_Bounding_Volume_Matches_Sample() {
using namespace MonteCarloTesting;
auto sample = createTestSample(TestSampleType::SolidSphere);
const auto sampleBox = sample.getShape().getBoundingBox();
MCInteractionVolume interactor(sample, sampleBox);
......@@ -42,8 +44,6 @@ public:
void test_Absorption_In_Solid_Sample_Gives_Expected_Answer() {
using Mantid::Kernel::V3D;
using namespace MonteCarloTesting;
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), endPos(0.7, 0.7, 1.4);
......@@ -62,8 +62,6 @@ public:
void test_Absorption_In_Sample_With_Hole_Container_Scatter_In_All_Segments() {
using Mantid::Kernel::V3D;
using namespace MonteCarloTesting;
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), endPos(2.0, 0.0, 0.0);
......@@ -95,8 +93,6 @@ public:
void
test_Absorption_In_Sample_And_Environment_Container_Scatter_In_All_Segments() {
using Mantid::Kernel::V3D;
using namespace MonteCarloTesting;
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), endPos(2.0, 0.0, 0.0);
......@@ -138,9 +134,7 @@ public:
// Failure cases
//----------------------------------------------------------------------------
void test_Construction_With_Invalid_Sample_Shape_Throws_Error() {
using Mantid::API::Sample;
Sample sample;
Mantid::API::Sample sample;
// nothing
TS_ASSERT_THROWS(
MCInteractionVolume mcv(sample, sample.getShape().getBoundingBox()),
......@@ -152,9 +146,6 @@ public:
}
void test_Throws_If_Point_Cannot_Be_Generated() {
using namespace Mantid::Kernel;
using namespace MonteCarloTesting;
using namespace ::testing;
// Testing inputs
const V3D startPos(-2.0, 0.0, 0.0), endPos(2.0, 0.0, 0.0);
......@@ -172,15 +163,11 @@ public:
}
void testGeneratePointConsidersAllComponents() {
using namespace ::testing;
using namespace Mantid::Kernel;
using namespace MonteCarloTesting;
auto kit = createTestKit();
size_t maxAttempts(1);
using Mantid::API::Sample;
Sample sample;
Mantid::API::Sample sample;
sample.setShape(ComponentCreationHelper::createSphere(1));
sample.setEnvironment(
std::make_unique<Mantid::Geometry::SampleEnvironment>(*kit));
......@@ -236,9 +223,6 @@ public:
}
void testGeneratePointRespectsActiveRegion() {
using namespace MonteCarloTesting;
using namespace ::testing;
using namespace Mantid::Kernel;
auto kit = createTestKit();
size_t maxAttempts(1);
......
......@@ -79,15 +79,16 @@ public:
int getPointInObject(Kernel::V3D &point) const override {
return m_shape->getPointInObject(point);
}
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t i,
Kernel::V3D &point) const override {
return m_shape->generatePointInObject(rng, i, point);
}
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion, const size_t i,
Kernel::V3D &point) const override {
return m_shape->generatePointInObject(rng, activeRegion, i, point);
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t i) const override {
return m_shape->generatePointInObject(rng, i);
}
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t i) const override {
return m_shape->generatePointInObject(rng, activeRegion, i);
}
detail::ShapeInfo::GeometryShape shape() const override {
......
......@@ -165,11 +165,13 @@ public:
int getPointInObject(Kernel::V3D &point) const override;
/// Select a random point within the object
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t, Kernel::V3D &point) const override;
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion, const size_t,
Kernel::V3D &point) const override;
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t) const override;
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t) const override;
// Rendering member functions
void draw() const override;
......
......@@ -9,6 +9,7 @@
#include "MantidGeometry/DllConfig.h"
#include "MantidGeometry/Rendering/ShapeInfo.h"
#include <boost/optional.hpp>
#include <boost/shared_ptr.hpp>
#include <map>
#include <vector>
......@@ -70,13 +71,13 @@ public:
virtual int getPointInObject(Kernel::V3D &point) const = 0;
virtual bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t,
Kernel::V3D &point) const = 0;
virtual bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t,
Kernel::V3D &point) const = 0;
virtual boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t) const = 0;
virtual boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t) const = 0;
virtual detail::ShapeInfo::GeometryShape shape() const = 0;
virtual const detail::ShapeInfo &shapeInfo() const = 0;
......
......@@ -110,11 +110,13 @@ public:
int getPointInObject(Kernel::V3D &point) const override;
/// Select a random point within the object
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t, Kernel::V3D &point) const override;
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion, const size_t,
Kernel::V3D &point) const override;
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t) const override;
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t) const override;
// Rendering member functions
void draw() const override;
......
......@@ -65,11 +65,13 @@ public:
double &ymin, double &zmin) const override;
int getPointInObject(Kernel::V3D &point) const override;
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t, Kernel::V3D &point) const override;
bool generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion, const size_t,
Kernel::V3D &point) const override;
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t) const override;
boost::optional<Kernel::V3D>
generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t) const override;
detail::ShapeInfo::GeometryShape shape() const override;
const detail::ShapeInfo &shapeInfo() const override;
void GetObjectGeom(detail::ShapeInfo::GeometryShape &type,
......
......@@ -1982,12 +1982,12 @@ int CSGObject::getPointInObject(Kernel::V3D &point) const {
* @param rng A reference to a PseudoRandomNumberGenerator where
* nextValue should return a flat random number between 0.0 & 1.0
* @param maxAttempts The maximum number of attempts at generating a point
* @param point The generated point
* @return whether a point was generated in the object or not
*/
bool CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
const size_t maxAttempts,
V3D &point) const {
boost::optional<Kernel::V3D>
CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
const size_t maxAttempts) const {
boost::optional<V3D> point{boost::none};
// If the shape fills its bounding box well enough then the most efficient
// way to get the point is just brute force. We'll try that first with
// just a few attempts.
......@@ -2000,7 +2000,7 @@ bool CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
boost::optional<V3D> maybePoint{
RandomPoint::inGenericShape(*this, rng, bruteForceAttempts)};
if (maybePoint) {
point = *maybePoint;
point = maybePoint;
} else {
switch (shape()) {
case detail::ShapeInfo::GeometryShape::CUBOID:
......@@ -2018,15 +2018,10 @@ bool CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
default:
maybePoint = RandomPoint::inGenericShape(
*this, rng, maxAttempts - bruteForceAttempts);
if (!maybePoint) {
return false;
} else {
point = *maybePoint;
break;
}
point = maybePoint;
}
}
return true;
return point;
}
/**
......@@ -2037,14 +2032,13 @@ bool CSGObject::generatePointInObject(PseudoRandomNumberGenerator &rng,
* @param activeRegion Restrict point generation to this sub-region of the
* object
* @param maxAttempts The maximum number of attempts at generating a point
* @param point The generated point
* @return whether a point was generated in the object or not
*/
bool CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t maxAttempts,
V3D &point) const {
boost::optional<V3D> pointGenerated{boost::none};
boost::optional<Kernel::V3D>
CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const BoundingBox &activeRegion,
const size_t maxAttempts) const {
boost::optional<V3D> point{boost::none};
// We'll first try brute force. If the shape fills its bounding box
// well enough, this should be the fastest method.
// Increasing the brute force attemps speeds up the shapes that fill
......@@ -2053,9 +2047,8 @@ bool CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
// shape, its dimension and orientation.
const size_t bruteForceAttempts{
std::min(static_cast<size_t>(5), maxAttempts)};
pointGenerated =
RandomPoint::bounded(*this, rng, activeRegion, bruteForceAttempts);
if (!pointGenerated) {
point = RandomPoint::bounded(*this, rng, activeRegion, bruteForceAttempts);
if (!point) {
detail::ShapeInfo::GeometryShape shape;
std::vector<Kernel::V3D> shapeVectors;
double radius;
......@@ -2064,37 +2057,32 @@ bool CSGObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
GetObjectGeom(shape, shapeVectors, innerRadius, radius, height);
switch (shape) {
case detail::ShapeInfo::GeometryShape::CUBOID:
pointGenerated = RandomPoint::bounded<RandomPoint::inCuboid>(
point = RandomPoint::bounded<RandomPoint::inCuboid>(
m_handler->shapeInfo(), rng, activeRegion,
maxAttempts - bruteForceAttempts);
break;
case detail::ShapeInfo::GeometryShape::CYLINDER:
pointGenerated = RandomPoint::bounded<RandomPoint::inCylinder>(
point = RandomPoint::bounded<RandomPoint::inCylinder>(
m_handler->shapeInfo(), rng, activeRegion,
maxAttempts - bruteForceAttempts);
break;
case detail::ShapeInfo::GeometryShape::HOLLOWCYLINDER:
pointGenerated = RandomPoint::bounded<RandomPoint::inHollowCylinder>(
point = RandomPoint::bounded<RandomPoint::inHollowCylinder>(
m_handler->shapeInfo(), rng, activeRegion,
maxAttempts - bruteForceAttempts);
break;
case detail::ShapeInfo::GeometryShape::SPHERE:
pointGenerated = RandomPoint::bounded<RandomPoint::inSphere>(
point = RandomPoint::bounded<RandomPoint::inSphere>(
m_handler->shapeInfo(), rng, activeRegion,
maxAttempts - bruteForceAttempts);
break;
default:
pointGenerated = RandomPoint::bounded(*this, rng, activeRegion,
maxAttempts - bruteForceAttempts);
point = RandomPoint::bounded(*this, rng, activeRegion,
maxAttempts - bruteForceAttempts);
break;
}
}
if (!pointGenerated) {
return false;
} else {
point = *pointGenerated;
return true;
}
return point;
}
/**
......
......@@ -383,18 +383,17 @@ int MeshObject::getPointInObject(Kernel::V3D &point) const {
* @param rng A reference to a PseudoRandomNumberGenerator where
* nextValue should return a flat random number between 0.0 & 1.0
* @param maxAttempts The maximum number of attempts at generating a point
* @param point The newly generated point
* @return Whether a point was generated or not
* @return The generated point
*/
bool MeshObject::generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng,
const size_t maxAttempts,
Kernel::V3D &point) const {
boost::optional<Kernel::V3D>