forked from npshub/mantid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
FakeMD.cpp
447 lines (387 loc) · 15.9 KB
/
FakeMD.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 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
//--------------------------------------------------------------------------------------------------
#include "MantidDataObjects/FakeMD.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDEventInserter.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/ThreadPool.h"
#include "MantidKernel/ThreadScheduler.h"
#include "MantidKernel/Utils.h"
#include "MantidKernel/normal_distribution.h"
#include "boost/math/distributions.hpp"
namespace Mantid::DataObjects {
using Kernel::ThreadPool;
using Kernel::ThreadSchedulerFIFO;
/**
* Constructor
* @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 std::vector<double> &ellipsoidParams, const int randomSeed, const bool randomizeSignal)
: m_uniformParams(uniformParams), m_peakParams(peakParams), 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");
}
}
/**
* Add the fake data to the given workspace
* @param workspace A pointer to MD event workspace to fill using the object
* parameters
*/
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
workspace->setFileNeedsUpdating(true);
}
/**
* Setup a detector cache for randomly picking IDs from the first
* instrument in the ExperimentInfo list.
* @param workspace The input workspace
*/
void FakeMD::setupDetectorCache(const API::IMDEventWorkspace &workspace) {
try {
auto inst = workspace.getExperimentInfo(0)->getInstrument();
m_detIDs = inst->getDetectorIDs(true); // true=skip monitors
size_t max = m_detIDs.size() - 1;
m_uniformDist = Kernel::uniform_int_distribution<size_t>(0, max); // Includes max
} catch (std::invalid_argument &) {
}
}
/**
* Function generates random uniformly distributed events within an nD-sphere
* to simulate a single-crystal peak and adds it to the workspace.
*
* @param ws A pointer to the workspace that receives the events
*/
template <typename MDE, size_t nd> void FakeMD::addFakePeak(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_peakParams.empty())
return;
if (m_peakParams.size() != nd + 2)
throw std::invalid_argument("PeakParams needs to have ndims+2 arguments.");
if (m_peakParams[0] <= 0)
throw std::invalid_argument("PeakParams: number_of_events needs to be > 0");
auto num = size_t(m_peakParams[0]);
// Width of the peak
auto desiredRadius = m_peakParams.back();
std::mt19937 rng(static_cast<unsigned int>(m_randomSeed));
std::uniform_real_distribution<coord_t> flat(0, 1.0);
Kernel::normal_distribution<coord_t> normal(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 i = 0; i < num; ++i) {
// generate point along radius with ^1/n scaling for uniformity (i.e. Prob(point < r) ~ r^nd)
auto radius_frac = flat(rng);
auto radius = static_cast<coord_t>(desiredRadius * pow(radius_frac, 1.0f / nd));
// to get direction generate uniformly distributed points on surface of a unit N sphere using
// uniformly dist. rand numbers as per http://corysimon.github.io/articles/uniformdistn-on-sphere/
coord_t centers[nd];
coord_t modulusSq = 0;
while (modulusSq < 1E-6) {
// small modulus may cause issues with floating point precision when dividing later
for (size_t d = 0; d < nd; d++) {
centers[d] = normal(rng);
modulusSq += centers[d] * centers[d];
}
}
// now normalise by modulus, scale by radius and translate to peak centre
auto inverseModulus = static_cast<coord_t>(1.0f / sqrt(modulusSq));
for (size_t d = 0; d < nd; d++) {
// Multiply by the scaling and the desired peak radius
centers[d] *= radius * inverseModulus;
// Also offset by the center of the peak, as taken in Params
centers[d] += static_cast<coord_t>(m_peakParams[d + 1]);
}
// Default or randomized error/signal
float signal = 1.0;
float errorSquared = 1.0;
if (m_randomizeSignal) {
signal = float(0.5 + flat(rng));
errorSquared = float(0.5 + flat(rng));
}
// Create and add the event.
eventHelper.insertMDEvent(signal, errorSquared, 0, 0, pickDetectorID(),
centers); // 0 = associated experiment-info index
}
ws->splitBox();
auto *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
ws->splitAllIfNeeded(ts);
tp.joinAll();
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));
Kernel::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, 0, pickDetectorID(),
eventCenter); // 0 = associated experiment-info 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
*/
template <typename MDE, size_t nd> void FakeMD::addFakeUniformData(typename MDEventWorkspace<MDE, nd>::sptr ws) {
if (m_uniformParams.empty())
return;
bool randomEvents = true;
if (m_uniformParams[0] < 0) {
randomEvents = false;
m_uniformParams[0] = -m_uniformParams[0];
}
if (m_uniformParams.size() == 1) {
if (randomEvents) {
for (size_t d = 0; d < nd; ++d) {
m_uniformParams.emplace_back(ws->getDimension(d)->getMinimum());
m_uniformParams.emplace_back(ws->getDimension(d)->getMaximum());
}
} else // regular events
{
auto nPoints = size_t(m_uniformParams[0]);
double Vol = 1;
for (size_t d = 0; d < nd; ++d)
Vol *= (ws->getDimension(d)->getMaximum() - ws->getDimension(d)->getMinimum());
if (Vol == 0 || Vol > std::numeric_limits<float>::max())
throw std::invalid_argument(" Domain ranges are not defined properly for workspace: " + ws->getName());
double dV = Vol / double(nPoints);
double delta0 = std::pow(dV, 1. / double(nd));
for (size_t d = 0; d < nd; ++d) {
double min = ws->getDimension(d)->getMinimum();
m_uniformParams.emplace_back(min * (1 + FLT_EPSILON) - min + FLT_EPSILON);
double extent = ws->getDimension(d)->getMaximum() - min;
auto nStrides = size_t(extent / delta0);
if (nStrides < 1)
nStrides = 1;
m_uniformParams.emplace_back(extent / static_cast<double>(nStrides));
}
}
}
if ((m_uniformParams.size() != 1 + nd * 2))
throw std::invalid_argument("UniformParams: needs to have ndims*2+1 arguments ");
if (randomEvents)
addFakeRandomData<MDE, nd>(m_uniformParams, ws);
else
addFakeRegularData<MDE, nd>(m_uniformParams, ws);
ws->splitBox();
auto *ts = new ThreadSchedulerFIFO();
ThreadPool tp(ts);
ws->splitAllIfNeeded(ts);
tp.joinAll();
ws->refreshCache();
}
/**
* Add fake randomized data to the workspace
* @param params A reference to the parameter vector
* @param ws The workspace to hold the data
*/
template <typename MDE, size_t nd>
void FakeMD::addFakeRandomData(const std::vector<double> ¶ms, typename MDEventWorkspace<MDE, nd>::sptr ws) {
auto num = size_t(params[0]);
if (num == 0)
throw std::invalid_argument(" number of distributed events can not be equal to 0");
// Inserter to help choose the correct event type
auto eventHelper = MDEventInserter<typename MDEventWorkspace<MDE, nd>::sptr>(ws);
// Array of distributions for each dimension
std::mt19937 rng(static_cast<unsigned int>(m_randomSeed));
std::array<std::uniform_real_distribution<double>, nd> gens;
for (size_t d = 0; d < nd; ++d) {
double min = params[d * 2 + 1];
double max = params[d * 2 + 2];
if (max <= min)
throw std::invalid_argument("UniformParams: min must be < max for all dimensions.");
gens[d] = std::uniform_real_distribution<double>(min, max);
}
// Unit-size randomizer
std::uniform_real_distribution<double> flat(0, 1.0); // Random from 0 to 1.0
// Create all the requested events
for (size_t i = 0; i < num; ++i) {
coord_t centers[nd];
for (size_t d = 0; d < nd; d++) {
centers[d] = static_cast<coord_t>(gens[d](rng));
}
// Default or randomized error/signal
float signal = 1.0;
float errorSquared = 1.0;
if (m_randomizeSignal) {
signal = float(0.5 + flat(rng));
errorSquared = float(0.5 + flat(rng));
}
// Create and add the event.
eventHelper.insertMDEvent(signal, errorSquared, 0, 0, pickDetectorID(),
centers); // 0 = associated experiment-info index
}
}
template <typename MDE, size_t nd>
void FakeMD::addFakeRegularData(const std::vector<double> ¶ms, typename MDEventWorkspace<MDE, nd>::sptr ws) {
// the parameters for regular distribution of events over the box
std::vector<double> startPoint(nd), delta(nd);
std::vector<size_t> indexMax(nd);
size_t gridSize(0);
auto num = size_t(params[0]);
if (num == 0)
throw std::invalid_argument(" number of distributed events can not be equal to 0");
// Inserter to help choose the correct event type
auto eventHelper = MDEventInserter<typename MDEventWorkspace<MDE, nd>::sptr>(ws);
gridSize = 1;
for (size_t d = 0; d < nd; ++d) {
double min = ws->getDimension(d)->getMinimum();
double max = ws->getDimension(d)->getMaximum();
double shift = params[d * 2 + 1];
double step = params[d * 2 + 2];
if (shift < 0)
shift = 0;
if (shift >= step)
shift = step * (1 - FLT_EPSILON);
startPoint[d] = min + shift;
if ((startPoint[d] < min) || (startPoint[d] >= max))
throw std::invalid_argument("RegularData: starting point must be within "
"the box for all dimensions.");
if (step <= 0)
throw(std::invalid_argument("Step of the regular grid is less or equal to 0"));
indexMax[d] = size_t((max - min) / step);
if (indexMax[d] == 0)
indexMax[d] = 1;
// deal with round-off errors
while ((startPoint[d] + double(indexMax[d] - 1) * step) >= max)
step *= (1 - FLT_EPSILON);
delta[d] = step;
gridSize *= indexMax[d];
}
// Create all the requested events
size_t cellCount(0);
for (size_t i = 0; i < num; ++i) {
coord_t centers[nd];
auto indexes = Kernel::Utils::getIndicesFromLinearIndex(cellCount, indexMax);
++cellCount;
if (cellCount >= gridSize)
cellCount = 0;
for (size_t d = 0; d < nd; d++) {
centers[d] = coord_t(startPoint[d] + delta[d] * double(indexes[d]));
}
// Default or randomized error/signal
float signal = 1.0;
float errorSquared = 1.0;
// Create and add the event.
eventHelper.insertMDEvent(signal, errorSquared, 0, 0, pickDetectorID(),
centers); // 0 = associated experiment-info index
}
}
/**
* Pick a detector ID for a particular event
* @returns A detector ID randomly selected from the instrument
*/
detid_t FakeMD::pickDetectorID() {
if (m_detIDs.empty()) {
return -1;
} else {
return m_detIDs[m_uniformDist(m_randGen)];
}
}
} // namespace Mantid::DataObjects