Unverified Commit a3d090d9 authored by Gigg, Martyn Anthony's avatar Gigg, Martyn Anthony Committed by GitHub
Browse files

Merge pull request #28648 from mantidproject/28339_ellipsoid_integration_in_IntegratePeaksMD

parents 3c7dcd7d f5bd37e7
......@@ -38,29 +38,36 @@ DECLARE_VECTOR_PARAMETER(DimensionsUsedVectorParam, bool)
*/
class DLLExport CoordTransformDistance : public Mantid::API::CoordTransform {
public:
CoordTransformDistance(const size_t inD, const coord_t *center,
const bool *dimensionsUsed, const size_t outD = 1);
CoordTransformDistance(
const size_t inD, const coord_t *center, const bool *dimensionsUsed,
const size_t outD = 1,
const std::vector<Kernel::V3D> &eigenvects = std::vector<Kernel::V3D>(0),
const std::vector<double> &eigenvals = std::vector<double>(0, 0.0));
CoordTransform *clone() const override;
~CoordTransformDistance() override;
std::string toXMLString() const override;
std::string id() const override;
void apply(const coord_t *inputVector, coord_t *outVector) const override;
/// Return the center coordinate array
const coord_t *getCenter() { return m_center; }
const std::vector<coord_t> &getCenter() { return m_center; }
/// Return the dimensions used bool array
const bool *getDimensionsUsed() { return m_dimensionsUsed; }
const std::vector<bool> &getDimensionsUsed() { return m_dimensionsUsed; }
protected:
/// Coordinates at the center
coord_t *m_center;
std::vector<coord_t> m_center;
/// Parmeter where True is set for those dimensions that are considered when
/// calculating distance
bool *m_dimensionsUsed;
};
std::vector<bool> m_dimensionsUsed;
// Eigenvectors and radii for nd ellipsoid
std::vector<Kernel::V3D> m_eigenvects;
std::vector<double> m_eigenvals;
double m_maxEigenval;
};
} // namespace DataObjects
} // namespace Mantid
......@@ -22,7 +22,8 @@ namespace DataObjects {
class MANTID_DATAOBJECTS_DLL FakeMD {
public:
FakeMD(const std::vector<double> &uniformParams,
const std::vector<double> &peakParams, const int randomSeed,
const std::vector<double> &peakParams,
const std::vector<double> &ellipsoidParams, const int randomSeed,
const bool randomizeSignal);
void fill(const API::IMDEventWorkspace_sptr &workspace);
......@@ -33,6 +34,8 @@ private:
template <typename MDE, size_t nd>
void addFakePeak(typename MDEventWorkspace<MDE, nd>::sptr ws);
template <typename MDE, size_t nd>
void addFakeEllipsoid(typename MDEventWorkspace<MDE, nd>::sptr ws);
template <typename MDE, size_t nd>
void addFakeUniformData(typename MDEventWorkspace<MDE, nd>::sptr ws);
template <typename MDE, size_t nd>
......@@ -47,6 +50,7 @@ private:
//------------------ Member variables ------------------------------------
std::vector<double> m_uniformParams;
std::vector<double> m_peakParams;
std::vector<double> m_ellipsoidParams;
const int m_randomSeed;
const bool m_randomizeSignal;
mutable std::vector<detid_t> m_detIDs;
......
......@@ -28,19 +28,29 @@ namespace DataObjects {
*dimensions that are considered when
* calculating distance.
* @param outD :: # of output dimensions
* @param eigenvects :: eigenvectors of an ellipsoid (if used)
* @param eigenvals :: variances along each eigenvector (if used)
* @return
*/
CoordTransformDistance::CoordTransformDistance(const size_t inD,
const coord_t *center,
const bool *dimensionsUsed,
const size_t outD)
CoordTransformDistance::CoordTransformDistance(
const size_t inD, const coord_t *center, const bool *dimensionsUsed,
const size_t outD, const std::vector<Kernel::V3D> &eigenvects,
const std::vector<double> &eigenvals)
: CoordTransform(inD, outD) {
// Create and copy the arrays.
m_center = new coord_t[inD];
m_dimensionsUsed = new bool[inD];
m_center.reserve(inD);
m_dimensionsUsed.reserve(inD);
m_eigenvals.reserve(inD);
m_maxEigenval = 0.0;
for (size_t d = 0; d < inD; d++) {
m_center[d] = center[d];
m_dimensionsUsed[d] = dimensionsUsed[d];
m_center.push_back(center[d]);
m_dimensionsUsed.push_back(dimensionsUsed[d]);
if (eigenvals.size() == inD && eigenvects.size() == inD) {
// coord transform for n-ellipsoid specified
m_eigenvals.push_back(eigenvals[d]);
m_eigenvects.push_back(eigenvects[d]);
m_maxEigenval = std::max(m_maxEigenval, eigenvals[d]);
}
}
}
......@@ -48,16 +58,7 @@ CoordTransformDistance::CoordTransformDistance(const size_t inD,
/** Virtual cloner
* @return a copy of this object */
CoordTransform *CoordTransformDistance::clone() const {
auto out = new CoordTransformDistance(inD, m_center, m_dimensionsUsed);
return out;
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
CoordTransformDistance::~CoordTransformDistance() {
delete[] m_center;
delete[] m_dimensionsUsed;
return new CoordTransformDistance(*this);
}
//----------------------------------------------------------------------------------------------
......@@ -74,10 +75,26 @@ void CoordTransformDistance::apply(const coord_t *inputVector,
coord_t *outVector) const {
if (outD == 1) {
coord_t distanceSquared = 0;
for (size_t d = 0; d < inD; d++) {
if (m_dimensionsUsed[d]) {
coord_t dist = inputVector[d] - m_center[d];
distanceSquared += (dist * dist);
if (m_eigenvals.size() == 3) {
// ellipsoid transfrom with ratio of axes given by standard deviation
// i.e. sqrt(eigenvals)
for (size_t d = 0; d < inD; d++) {
// do dot prod with eigenvector
coord_t dist = 0.0;
for (size_t dd = 0; dd < inD; dd++) {
dist += static_cast<coord_t>(m_eigenvects[d][dd]) *
(inputVector[dd] - m_center[dd]);
}
distanceSquared += (dist * dist) *
static_cast<coord_t>(m_maxEigenval / m_eigenvals[d]);
}
} else {
// nd spherical transform
for (size_t d = 0; d < inD; d++) {
if (m_dimensionsUsed[d]) {
coord_t dist = inputVector[d] - m_center[d];
distanceSquared += (dist * dist);
}
}
}
/// Return the only output dimension
......
......@@ -17,6 +17,8 @@
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/Utils.h"
#include "boost/math/distributions.hpp"
namespace Mantid {
namespace DataObjects {
......@@ -28,19 +30,22 @@ using Kernel::ThreadSchedulerFIFO;
* @param uniformParams Add a uniform, randomized distribution of events
* @param peakParams Add a peak with a normal distribution around a central
point
* @param ellipsoidParams Add a multivariate gaussian peak (ellipsoid)
* @param randomSeed Seed int for the random number generator
* @param randomizeSignal If true, the events' signal and error values will be "
randomized around 1.0+-0.5
*/
FakeMD::FakeMD(const std::vector<double> &uniformParams,
const std::vector<double> &peakParams, const int randomSeed,
const std::vector<double> &peakParams,
const std::vector<double> &ellipsoidParams, const int randomSeed,
const bool randomizeSignal)
: m_uniformParams(uniformParams), m_peakParams(peakParams),
m_randomSeed(randomSeed), m_randomizeSignal(randomizeSignal), m_detIDs(),
m_randGen(1), m_uniformDist() {
if (uniformParams.empty() && peakParams.empty()) {
throw std::invalid_argument(
"You must specify at least one of peakParams or uniformParams");
m_ellipsoidParams(ellipsoidParams), m_randomSeed(randomSeed),
m_randomizeSignal(randomizeSignal), m_detIDs(), m_randGen(1),
m_uniformDist() {
if (uniformParams.empty() && peakParams.empty() && ellipsoidParams.empty()) {
throw std::invalid_argument("You must specify at least one of peakParams, "
"ellipsoidParams or uniformParams");
}
}
......@@ -53,6 +58,7 @@ void FakeMD::fill(const API::IMDEventWorkspace_sptr &workspace) {
setupDetectorCache(*workspace);
CALL_MDEVENT_FUNCTION(this->addFakePeak, workspace)
CALL_MDEVENT_FUNCTION(this->addFakeEllipsoid, workspace)
CALL_MDEVENT_FUNCTION(this->addFakeUniformData, workspace)
// Mark that events were added, so the file back end (if any) needs updating
......@@ -150,6 +156,120 @@ void FakeMD::addFakePeak(typename MDEventWorkspace<MDE, nd>::sptr ws) {
ws->refreshCache();
}
/**
* Function that adds a fake single-crystal peak with a multivariate normal
* distribution (an ellipsoid) around a central point (x1,...x_N).
* The ellipsoid is defined by N eigenvectors with N elements and
* N eigenvalues which correpsond to the variance along the principal axes.
* The peak is generated from an affine transformation of a uniform normal
* distirbution of variance = 1.
*
* @param ws A pointer to the workspace that receives the events
*/
template <typename MDE, size_t nd>
void FakeMD::addFakeEllipsoid(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_ellipsoidParams.empty())
return;
if (m_ellipsoidParams.size() != 2 + 2 * nd + nd * nd)
throw std::invalid_argument(
"EllipsoidParams: incorrect number of parameters.");
if (m_ellipsoidParams[0] <= 0)
throw std::invalid_argument(
"EllipsoidParams: number_of_events needs to be > 0");
// extract input parameters
auto numEvents = size_t(m_ellipsoidParams[0]);
coord_t center[nd];
Kernel::Matrix<double> Evec(nd, nd); // hold eigenvectors
Kernel::Matrix<double> stds(nd,
nd); // hold sqrt(eigenvals) standard devs on diag
for (size_t n = 0; n < nd; n++) {
center[n] = static_cast<coord_t>(m_ellipsoidParams[n + 1]);
// get row/col index for eigenvector matrix
for (size_t d = 0; d < nd; d++) {
Evec[d][n] = m_ellipsoidParams[1 + nd + n * nd + d];
}
stds[n][n] =
sqrt(m_ellipsoidParams[m_ellipsoidParams.size() - (1 + nd) + n]);
}
auto doCounts = m_ellipsoidParams[m_ellipsoidParams.size() - 1];
// get affine transformation that maps unit variance spherical
// normal dist to ellipsoid
auto A = Evec * stds;
// calculate inverse of covariance matrix (if necassary)
Kernel::Matrix<double> invCov(nd, nd);
if (doCounts > 0) {
auto var = stds * stds;
// copy Evec to a matrix to hold inverse
Kernel::Matrix<double> invEvec(Evec.getVector()); // hold eigenvectors
// invert Evec matrix
invEvec.Invert();
// find covariance matrix to invert
invCov = Evec * var * invEvec; // covar mat
// invert covar matrix
invCov.Invert();
}
// get chi-squared boost function
boost::math::chi_squared chisq(nd);
// prepare random number generator
std::mt19937 rng(static_cast<unsigned int>(m_randomSeed));
std::normal_distribution<double> d(0.0, 1.0); // mean = 0, std = 1
// Inserter to help choose the correct event type
auto eventHelper =
MDEventInserter<typename MDEventWorkspace<MDE, nd>::sptr>(ws);
for (size_t iEvent = 0; iEvent < numEvents; ++iEvent) {
// sample uniform normal distribution
std::vector<double> pos;
for (size_t n = 0; n < nd; n++) {
pos.push_back(d(rng));
}
// apply affine transformation
pos = A * pos;
// calculate counts
float signal = 1.0;
float errorSquared = 1.0;
if (doCounts > 0) {
// calculate Mahalanobis distance
// https://en.wikipedia.org/wiki/Mahalanobis_distance
// md = sqrt(pos.T * invCov * pos)
auto tmp = invCov * pos;
double mdsq = 0.0;
for (size_t n = 0; n < nd; n++) {
mdsq += pos[n] * tmp[n];
}
// for a multivariate normal dist m-dist is distribute
// as chi-squared pdf with nd degrees of freedom
signal = static_cast<float>(boost::math::pdf(chisq, sqrt(mdsq)));
errorSquared = signal;
}
// convert pos to coord_t and offset by center
coord_t eventCenter[nd];
for (size_t n = 0; n < nd; n++) {
eventCenter[n] = static_cast<coord_t>(pos[n] + center[n]);
}
// add event (need to convert pos to coord_t)
eventHelper.insertMDEvent(signal, errorSquared, 0, pickDetectorID(),
eventCenter); // 0 = run index
}
ws->splitBox();
auto *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
ws->splitAllIfNeeded(ts);
tp.joinAll();
ws->refreshCache();
}
/**
* Function makes up a fake uniform event data and adds it to the workspace.
* @param ws
......
......@@ -21,7 +21,8 @@ using Mantid::API::CoordTransform;
class CoordTransformDistanceTest : public CxxTest::TestSuite {
public:
/** Helper to compare two "vectors" (bare float arrays) */
void compare(size_t numdims, coord_t *value, const coord_t *expected) {
void compare(size_t numdims, coord_t *value,
const std::vector<coord_t> expected) {
for (size_t i = 0; i < numdims; i++)
TS_ASSERT_DELTA(value[i], expected[i], 1e-5);
}
......@@ -31,14 +32,14 @@ public:
bool used[4] = {true, false, true, true};
CoordTransformDistance ct(4, center, used);
// A copy was made
const coord_t *transformCentres = ct.getCenter();
TS_ASSERT_DIFFERS(transformCentres, center);
const bool *usedDims = ct.getDimensionsUsed();
TS_ASSERT_DIFFERS(ct.getDimensionsUsed(), used);
const std::vector<coord_t> transformCentres = ct.getCenter();
const std::vector<bool> usedDims = ct.getDimensionsUsed();
// Contents are good
compare(4, center, ct.getCenter());
for (size_t i = 0; i < 4; i++)
TS_ASSERT_EQUALS(used[i], usedDims[i]);
for (size_t i = 0; i < 4; i++) {
const bool dimUsed = usedDims[i];
TS_ASSERT_EQUALS(used[i], dimUsed);
}
}
/** Clone then apply */
......@@ -91,6 +92,52 @@ public:
TS_ASSERT_DELTA(out, 4.0, 1e-5);
}
/** Calculate the distance (squared) for ellipsoid*/
void test_distance_ellipsoid() {
// Build it
coord_t center[3] = {1, 2, 3};
bool used[3] = {true, true, true};
std::vector<Kernel::V3D> eigenvects;
std::vector<double> eigenvals;
eigenvects.push_back(Kernel::V3D(1.0, 0.0, 0.0));
eigenvals.push_back(4);
eigenvects.push_back(Kernel::V3D(0.0, 1.0, 0.0));
eigenvals.push_back(1);
eigenvects.push_back(Kernel::V3D(0.0, 0.0, 1.0));
eigenvals.push_back(1);
CoordTransformDistance ct(3, center, used, 1, /*outD*/
eigenvects, eigenvals);
coord_t out = 0;
coord_t in1[3] = {1, 2, 3};
TS_ASSERT_THROWS_NOTHING(ct.apply(in1, &out));
TS_ASSERT_DELTA(out, 0.0, 1e-5);
coord_t in2[3] = {1, 2, 4};
TS_ASSERT_THROWS_NOTHING(ct.apply(in2, &out));
TS_ASSERT_DELTA(out, 4.0, 1e-5);
coord_t in3[3] = {1, 3, 3};
TS_ASSERT_THROWS_NOTHING(ct.apply(in3, &out));
TS_ASSERT_DELTA(out, 4.0, 1e-5);
coord_t in4[3] = {3, 2, 3};
TS_ASSERT_THROWS_NOTHING(ct.apply(in4, &out));
TS_ASSERT_DELTA(out, 4.0, 1e-5);
coord_t in5[3] = {1, static_cast<coord_t>(2 + M_SQRT1_2),
static_cast<coord_t>(3 - M_SQRT1_2)};
TS_ASSERT_THROWS_NOTHING(ct.apply(in5, &out));
TS_ASSERT_DELTA(out, 4.0, 1e-5);
coord_t in6[3] = {1, static_cast<coord_t>(2 + M_SQRT2),
static_cast<coord_t>(3 - M_SQRT2)};
TS_ASSERT_THROWS_NOTHING(ct.apply(in6, &out));
TS_ASSERT_DELTA(out, 16.0, 1e-5);
}
/** Test serialization */
void test_to_xml_string() {
std::string expectedResult =
......
......@@ -31,12 +31,13 @@ public:
const std::vector<double> peakParams;
const std::vector<double> uniformParams;
const std::vector<double> ellipsoidParams;
const int randomSeed(0);
const bool randomizeSignal(false);
TS_ASSERT_THROWS(
FakeMD(uniformParams, peakParams, randomSeed, randomizeSignal),
const std::invalid_argument &);
TS_ASSERT_THROWS(FakeMD(uniformParams, peakParams, ellipsoidParams,
randomSeed, randomizeSignal),
const std::invalid_argument &);
}
//---------------------------------------------------------------------------
......@@ -54,13 +55,44 @@ public:
const std::vector<double> peakParams = {1000.0, 5.0, 5.0, 5.0, 1.0};
const std::vector<double> uniformParams = {10000.0};
const std::vector<double> ellipsoidParams = {500.0, 5.0, 5.0, 5.0, 1.0, 0.0,
0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
1.0, 0.5, 0.5, 0.5, -1.0};
const int randomSeed(0);
const bool randomizeSignal(false);
FakeMD faker(uniformParams, peakParams, randomSeed, randomizeSignal);
FakeMD faker(uniformParams, peakParams, ellipsoidParams, randomSeed,
randomizeSignal);
faker.fill(fakeData);
// Now there are 11000 more points.
TS_ASSERT_EQUALS(fakeData->getNPoints(), 12000);
TS_ASSERT_EQUALS(fakeData->getNPoints(), 12500);
}
void test_ellipsoid_counts() {
using Mantid::DataObjects::FakeMD;
using Mantid::DataObjects::MDEventsTestHelper::makeMDEW;
// Destination workspace
auto fakeData = makeMDEW<3>(10, 0.0, 10.0, 0);
const std::vector<double> peakParams;
const std::vector<double> uniformParams;
const std::vector<double> ellipsoidParams = {
2000.0, 5.0, 5.0, 5.0, 1.0, 0.0, 0.0, 0.0, 1.0,
0.0, 0.0, 0.0, 1.0, 0.5, 0.5, 0.5, 1.0};
const int randomSeed(0);
const bool randomizeSignal(false);
FakeMD faker(uniformParams, peakParams, ellipsoidParams, randomSeed,
randomizeSignal);
faker.fill(fakeData);
auto Npts = fakeData->getNPoints();
TS_ASSERT_EQUALS(Npts, ellipsoidParams[0]);
// avg of counts converges to 0.2175 for 3D multivariate gaussian
TS_ASSERT_DELTA(fakeData->getBox()->getSignal(),
static_cast<double>(Npts) * 0.2175,
static_cast<double>(Npts) * 0.0015);
}
void test_exec_randomizeSignal() {
......@@ -73,10 +105,12 @@ public:
const std::vector<double> peakParams = {100.0, 5.0, 5.0, 5.0, 1.0};
const std::vector<double> uniformParams = {100.0};
const std::vector<double> ellipsoidParams;
const int randomSeed(0);
const bool randomizeSignal(true);
FakeMD faker(uniformParams, peakParams, randomSeed, randomizeSignal);
FakeMD faker(uniformParams, peakParams, ellipsoidParams, randomSeed,
randomizeSignal);
faker.fill(fakeData);
// Now there are 200 more points.
......@@ -104,10 +138,12 @@ public:
const std::vector<double> peakParams;
const std::vector<double> uniformParams = {-1000.0};
const std::vector<double> ellipsoidParams;
const int randomSeed(0);
const bool randomizeSignal(false);
FakeMD faker(uniformParams, peakParams, randomSeed, randomizeSignal);
FakeMD faker(uniformParams, peakParams, ellipsoidParams, randomSeed,
randomizeSignal);
faker.fill(fakeData);
// Now there are 1000 more points.
......@@ -137,10 +173,12 @@ public:
const std::vector<double> peakParams;
const std::vector<double> uniformParams = {-1000.0};
const std::vector<double> ellipsoidParams;
const int randomSeed(0);
const bool randomizeSignal(false);
FakeMD faker(uniformParams, peakParams, randomSeed, randomizeSignal);
FakeMD faker(uniformParams, peakParams, ellipsoidParams, randomSeed,
randomizeSignal);
faker.fill(fakeData);
TS_ASSERT_EQUALS(1000, fakeData->getNEvents());
......
......@@ -189,6 +189,7 @@ set(
inc/MantidMDAlgorithms/LoadSQW.h
inc/MantidMDAlgorithms/LoadSQW2.h
inc/MantidMDAlgorithms/LogarithmMD.h
inc/MantidMDAlgorithms/MDBoxMaskFunction.h
inc/MantidMDAlgorithms/MDEventTreeBuilder.h
inc/MantidMDAlgorithms/MDEventWSWrapper.h
inc/MantidMDAlgorithms/MDNorm.h
......@@ -305,6 +306,7 @@ set(TEST_FILES
LoadSQW2Test.h
LoadSQWTest.h
LogarithmMDTest.h
MDBoxMaskFunctionTest.h
MDEventWSWrapperTest.h
MDNormDirectSCTest.h
MDNormSCDTest.h
......
......@@ -58,6 +58,21 @@ private:
/// Input MDEventWorkspace
Mantid::API::IMDEventWorkspace_sptr inWS;
// find the eigenvectors and eigenvalues that diagonalise the covariance
// matrix that defines an ellipsoid.
template <typename MDE, size_t nd>
void findEllipsoid(typename DataObjects::MDEventWorkspace<MDE, nd>::sptr ws,
const Mantid::API::CoordTransform &getRadiusSq,
const Mantid::Kernel::V3D &pos,
const coord_t &radiusSquared, const bool &qAxisBool,
const double &bgDensity,
std::vector<Mantid::Kernel::V3D> &eigenvects,
std::vector<double> &eigenvals);
// get matrix to transform from Qlab to plane perp to Q
void getPinv(const Mantid::Kernel::V3D &q,
Mantid::Kernel::Matrix<double> &Pinv);
/// Calculate if this Q is on a detector
void calculateE1(const Geometry::DetectorInfo &detectorInfo);
double detectorQ(Mantid::Kernel::V3D QLabFrame, double r);
......
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 2020 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//--------------------------------------------------------------------------------------------------
// Includes
//--------------------------------------------------------------------------------------------------
#pragma once
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include <MantidKernel/V3D.h>
namespace Mantid::Geometry::MDAlgorithms {
class MDBoxMaskFunction : public Mantid::Geometry::MDImplicitFunction {
private:
Mantid::Kernel::V3D m_pos;
double m_radiusSquared;
public:
// constructor
MDBoxMaskFunction(const Mantid::Kernel::V3D &pos,
const double &radiusSquared) {
m_pos = pos;
m_radiusSquared = radiusSquared;
}
using MDImplicitFunction::isPointContained; // Avoids Intel compiler
// warning.
bool isPointContained(const coord_t *coords) override {
double sum = 0;
for (size_t i = 0; i < 3; i++) {