Commit fad82185 authored by Antti Soininen's avatar Antti Soininen
Browse files

Fix cuboid volume calculation and rework the Monte Carlo method.

The Monte Carlo method now has two methods: one runs the simulation a
fixed number of iterations, the other manages the first and checks for
convergence.

Re #19154
parent c268ee4f
...@@ -217,6 +217,8 @@ private: ...@@ -217,6 +217,8 @@ private:
/// Returns the volume. /// Returns the volume.
double monteCarloVolume() const; double monteCarloVolume() const;
/// Returns the volume.
double singleShotMonteCarloVolume(const int shotSize, const size_t seed) const;
/// Top rule [ Geometric scope of object] /// Top rule [ Geometric scope of object]
std::unique_ptr<Rule> TopRule; std::unique_ptr<Rule> TopRule;
......
#include "MantidGeometry/Objects/Object.h" #include "MantidGeometry/Objects/Object.h"
#include "MantidGeometry/Objects/Rules.h" #include "MantidGeometry/Objects/Rules.h"
#include "MantidGeometry/Objects/Track.h" #include "MantidGeometry/Objects/Track.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/PseudoRandomNumberGenerator.h"
#include "MantidKernel/Strings.h"
#include "MantidGeometry/Surfaces/Cone.h"
#include "MantidGeometry/Surfaces/Cylinder.h"
#include "MantidGeometry/Surfaces/LineIntersectVisit.h"
#include "MantidGeometry/Surfaces/Surface.h"
#include "MantidGeometry/Rendering/CacheGeometryHandler.h" #include "MantidGeometry/Rendering/CacheGeometryHandler.h"
#include "MantidGeometry/Rendering/GeometryHandler.h" #include "MantidGeometry/Rendering/GeometryHandler.h"
#include "MantidGeometry/Rendering/GluGeometryHandler.h" #include "MantidGeometry/Rendering/GluGeometryHandler.h"
#include "MantidGeometry/Rendering/vtkGeometryCacheReader.h" #include "MantidGeometry/Rendering/vtkGeometryCacheReader.h"
#include "MantidGeometry/Rendering/vtkGeometryCacheWriter.h" #include "MantidGeometry/Rendering/vtkGeometryCacheWriter.h"
#include "MantidGeometry/Surfaces/Cone.h"
#include "MantidGeometry/Surfaces/Cylinder.h"
#include "MantidGeometry/Surfaces/LineIntersectVisit.h"
#include "MantidGeometry/Surfaces/Surface.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/make_unique.h"
#include "MantidKernel/Material.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/PseudoRandomNumberGenerator.h"
#include "MantidKernel/Quat.h" #include "MantidKernel/Quat.h"
#include "MantidKernel/RegexStrings.h" #include "MantidKernel/RegexStrings.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/Tolerance.h" #include "MantidKernel/Tolerance.h"
#include "MantidKernel/make_unique.h"
#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/error_of_mean.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/make_shared.hpp> #include <boost/make_shared.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/ranlux.hpp>
#include <boost/random/uniform_01.hpp>
#include <array> #include <array>
#include <deque> #include <deque>
...@@ -1490,20 +1496,33 @@ double Object::volume() const { ...@@ -1490,20 +1496,33 @@ double Object::volume() const {
static_cast<GluGeometryHandler::GeometryType>(type); static_cast<GluGeometryHandler::GeometryType>(type);
switch (gluType) { switch (gluType) {
case GluGeometryHandler::GeometryType::CUBOID: { case GluGeometryHandler::GeometryType::CUBOID: {
// Here, the volume is calculated by the triangular method.
// We use one of the vertices (vectors[0]) as the reference
// point.
double volume = 0.0; double volume = 0.0;
const Kernel::V3D vertex12 = vectors[1] + vectors[2]; // Vertices. Name codes follow flb = front-left-bottom etc.
const Kernel::V3D vertex13 = vectors[1] + vectors[3]; const Kernel::V3D &flb = vectors[0];
const Kernel::V3D vertex23 = vectors[2] + vectors[3]; const Kernel::V3D &flt = vectors[1];
const Kernel::V3D vertex123 = vectors[1] + vertex23; const Kernel::V3D &frb = vectors[3];
volume += vertex23.scalar_prod(vectors[3].cross_prod(vectors[2])); const Kernel::V3D frt = frb + flt - flb;
volume += vertex23.scalar_prod(vectors[2].cross_prod(vertex123)); const Kernel::V3D &blb = vectors[2];
volume += vertex23.scalar_prod(vertex123.cross_prod(vectors[3])); const Kernel::V3D blt = blb + flt - flb;
volume += vertex12.scalar_prod(vectors[2].cross_prod(vectors[1])); const Kernel::V3D brb = blb + frb - flb;
volume += vertex12.scalar_prod(vectors[1].cross_prod(vertex13)); const Kernel::V3D brt = frt + blb - flb;
volume += vertex12.scalar_prod(vertex123.cross_prod(vectors[2])); // Normals point out, follow right-handed rule when
volume += vectors[1].scalar_prod(vectors[3].cross_prod(vertex13)); // defining the triangle faces.
volume += vertex13.scalar_prod(vectors[3].cross_prod(vertex123)); volume += flb.scalar_prod(flt.cross_prod(blb));
volume += vertex13.scalar_prod(vertex123.cross_prod(vertex23)); volume += blb.scalar_prod(flt.cross_prod(blt));
volume += flb.scalar_prod(frb.cross_prod(flt));
volume += frb.scalar_prod(frt.cross_prod(flt));
volume += flb.scalar_prod(blb.cross_prod(frb));
volume += blb.scalar_prod(brb.cross_prod(frb));
volume += frb.scalar_prod(brb.cross_prod(frt));
volume += brb.scalar_prod(brt.cross_prod(frt));
volume += flt.scalar_prod(frt.cross_prod(blt));
volume += frt.scalar_prod(brt.cross_prod(blt));
volume += blt.scalar_prod(brt.cross_prod(blb));
volume += brt.scalar_prod(brb.cross_prod(blb));
return volume / 6; return volume / 6;
} }
case GluGeometryHandler::GeometryType::SPHERE: case GluGeometryHandler::GeometryType::SPHERE:
...@@ -1517,44 +1536,72 @@ double Object::volume() const { ...@@ -1517,44 +1536,72 @@ double Object::volume() const {
} }
/** /**
* The volume is calculated using the Monte Carlo method. * Calculates the volume of this object by the Monte Carlo method.
* @returns The volume of this object. * This method manages singleShotMonteCarloVolume() and uses the
* variance as the convergence criteria.
* @returns The simulated volume of this object.
*/ */
double Object::monteCarloVolume() const { double Object::monteCarloVolume() const {
using namespace boost::accumulators;
const int singleShotIterations = 10000;
accumulator_set<double, features<tag::mean, tag::error_of<tag::mean>>> accumulate;
// For seeding the single shot runs.
boost::random::ranlux48 rnEngine;
// Warm up statistics.
for (int i = 0; i < 10; ++i) {
const auto seed = rnEngine();
const double volume = singleShotMonteCarloVolume(singleShotIterations, seed);
accumulate(volume);
}
const double relativeErrorTolerance = 1e-3;
double currentMean;
double currentError;
do {
const auto seed = rnEngine();
const double volume = singleShotMonteCarloVolume(singleShotIterations, seed);
accumulate(volume);
currentMean = mean(accumulate);
currentError = error_of<tag::mean>(accumulate);
if (std::isnan(currentError)) {
currentError = 0;
}
} while (currentError / currentMean > relativeErrorTolerance);
return currentMean;
}
/**
* Calculates the volume using the Monte Carlo method.
* @param shotSize Number of iterations.
* @param seed A number to seed the random number generator.
* @returns The simulated volume of this object.
*/
double Object::singleShotMonteCarloVolume(const int shotSize, const size_t seed) const {
const auto &boundingBox = getBoundingBox(); const auto &boundingBox = getBoundingBox();
if (boundingBox.isNull()) { if (boundingBox.isNull()) {
return 0; return 0;
} }
int total = 0;
int totalHits = 0; int totalHits = 0;
// The maximum number of iterations is chosen by the
// Stetson-Harrison method. This means that the following
// number may or may not work in every scenario.
const int maxIterations = 1000000;
const double boundingDx = boundingBox.xMax() - boundingBox.xMin(); const double boundingDx = boundingBox.xMax() - boundingBox.xMin();
const double boundingDy = boundingBox.yMax() - boundingBox.yMin(); const double boundingDy = boundingBox.yMax() - boundingBox.yMin();
const double boundingDz = boundingBox.yMax() - boundingBox.yMin(); const double boundingDz = boundingBox.zMax() - boundingBox.zMin();
PARALLEL { PARALLEL {
const auto threadCount = PARALLEL_NUMBER_OF_THREADS; const auto threadCount = PARALLEL_NUMBER_OF_THREADS;
const auto currentThreadNum = PARALLEL_THREAD_NUMBER; const auto currentThreadNum = PARALLEL_THREAD_NUMBER;
size_t blocksize = maxIterations / threadCount; size_t blocksize = shotSize / threadCount;
if (currentThreadNum == threadCount - 1) { if (currentThreadNum == threadCount - 1) {
// Last thread may do some extra work. // Last thread may do some extra work.
blocksize = maxIterations - (threadCount - 1) * blocksize; blocksize = shotSize - (threadCount - 1) * blocksize;
} }
// The Mersenne twister randon number engine has poor performance boost::random::mt19937 rnEngine(static_cast<boost::random::mt19937::result_type>(seed));
// with discard(). Use more suitable engine for parallelization. // All threads init their engine with the same seed.
std::ranlux48 rnEngine;
// All threads init their engine with the same (defualt) seed.
// We discard the random numbers used by the other threads. // We discard the random numbers used by the other threads.
// This ensures reproducible results independent of the number // This ensures reproducible results independent of the number
// of threads. // of threads.
// We need three random numbers for each iteration. // We need three random numbers for each iteration.
rnEngine.discard(currentThreadNum * 3 * blocksize); rnEngine.discard(currentThreadNum * 3 * blocksize);
std::uniform_real_distribution<double> rnDistribution(0, 1); boost::random::uniform_01<double> rnDistribution;
int hits = 0; int hits = 0;
int localTotal = 0; for (int i = 0; i < static_cast<int>(blocksize); ++i) {
for (; localTotal < static_cast<int>(blocksize); ++localTotal) {
double rnd = rnDistribution(rnEngine); double rnd = rnDistribution(rnEngine);
const double x = boundingBox.xMin() + rnd * boundingDx; const double x = boundingBox.xMin() + rnd * boundingDx;
rnd = rnDistribution(rnEngine); rnd = rnDistribution(rnEngine);
...@@ -1566,13 +1613,11 @@ double Object::monteCarloVolume() const { ...@@ -1566,13 +1613,11 @@ double Object::monteCarloVolume() const {
} }
} }
// Collect results. // Collect results.
PARALLEL_CRITICAL(monteCarloVolume) { PARALLEL_ATOMIC
total += localTotal; totalHits += hits;
totalHits += hits;
}
} }
const double ratio = const double ratio =
static_cast<double>(totalHits) / static_cast<double>(total); static_cast<double>(totalHits) / static_cast<double>(shotSize);
const double boundingVolume = boundingDx * boundingDy * boundingDz; const double boundingVolume = boundingDx * boundingDy * boundingDz;
return ratio * boundingVolume; return ratio * boundingVolume;
} }
......
#ifndef MANTID_TESTOBJECT__ #ifndef MANTID_TESTOBJECT__
#define MANTID_TESTOBJECT__ #define MANTID_TESTOBJECT__
#include <cxxtest/TestSuite.h>
#include <cmath>
#include <ostream>
#include <vector>
#include <algorithm>
#include <ctime>
#include "boost/shared_ptr.hpp"
#include "boost/make_shared.hpp"
#include "MantidGeometry/Objects/Object.h" #include "MantidGeometry/Objects/Object.h"
#include "MantidGeometry/Surfaces/Cylinder.h" #include "MantidGeometry/Surfaces/Cylinder.h"
#include "MantidGeometry/Surfaces/Sphere.h" #include "MantidGeometry/Surfaces/Sphere.h"
#include "MantidGeometry/Surfaces/Plane.h" #include "MantidGeometry/Surfaces/Plane.h"
...@@ -21,12 +12,25 @@ ...@@ -21,12 +12,25 @@
#include "MantidGeometry/Objects/Track.h" #include "MantidGeometry/Objects/Track.h"
#include "MantidGeometry/Rendering/GluGeometryHandler.h" #include "MantidGeometry/Rendering/GluGeometryHandler.h"
#include "MantidGeometry/Objects/ShapeFactory.h" #include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidKernel/make_unique.h"
#include "MantidKernel/Material.h" #include "MantidKernel/Material.h"
#include "MantidKernel/MersenneTwister.h" #include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/WarningSuppressions.h" #include "MantidKernel/WarningSuppressions.h"
#include "MantidTestHelpers/ComponentCreationHelper.h" #include "MantidTestHelpers/ComponentCreationHelper.h"
#include <cxxtest/TestSuite.h>
#include <cmath>
#include <ostream>
#include <vector>
#include <algorithm>
#include <ctime>
#include "boost/shared_ptr.hpp"
#include "boost/make_shared.hpp"
#include <gmock/gmock.h> #include <gmock/gmock.h>
#include <Poco/DOM/AutoPtr.h>
#include <Poco/DOM/Document.h>
using namespace Mantid; using namespace Mantid;
using namespace Geometry; using namespace Geometry;
...@@ -859,15 +863,83 @@ public: ...@@ -859,15 +863,83 @@ public:
expected, satol); expected, satol);
} }
void testVolumeUnitCube() { void testExactVolumeCuboid() {
Object_sptr cube = createUnitCube(); using namespace Poco::XML;
TS_ASSERT_DELTA(cube->volume(), 1.0, 1e-6) const double width = 1.23;
const double height = 4.98;
const double thickness = 8.14;
AutoPtr<Document> shapeDescription = new Document;
AutoPtr<Element> typeElement = shapeDescription->createElement("type");
typeElement->setAttribute("name", "testCuboid");
AutoPtr<Element> shapeElement = createCuboidTypeElement("cuboid-shape", width, height, thickness, shapeDescription);
typeElement->appendChild(shapeElement);
AutoPtr<Element> algebraElement = shapeDescription->createElement("algebra");
algebraElement->setAttribute("val", "cuboid-shape");
typeElement->appendChild(algebraElement);
ShapeFactory shapeFactory;
auto cuboid = shapeFactory.createShape(typeElement);
const double cuboidVolume = width * height * thickness;
TS_ASSERT_DELTA(cuboid->volume(), cuboidVolume, 1e-6)
}
void testExactVolumeSphere() {
using namespace Poco::XML;
const double radius = 99.9;
AutoPtr<Document> shapeDescription = new Document;
AutoPtr<Element> typeElement = shapeDescription->createElement("type");
typeElement->setAttribute("name", "testSphere");
AutoPtr<Element> shapeElement = createSphereTypeElement("sphere-shape", radius, shapeDescription);
typeElement->appendChild(shapeElement);
AutoPtr<Element> algebraElement = shapeDescription->createElement("algebra");
algebraElement->setAttribute("val", "sphere-shape");
typeElement->appendChild(algebraElement);
ShapeFactory shapeFactory;
auto cuboid = shapeFactory.createShape(typeElement);
const double sphereVolume = 4.0 / 3.0 * M_PI * radius * radius * radius;
TS_ASSERT_DELTA(cuboid->volume(), sphereVolume, 1e-6)
} }
void testVolumeSmallCappedCylinder() { void testExactVolumeCylinder() {
Object_sptr cylinder = createSmallCappedCylinder(); using namespace Poco::XML;
const double correctVolume = M_PI * 0.005 * 0.005 * 0.003; const double radius = 0.99;
TS_ASSERT_DELTA(cylinder->volume(), correctVolume, 1e-6) const double height = 88;
AutoPtr<Document> shapeDescription = new Document;
AutoPtr<Element> typeElement = shapeDescription->createElement("type");
typeElement->setAttribute("name", "testCylinder");
AutoPtr<Element> shapeElement = createCylinderTypeElement("cylinder-shape", height, radius, shapeDescription);
typeElement->appendChild(shapeElement);
AutoPtr<Element> algebraElement = shapeDescription->createElement("algebra");
algebraElement->setAttribute("val", "cylinder-shape");
typeElement->appendChild(algebraElement);
ShapeFactory shapeFactory;
auto cuboid = shapeFactory.createShape(typeElement);
const double cylinderVolume = height * M_PI * radius * radius;
TS_ASSERT_DELTA(cuboid->volume(), cylinderVolume, 1e-6)
}
void testMonteCarloVolume() {
// We use a cuboid with spherical void here.
using namespace Poco::XML;
const double width = 71.99;
const double height = 11.87;
const double thickness = 74.1;
AutoPtr<Document> shapeDescription = new Document;
AutoPtr<Element> typeElement = shapeDescription->createElement("type");
typeElement->setAttribute("name", "testShape");
AutoPtr<Element> shapeElement = createCuboidTypeElement("solid-cuboid", width, height, thickness, shapeDescription);
typeElement->appendChild(shapeElement);
const double radius = 0.47 * std::min(std::min(width, height), thickness);
shapeElement = createSphereTypeElement("void-sphere", radius, shapeDescription);
typeElement->appendChild(shapeElement);
AutoPtr<Element> algebraElement = shapeDescription->createElement("algebra");
algebraElement->setAttribute("val", "solid-cuboid (# void-sphere)");
typeElement->appendChild(algebraElement);
ShapeFactory shapeFactory;
auto cuboid = shapeFactory.createShape(typeElement);
const double cuboidVolume = width * height * thickness;
const double sphereVolume = 4.0 / 3.0 * M_PI * radius * radius * radius;
const double correctVolume = cuboidVolume - sphereVolume;
TS_ASSERT_DELTA(cuboid->volume(), correctVolume, 1e-3 * correctVolume)
} }
void testGetBoundingBoxForCylinder() void testGetBoundingBoxForCylinder()
...@@ -1351,6 +1423,71 @@ private: ...@@ -1351,6 +1423,71 @@ private:
retVal->populate(HexSurMap); retVal->populate(HexSurMap);
return retVal; return retVal;
} }
static Poco::XML::AutoPtr<Poco::XML::Element> createCuboidTypeElement(const std::string &id, const double width, const double height, const double thickness, Poco::XML::AutoPtr<Poco::XML::Document> &document) {
using namespace Poco::XML;
AutoPtr<Element> shapeElement = document->createElement("cuboid");
shapeElement->setAttribute("id", id);
AutoPtr<Element> element = document->createElement("left-front-bottom-point");
element->setAttribute("x", std::to_string(-width / 2));
element->setAttribute("y", std::to_string(-height / 2));
element->setAttribute("z", std::to_string(thickness / 2));
shapeElement->appendChild(element);
element = document->createElement("left-front-top-point");
element->setAttribute("x", std::to_string(-width / 2));
element->setAttribute("y", std::to_string(height / 2));
element->setAttribute("z", std::to_string(thickness / 2));
shapeElement->appendChild(element);
element = document->createElement("left-back-bottom-point");
element->setAttribute("x", std::to_string(-width / 2));
element->setAttribute("y", std::to_string(-height / 2));
element->setAttribute("z", std::to_string(-thickness / 2));
shapeElement->appendChild(element);
element = document->createElement("right-front-bottom-point");
element->setAttribute("x", std::to_string(width / 2));
element->setAttribute("y", std::to_string(-height / 2));
element->setAttribute("z", std::to_string(thickness / 2));
shapeElement->appendChild(element);
return shapeElement;
}
static Poco::XML::AutoPtr<Poco::XML::Element> createSphereTypeElement(const std::string &id, const double radius, Poco::XML::AutoPtr<Poco::XML::Document> &document) {
using namespace Poco::XML;
AutoPtr<Element>shapeElement = document->createElement("sphere");
shapeElement->setAttribute("id", id);
AutoPtr<Element> element = document->createElement("centre");
element->setAttribute("x", "0.0");
element->setAttribute("y", "0.0");
element->setAttribute("z", "0.0");
shapeElement->appendChild(element);
element = document->createElement("radius");
element->setAttribute("val", std::to_string(radius));
shapeElement->appendChild(element);
return shapeElement;
}
static Poco::XML::AutoPtr<Poco::XML::Element> createCylinderTypeElement(const std::string &id, const double height, const double radius, Poco::XML::AutoPtr<Poco::XML::Document> &document) {
using namespace Poco::XML;
AutoPtr<Element>shapeElement = document->createElement("cylinder");
shapeElement->setAttribute("id", id);
AutoPtr<Element> element = document->createElement("centre-of-bottom-base");
element->setAttribute("x", std::to_string(-height / 2));
element->setAttribute("y", "0.0");
element->setAttribute("z", "0.0");
shapeElement->appendChild(element);
element = document->createElement("axis");
element->setAttribute("x", "1.0");
element->setAttribute("y", "0.0");
element->setAttribute("z", "0.0");
shapeElement->appendChild(element);
element = document->createElement("radius");
element->setAttribute("val", std::to_string(radius));
shapeElement->appendChild(element);
element = document->createElement("height");
element->setAttribute("val", std::to_string(height));
shapeElement->appendChild(element);
return shapeElement;
}
}; };
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment