19#include <boost/optional/optional.hpp>
27using API::IPeaksWorkspace;
29using Geometry::IndexingUtils;
31using Geometry::OrientedLattice;
37const auto OPTIMIZE_UB_ATTEMPTS{4};
40const std::string PEAKSWORKSPACE{
"PeaksWorkspace"};
42const std::string SATE_TOLERANCE{
"ToleranceForSatellite"};
43const std::string ROUNDHKLS{
"RoundHKLs"};
44const std::string COMMONUB{
"CommonUBForAll"};
45const std::string SAVEMODINFO{
"SaveModulationInfo"};
46const std::string AVERAGE_ERR{
"AverageError"};
47const std::string NUM_INDEXED{
"NumIndexed"};
48const std::string MAIN_NUM_INDEXED{
"MainNumIndexed"};
49const std::string SATE_NUM_INDEXED{
"SateNumIndexed"};
50const std::string MAIN_ERR{
"MainError"};
51const std::string SATE_ERR{
"SatelliteError"};
53struct SatelliteIndexingArgs {
60struct IndexPeaksArgs {
66 static IndexPeaksArgs parse(
const API::Algorithm &alg) {
72 std::vector<V3D> modVectorsToUse;
73 modVectorsToUse.reserve(3);
74 bool crossTermToUse{
false};
84 if (maxOrderFromAlg > 0 && modVectorsToUse.size() == 0) {
87 maxOrderToUse = maxOrderFromAlg;
88 }
else if (maxOrderFromAlg == 0 && modVectorsToUse.size() == 0) {
90 const auto &lattice = peaksWS->sample().getOrientedLattice();
91 crossTermToUse = lattice.getCrossTerm();
92 maxOrderToUse = lattice.getMaxOrder();
96 if (maxOrderToUse > 0) {
97 modVectorsToUse =
validModulationVectors(lattice.getModVec(0), lattice.getModVec(1), lattice.getModVec(2));
102 maxOrderToUse = maxOrderFromAlg;
108 alg.getProperty(ROUNDHKLS),
109 alg.getProperty(COMMONUB),
110 alg.getProperty(SAVEMODINFO),
111 SatelliteIndexingArgs{alg.getProperty(SATE_TOLERANCE), maxOrderToUse, modVectorsToUse, crossTermToUse}};
126struct PeakIndexingStats {
127 PeakIndexingStats &
operator+=(
const PeakIndexingStats &
rhs) {
140struct CombinedIndexingStats {
141 CombinedIndexingStats &
operator+=(
const CombinedIndexingStats &
rhs) {
148 inline int totalNumIndexed()
const {
return main.numIndexed +
satellites.numIndexed; }
150 inline double mainError()
const {
151 if (
main.numIndexed == 0)
153 return main.error /
main.numIndexed;
156 inline double satelliteError()
const {
162 inline double averageError()
const {
return (
main.error +
satellites.error) / totalNumIndexed(); }
177 double errorAtStart{0.0};
178 std::vector<V3D> millerIndices;
179 millerIndices.reserve(qSample.size());
180 const int numIndexedAtStart =
183 if (numIndexedAtStart < 3) {
188 for (
auto i = 0; i < OPTIMIZE_UB_ATTEMPTS; ++i) {
196 optimizedUB = ubOrig;
199 double errorInLoop{0.0};
200 const int numIndexedInLoop =
202 if (numIndexedInLoop < numIndexedAtStart)
209using IndexedSatelliteInfo = std::tuple<V3D, V3D, double>;
234boost::optional<IndexedSatelliteInfo> indexSatellite(
const V3D &mainHKL,
const int maxOrder,
238 bool foundSatellite{
false};
239 V3D indexedIntHKL, indexedMNP;
240 for (
const auto &mnpOffset : offsets) {
241 const auto candidateIntHKL = mainHKL - std::get<3>(mnpOffset);
242 const V3D candidateMNP{std::get<0>(mnpOffset), std::get<1>(mnpOffset), std::get<2>(mnpOffset)};
244 indexedIntHKL = candidateIntHKL;
245 indexedMNP = candidateMNP;
246 foundSatellite =
true;
252 return std::make_tuple(indexedIntHKL, indexedMNP, indexedIntHKL.hklError());
268CombinedIndexingStats indexPeaks(
const std::vector<IPeak *> &peaks,
DblMatrix ub,
const double mainTolerance,
269 const bool roundHKLs,
const bool optimizeUB,
270 const Prop::SatelliteIndexingArgs &satelliteArgs) {
271 const auto nPeaks = peaks.size();
272 std::vector<V3D> qSample(nPeaks);
273 std::generate(std::begin(qSample), std::end(qSample),
274 [&peaks, i = 0u]()
mutable {
return peaks[i++]->getQSampleFrame(); });
279 CombinedIndexingStats stats;
281 for (
auto i = 0u; i < peaks.size(); ++i) {
282 const auto peak = peaks[i];
285 stats.main.numIndexed++;
286 stats.main.error += nominalHKL.hklError() / 3.0;
290 peak->setHKL(nominalHKL);
291 peak->setIntHKL(nominalHKL);
292 peak->setIntMNP(V3D(0, 0, 0));
293 }
else if (satelliteArgs.maxOrder > 0) {
294 auto result = indexSatellite(nominalHKL, satelliteArgs.maxOrder, satelliteArgs.modVectors,
295 satelliteArgs.tolerance, satelliteArgs.crossTerms);
297 const auto &satelliteInfo = result.get();
300 peak->setHKL(nominalHKL);
301 peak->setIntHKL(std::get<0>(satelliteInfo));
302 peak->setIntMNP(std::get<1>(satelliteInfo));
303 stats.satellites.numIndexed++;
304 stats.satellites.error += std::get<2>(satelliteInfo) / 3.;
308 peak->setHKL(V3D(0, 0, 0));
309 peak->setIntHKL(V3D(0, 0, 0));
310 peak->setIntMNP(V3D(0, 0, 0));
313 peak->setHKL(V3D(0, 0, 0));
314 peak->setIntHKL(V3D(0, 0, 0));
315 peak->setIntMNP(V3D(0, 0, 0));
330void logIndexingResults(std::ostream &out,
const CombinedIndexingStats &indexingInfo,
const int runNo,
331 const size_t nPeaksTotal,
const Prop::IndexPeaksArgs &args) {
333 out <<
"Run " << runNo;
336 out <<
" indexed " << indexingInfo.totalNumIndexed() <<
" peaks out of " << nPeaksTotal;
337 if (args.satellites.maxOrder > 0) {
338 out <<
" of which, " << indexingInfo.main.numIndexed <<
" main Bragg peaks are indexed with tolerance of "
339 << args.mainTolerance <<
", " << indexingInfo.satellites.numIndexed
340 <<
" satellite peaks are indexed with tolerance of " << args.satellites.tolerance <<
'\n';
341 out <<
" Average error in h,k,l for indexed peaks = " << indexingInfo.averageError() <<
'\n';
342 out <<
" Average error in h,k,l for indexed main peaks = " << indexingInfo.main.error <<
'\n';
343 out <<
" Average error in h,k,l for indexed satellite peaks = " << indexingInfo.satellites.error <<
'\n';
345 out <<
" with tolerance of " << args.mainTolerance <<
'\n';
346 out <<
" Average error in h,k,l for indexed peaks = " << indexingInfo.mainError() <<
'\n';
364 "Input Peaks Workspace");
366 auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
367 mustBePositiveDbl->setLower(0.0);
369 this->
declareProperty(Prop::SATE_TOLERANCE, 0.15, mustBePositiveDbl,
"Satellite peak indexing tolerance",
371 this->
declareProperty(Prop::ROUNDHKLS,
true,
"Round H, K and L values to integers");
372 this->
declareProperty(Prop::COMMONUB,
false,
"Index all orientations with a common UB");
375 "If true, update the OrientedLattice with the maxOrder, "
376 "modulation vectors & cross terms values input to the algorithm");
381 this->
declareProperty(Prop::MAIN_NUM_INDEXED, 0,
"Gets set with the number of indexed main peaks.",
383 this->
declareProperty(Prop::SATE_NUM_INDEXED, 0,
"Gets set with the number of indexed main peaks.",
385 this->
declareProperty(Prop::MAIN_ERR, 0.0,
"Gets set with the average HKL indexing error of Main Peaks.",
387 this->
declareProperty(Prop::SATE_ERR, 0.0,
"Gets set with the average HKL indexing error of Satellite Peaks.",
396 std::map<std::string, std::string> helpMsgs;
400 ws->sample().getOrientedLattice();
401 }
catch (std::runtime_error &) {
402 helpMsgs[Prop::PEAKSWORKSPACE] =
"No UB Matrix defined in the lattice.";
406 const auto args = Prop::IndexPeaksArgs::parse(*
this);
407 const bool isSave = this->
getProperty(Prop::SAVEMODINFO);
408 const bool isMOZero = (args.satellites.maxOrder == 0);
409 bool isAllVecZero =
true;
412 for (
size_t vecNo = 0; vecNo < args.satellites.modVectors.size(); vecNo++) {
413 if (args.satellites.modVectors[vecNo] !=
V3D(0.0, 0.0, 0.0)) {
414 isAllVecZero =
false;
416 g_log.
warning() <<
"Mod vector " << vecNo <<
" is invalid (0, 0, 0)"
420 if (isMOZero && !isAllVecZero) {
421 helpMsgs[
"MaxOrder"] =
"Max Order cannot be zero if a Modulation Vector has been supplied.";
423 if (!isMOZero && isAllVecZero) {
424 helpMsgs[
"ModVector1"] =
"At least one Modulation Vector must be supplied if Max Order set.";
426 if (isSave && isAllVecZero) {
427 helpMsgs[Prop::SAVEMODINFO] =
"Modulation info cannot be saved with no "
428 "valid Modulation Vectors supplied.";
430 if (isSave && isMOZero) {
431 helpMsgs[
"MaxOrder"] =
"Modulation info cannot be saved with Max Order = 0.";
440 const auto args = Prop::IndexPeaksArgs::parse(*
this);
443 if (args.workspace->getNumberPeaks() == 0) {
444 g_log.
warning(
"Empty peaks workspace. Nothing to index");
449 if (args.storeModulationInfo) {
450 auto &lattice = args.workspace->mutableSample().getOrientedLattice();
451 lattice.setMaxOrder(args.satellites.maxOrder);
452 lattice.setCrossTerm(args.satellites.crossTerms);
454 if (args.satellites.modVectors.size() >= 1) {
455 lattice.setModVec1(args.satellites.modVectors[0]);
458 lattice.setModVec1(
V3D(0.0, 0.0, 0.0));
461 if (args.satellites.modVectors.size() >= 2) {
462 lattice.setModVec2(args.satellites.modVectors[1]);
465 lattice.setModVec2(
V3D(0.0, 0.0, 0.0));
468 if (args.satellites.modVectors.size() >= 3) {
469 lattice.setModVec3(args.satellites.modVectors[2]);
472 lattice.setModVec3(
V3D(0.0, 0.0, 0.0));
476 CombinedIndexingStats indexingInfo;
477 const auto &lattice = args.workspace->sample().getOrientedLattice();
478 const auto &sampleUB = lattice.getUB();
481 std::vector<IPeak *> allPeaksRef;
482 allPeaksRef.reserve(args.workspace->getNumberPeaks());
483 for (
int i = 0; i < args.workspace->getNumberPeaks(); i++) {
484 allPeaksRef.emplace_back(args.workspace->getPeakPtr(i));
486 const bool optimizeUB{
false};
487 indexingInfo = indexPeaks(allPeaksRef, sampleUB, args.mainTolerance, args.roundHKLs, optimizeUB, args.satellites);
490 std::unordered_map<int, std::vector<IPeak *>> peaksPerRun;
491 for (
int i = 0; i < args.workspace->getNumberPeaks(); i++)
492 peaksPerRun[args.workspace->getPeak(i).getRunNumber()].emplace_back(args.workspace->getPeakPtr(i));
493 if (peaksPerRun.size() < 2) {
494 g_log.
warning(
"Peaks from only one run exist but CommonUBForAll=True so peaks will be indexed with an optimised "
495 "UB which will not be saved in the workspace.");
497 const bool optimizeUB{
true};
498 for (
const auto &runPeaks : peaksPerRun) {
499 const auto &peaks = runPeaks.second;
500 const auto indexedInRun =
501 indexPeaks(peaks, sampleUB, args.mainTolerance, args.roundHKLs, optimizeUB, args.satellites);
502 logIndexingResults(
g_log.
notice(), indexedInRun, runPeaks.first, peaks.size(), args);
503 indexingInfo += indexedInRun;
507 setProperty(
"NumIndexed", indexingInfo.totalNumIndexed());
508 setProperty(
"MainNumIndexed", indexingInfo.main.numIndexed);
509 setProperty(
"SateNumIndexed", indexingInfo.satellites.numIndexed);
510 setProperty(
"AverageError", indexingInfo.averageError());
511 setProperty(
"MainError", indexingInfo.mainError());
512 setProperty(
"SatelliteError", indexingInfo.satelliteError());
515 logIndexingResults(
g_log.
notice(), indexingInfo, -1, args.workspace->getNumberPeaks(), args);
517 g_log.
notice() << args.workspace->sample().getOrientedLattice() <<
"\n";
#define DECLARE_ALGORITHM(classname)
const std::vector< double > & rhs
const int64_t TOLERANCE(1000000)
const std::vector< V3D > modVectors
const bool storeModulationInfo
const SatelliteIndexingArgs satellites
IPeaksWorkspace_sptr workspace
const double mainTolerance
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
A property class for workspaces.
void init() override
Initialize the algorithm's properties.
std::map< std::string, std::string > validateInputs() override
Validate all inputs once set.
void exec() override
Execute the algorithm.
static void RoundHKL(Kernel::V3D &hkl)
Round all the components of a HKL objects to the nearest integer.
static void RoundHKLs(std::vector< Kernel::V3D > &hkl_list)
Round all the components of a list of V3D objects, to the nearest integer.
static double Optimize_UB(Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &hkl_vectors, const std::vector< Kernel::V3D > &q_vectors, std::vector< double > &sigabc)
Find the UB matrix that most nearly maps hkl to qxyz for 3 or more peaks.
static bool ValidIndex(const Kernel::V3D &hkl, double tolerance)
Check is hkl is within tolerance of integer (h,k,l) non-zero values.
static int CalculateMillerIndices(const Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &q_vectors, double tolerance, std::vector< Kernel::V3D > &miller_indices, double &ave_error)
Given a UB, get list of Miller indices for specifed Qs and tolerance.
ArrayLenghtValidator : Validate length of an array property.
Support for a property that holds an array of values.
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
void warning(const std::string &msg)
Logs at warning level.
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
MatrixWorkspace_sptr MANTID_API_DLL operator+=(const MatrixWorkspace_sptr &lhs, const MatrixWorkspace_sptr &rhs)
Adds two workspaces.
std::vector< Kernel::V3D > addModulationVectors(const std::vector< double > &modVector1, const std::vector< double > &modVector2, const std::vector< double > &modVector3)
Create a list of valid modulation vectors from the input.
std::vector< MNPOffset > generateOffsetVectors(const std::vector< Kernel::V3D > &modVectors, const int maxOrder, const bool crossTerms)
Calculate a list of HKL offsets from the given modulation vectors.
std::vector< Kernel::V3D > validModulationVectors(const std::vector< double > &modVector1, const std::vector< double > &modVector2, const std::vector< double > &modVector3)
Create a list of valid modulation vectors from the input.
Mantid::Kernel::Matrix< double > DblMatrix
String constants for algorithm's properties.
static const std::string ModVector3
static const std::string MaxOrder
static const std::string ModVector2
static const std::string CrossTerms
static const std::string ModVector1
static void appendTo(API::IAlgorithm *alg)
Append the common set of properties that relate to modulation vectors to the given algorithm.
Describes the direction (within an algorithm) of a Property.
@ InOut
Both an input & output workspace.
@ Input
An input workspace.
@ Output
An output workspace.