IndexPeaks.cpp

Go to the documentation of this file.

// 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 +
#include "MantidCrystal/IndexPeaks.h"
#include "MantidAPI/IPeaksWorkspace.h"
#include "MantidAPI/Sample.h"
#include "MantidCrystal/PeakAlgorithmHelpers.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
 
#include <optional>
 
#include <algorithm>
 
namespace Mantid::Crystal {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(IndexPeaks)
 
using API::IPeaksWorkspace;
using API::IPeaksWorkspace_sptr;
using Geometry::IndexingUtils;
using Geometry::IPeak;
using Geometry::OrientedLattice;
using Kernel::DblMatrix;
using Kernel::Logger;
using Kernel::V3D;
 
namespace {
const auto OPTIMIZE_UB_ATTEMPTS{4};
 
namespace Prop {
const std::string PEAKSWORKSPACE{"PeaksWorkspace"};
const std::string TOLERANCE{"Tolerance"};
const std::string SATE_TOLERANCE{"ToleranceForSatellite"};
const std::string ROUNDHKLS{"RoundHKLs"};
const std::string COMMONUB{"CommonUBForAll"};
const std::string SAVEMODINFO{"SaveModulationInfo"};
const std::string AVERAGE_ERR{"AverageError"};
const std::string NUM_INDEXED{"NumIndexed"};
const std::string MAIN_NUM_INDEXED{"MainNumIndexed"};
const std::string SATE_NUM_INDEXED{"SateNumIndexed"};
const std::string MAIN_ERR{"MainError"};
const std::string SATE_ERR{"SatelliteError"};
const std::string UPDATE_UB{"UpdateUB"};
 
struct SatelliteIndexingArgs {
  const double tolerance;
  const int maxOrder;
  const std::vector<V3D> modVectors;
  const bool crossTerms;
};
 
struct IndexPeaksArgs {
  static IndexPeaksArgs parse(const API::Algorithm &alg) {
    const IPeaksWorkspace_sptr peaksWS = alg.getProperty(PEAKSWORKSPACE);
    const int maxOrderFromAlg = alg.getProperty(ModulationProperties::MaxOrder);
 
    // Init variables
    int maxOrderToUse{0};
    std::vector<V3D> modVectorsToUse;
    modVectorsToUse.reserve(3);
    bool crossTermToUse{false};
 
    // Parse mod vectors
    modVectorsToUse = addModulationVectors(alg.getProperty(ModulationProperties::ModVector1),
                                           alg.getProperty(ModulationProperties::ModVector2),
                                           alg.getProperty(ModulationProperties::ModVector3));
    // check the 3 mod vectors added from properties
    modVectorsToUse = validModulationVectors(modVectorsToUse[0], modVectorsToUse[1], modVectorsToUse[2]);
 
    // deal with case: max order > 0 and no mod vector is specified
    if (maxOrderFromAlg > 0 && modVectorsToUse.size() == 0) {
      // Max Order is larger than zero but no modulated vector specified
      // Assume that the caller method will handle this
      maxOrderToUse = maxOrderFromAlg;
    } else if (maxOrderFromAlg == 0 && modVectorsToUse.size() == 0) {
      // Use lattice definitions if they exist
      const auto &lattice = peaksWS->sample().getOrientedLattice();
      crossTermToUse = lattice.getCrossTerm();
      maxOrderToUse = lattice.getMaxOrder(); // the lattice can return a 0 here
 
      // if lattice has maxOrder, we will use the modVec from it, otherwise
      // stick to the input got from previous assignment
      if (maxOrderToUse > 0) {
        modVectorsToUse = validModulationVectors(lattice.getModVec(0), lattice.getModVec(1), lattice.getModVec(2));
      }
    } else {
      // Use user specified
      // default behavior: map everything automatically
      maxOrderToUse = maxOrderFromAlg;
      crossTermToUse = alg.getProperty(ModulationProperties::CrossTerms);
    }
 
    return {peaksWS,
            alg.getProperty(TOLERANCE),
            alg.getProperty(ROUNDHKLS),
            alg.getProperty(COMMONUB),
            alg.getProperty(SAVEMODINFO),
            SatelliteIndexingArgs{alg.getProperty(SATE_TOLERANCE), maxOrderToUse, modVectorsToUse, crossTermToUse},
            alg.getProperty(UPDATE_UB)};
  }
 
  IPeaksWorkspace_sptr workspace;
  const double mainTolerance;
  const bool roundHKLs;
  const bool commonUB;
  const bool storeModulationInfo;
  const SatelliteIndexingArgs satellites;
  const bool updateUB;
};
} // namespace Prop
 
struct PeakIndexingStats {
  PeakIndexingStats &operator+=(const PeakIndexingStats &rhs) {
    numIndexed += rhs.numIndexed;
    error += rhs.error;
    return *this;
  }
  int numIndexed{0};
  double error{0.0};
};
 
struct CombinedIndexingStats {
  CombinedIndexingStats &operator+=(const CombinedIndexingStats &rhs) {
    main += rhs.main;
    satellites += rhs.satellites;
    return *this;
  }
 
  inline int totalNumIndexed() const { return main.numIndexed + satellites.numIndexed; }
  inline double mainError() const {
    if (main.numIndexed == 0)
      return 0.0;
    return main.error / main.numIndexed;
  }
  inline double satelliteError() const {
    if (satellites.numIndexed == 0)
      return 0.0;
    return satellites.error / satellites.numIndexed;
  }
  inline double averageError() const { return (main.error + satellites.error) / totalNumIndexed(); }
 
  PeakIndexingStats main;
  PeakIndexingStats satellites;
};
 
DblMatrix optimizeUBMatrix(const DblMatrix &ubOrig, const std::vector<V3D> &qSample, const double tolerance) {
  DblMatrix optimizedUB(ubOrig);
 
  double errorAtStart{0.0};
  std::vector<V3D> millerIndices;
  millerIndices.reserve(qSample.size());
  const int numIndexedAtStart =
      IndexingUtils::CalculateMillerIndices(optimizedUB, qSample, tolerance, millerIndices, errorAtStart);
 
  if (numIndexedAtStart < 3) {
    // can't optimize without at least 3 indexed peaks
    return optimizedUB;
  }
 
  for (auto i = 0; i < OPTIMIZE_UB_ATTEMPTS; ++i) {
    try {
      // optimization requires rounded indices
      IndexingUtils::RoundHKLs(millerIndices);
      IndexingUtils::Optimize_UB(optimizedUB, millerIndices, qSample);
    } catch (...) {
      // If there is any problem, such as too few
      // independent peaks, just use the original UB
      optimizedUB = ubOrig;
      break;
    }
    double errorInLoop{0.0};
    const int numIndexedInLoop =
        IndexingUtils::CalculateMillerIndices(optimizedUB, qSample, tolerance, millerIndices, errorInLoop);
    if (numIndexedInLoop < numIndexedAtStart) // use the original UB
      break;
  }
  return optimizedUB;
}
 
using IndexedSatelliteInfo = std::tuple<V3D, V3D, V3D, double>;
 
std::optional<IndexedSatelliteInfo> indexSatellite(const V3D &mainHKL, const int maxOrder,
                                                   const std::vector<V3D> &modVectors, const double tolerance,
                                                   const bool crossTerms) {
  const auto offsets = generateOffsetVectors(modVectors, maxOrder, crossTerms);
  bool foundSatellite{false};
  V3D indexedIntHKL, indexedMNP, fractionalOffset;
  for (const auto &mnpOffset : offsets) {
    const auto candidateIntHKL = mainHKL - std::get<3>(mnpOffset);
    const V3D candidateMNP{std::get<0>(mnpOffset), std::get<1>(mnpOffset), std::get<2>(mnpOffset)};
    if (IndexingUtils::ValidIndex(candidateIntHKL, tolerance)) {
      indexedIntHKL = candidateIntHKL;
      indexedMNP = candidateMNP;
      fractionalOffset = std::get<3>(mnpOffset);
      foundSatellite = true;
      // we deliberately don't break and use the last valid
      // reflection we find.
    }
  }
  if (foundSatellite)
    return std::make_tuple(fractionalOffset, indexedIntHKL, indexedMNP, indexedIntHKL.hklError());
  else
    return std::nullopt;
}
 
CombinedIndexingStats indexPeaks(const std::vector<IPeak *> &peaks, DblMatrix ub, const double mainTolerance,
                                 const bool roundHKLs, const bool optimizeUB,
                                 const Prop::SatelliteIndexingArgs &satelliteArgs) {
  const auto nPeaks = peaks.size();
  std::vector<V3D> qSample(nPeaks);
  std::generate(std::begin(qSample), std::end(qSample),
                [&peaks, i = 0u]() mutable { return peaks[i++]->getQSampleFrame(); });
  if (optimizeUB) {
    ub = optimizeUBMatrix(ub, qSample, mainTolerance);
  }
 
  CombinedIndexingStats stats;
  ub.Invert();
  for (auto i = 0u; i < peaks.size(); ++i) {
    const auto peak = peaks[i];
    V3D nominalHKL = IndexingUtils::CalculateMillerIndices(ub, qSample[i]);
    if (IndexingUtils::ValidIndex(nominalHKL, mainTolerance)) {
      stats.main.numIndexed++;
      stats.main.error += nominalHKL.hklError() / 3.0;
      if (roundHKLs) {
        IndexingUtils::RoundHKL(nominalHKL);
      }
      peak->setHKL(nominalHKL);
      peak->setIntHKL(nominalHKL);
      peak->setIntMNP(V3D(0, 0, 0));
    } else if (satelliteArgs.maxOrder > 0) {
      auto result = indexSatellite(nominalHKL, satelliteArgs.maxOrder, satelliteArgs.modVectors,
                                   satelliteArgs.tolerance, satelliteArgs.crossTerms);
      if (result) {
        const auto &satelliteInfo = result.value();
        V3D hkl;
        ;
        if (roundHKLs) {
          hkl = std::get<1>(satelliteInfo);
          IndexingUtils::RoundHKL(hkl);
          hkl += std::get<0>(satelliteInfo);
        } else {
          hkl = nominalHKL;
        }
        peak->setHKL(hkl);
        peak->setIntHKL(std::get<1>(satelliteInfo));
        peak->setIntMNP(std::get<2>(satelliteInfo));
        stats.satellites.numIndexed++;
        stats.satellites.error += std::get<3>(satelliteInfo) / 3.;
      } else {
        // clear these to make sure leftover values from previous index peaks
        // run are not used
        peak->setHKL(V3D(0, 0, 0));
        peak->setIntHKL(V3D(0, 0, 0));
        peak->setIntMNP(V3D(0, 0, 0));
      }
    } else {
      peak->setHKL(V3D(0, 0, 0));
      peak->setIntHKL(V3D(0, 0, 0));
      peak->setIntMNP(V3D(0, 0, 0));
    }
  }
  return stats;
}
 
void logIndexingResults(std::ostream &out, const CombinedIndexingStats &indexingInfo, const int runNo,
                        const size_t nPeaksTotal, const Prop::IndexPeaksArgs &args) {
  if (runNo >= 0)
    out << "Run " << runNo;
  else
    out << "All runs";
  out << " indexed " << indexingInfo.totalNumIndexed() << " peaks out of " << nPeaksTotal;
  if (args.satellites.maxOrder > 0) {
    out << " of which, " << indexingInfo.main.numIndexed << " main Bragg peaks are indexed with tolerance of "
        << args.mainTolerance << ", " << indexingInfo.satellites.numIndexed
        << " satellite peaks are indexed with tolerance of " << args.satellites.tolerance << '\n';
    out << "  Average error in h,k,l for indexed peaks =  " << indexingInfo.averageError() << '\n';
    out << "  Average error in h,k,l for indexed main peaks =  " << indexingInfo.main.error << '\n';
    out << "  Average error in h,k,l for indexed satellite peaks =  " << indexingInfo.satellites.error << '\n';
  } else {
    out << " with tolerance of " << args.mainTolerance << '\n';
    out << "  Average error in h,k,l for indexed peaks =  " << indexingInfo.mainError() << '\n';
  }
}
 
} // namespace
 
void IndexPeaks::init() {
  using API::WorkspaceProperty;
  using Kernel::ArrayLengthValidator;
  using Kernel::ArrayProperty;
  using Kernel::BoundedValidator;
  using Kernel::Direction;
 
  // -- inputs --
  this->declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace_sptr::element_type>>(Prop::PEAKSWORKSPACE,
                                                                                                "", Direction::InOut),
                        "Input Peaks Workspace");
 
  auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
  mustBePositiveDbl->setLower(0.0);
  this->declareProperty(Prop::TOLERANCE, 0.15, mustBePositiveDbl, "Main peak indexing tolerance", Direction::Input);
  this->declareProperty(Prop::SATE_TOLERANCE, 0.15, mustBePositiveDbl, "Satellite peak indexing tolerance",
                        Direction::Input);
  this->declareProperty(Prop::ROUNDHKLS, true, "Round H, K and L values to integers");
  this->declareProperty(Prop::COMMONUB, false, "Index all orientations with a common UB");
  ModulationProperties::appendTo(this);
  this->declareProperty(Prop::SAVEMODINFO, false,
                        "If true, update the OrientedLattice with the maxOrder, "
                        "modulation vectors & cross terms values input to the algorithm");
  this->declareProperty(Prop::UPDATE_UB, false,
                        "Saves the optimized UB matrix to the workspace if CommonUBForAll=False. This option works "
                        "only when the peak workspace contains a single run.");
 
  // -- outputs --
  this->declareProperty(Prop::NUM_INDEXED, 0, "Gets set with the number of indexed peaks.", Direction::Output);
  this->declareProperty(Prop::AVERAGE_ERR, 0.0, "Gets set with the average HKL indexing error.", Direction::Output);
  this->declareProperty(Prop::MAIN_NUM_INDEXED, 0, "Gets set with the number of indexed main peaks.",
                        Direction::Output);
  this->declareProperty(Prop::SATE_NUM_INDEXED, 0, "Gets set with the number of indexed main peaks.",
                        Direction::Output);
  this->declareProperty(Prop::MAIN_ERR, 0.0, "Gets set with the average HKL indexing error of Main Peaks.",
                        Direction::Output);
  this->declareProperty(Prop::SATE_ERR, 0.0, "Gets set with the average HKL indexing error of Satellite Peaks.",
                        Direction::Output);
}
 
std::map<std::string, std::string> IndexPeaks::validateInputs() {
  std::map<std::string, std::string> helpMsgs;
 
  IPeaksWorkspace_sptr ws = this->getProperty(Prop::PEAKSWORKSPACE);
  try {
    ws->sample().getOrientedLattice();
  } catch (std::runtime_error &) {
    helpMsgs[Prop::PEAKSWORKSPACE] = "No UB Matrix defined in the lattice.";
    return helpMsgs;
  }
 
  const auto args = Prop::IndexPeaksArgs::parse(*this);
  const bool isSave = this->getProperty(Prop::SAVEMODINFO);
  const bool isMOZero = (args.satellites.maxOrder == 0);
  bool isAllVecZero = true;
  // parse() validates all the mod vectors. There should not be any modulated
  // vector in modVectors is equal to (0, 0, 0)
  for (size_t vecNo = 0; vecNo < args.satellites.modVectors.size(); vecNo++) {
    if (args.satellites.modVectors[vecNo] != V3D(0.0, 0.0, 0.0)) {
      isAllVecZero = false;
    } else {
      g_log.warning() << "Mod vector " << vecNo << " is invalid (0, 0, 0)"
                      << "\n";
    }
  }
  if (isMOZero && !isAllVecZero) {
    helpMsgs["MaxOrder"] = "Max Order cannot be zero if a Modulation Vector has been supplied.";
  }
  if (!isMOZero && isAllVecZero) {
    helpMsgs["ModVector1"] = "At least one Modulation Vector must be supplied if Max Order set.";
  }
  if (isSave && isAllVecZero) {
    helpMsgs[Prop::SAVEMODINFO] = "Modulation info cannot be saved with no "
                                  "valid Modulation Vectors supplied.";
  }
  if (isSave && isMOZero) {
    helpMsgs["MaxOrder"] = "Modulation info cannot be saved with Max Order = 0.";
  }
  return helpMsgs;
}
 
void IndexPeaks::exec() {
  const auto args = Prop::IndexPeaksArgs::parse(*this);
 
  // quick exit
  if (args.workspace->getNumberPeaks() == 0) {
    g_log.warning("Empty peaks workspace. Nothing to index");
    return;
  }
 
  // save modulation input if asked
  if (args.storeModulationInfo) {
    auto &lattice = args.workspace->mutableSample().getOrientedLattice();
    lattice.setMaxOrder(args.satellites.maxOrder);
    lattice.setCrossTerm(args.satellites.crossTerms);
 
    if (args.satellites.modVectors.size() >= 1) {
      lattice.setModVec1(args.satellites.modVectors[0]);
    } else {
      g_log.warning("empty modVector 1, skipping saving");
      lattice.setModVec1(V3D(0.0, 0.0, 0.0));
    }
 
    if (args.satellites.modVectors.size() >= 2) {
      lattice.setModVec2(args.satellites.modVectors[1]);
    } else {
      g_log.warning("empty modVector 2, skipping saving");
      lattice.setModVec2(V3D(0.0, 0.0, 0.0));
    }
 
    if (args.satellites.modVectors.size() >= 3) {
      lattice.setModVec3(args.satellites.modVectors[2]);
    } else {
      g_log.warning("empty modVector 3, skipping saving");
      lattice.setModVec3(V3D(0.0, 0.0, 0.0));
    }
    // set modUB now mod vectors populated
    lattice.setModUB(lattice.getUB() * lattice.getModHKL());
  }
 
  CombinedIndexingStats indexingInfo;
  const auto &lattice = args.workspace->sample().getOrientedLattice();
  const auto &sampleUB = lattice.getUB();
  if (args.commonUB) {
    // Use sample UB an all peaks regardless of run
    std::vector<IPeak *> allPeaksRef;
    allPeaksRef.reserve(args.workspace->getNumberPeaks());
    for (int i = 0; i < args.workspace->getNumberPeaks(); i++) {
      allPeaksRef.emplace_back(args.workspace->getPeakPtr(i));
    }
    const bool optimizeUB{false};
    indexingInfo = indexPeaks(allPeaksRef, sampleUB, args.mainTolerance, args.roundHKLs, optimizeUB, args.satellites);
  } else {
    // Use a UB optimized for each run
    std::unordered_map<int, std::vector<IPeak *>> peaksPerRun;
    for (int i = 0; i < args.workspace->getNumberPeaks(); i++)
      peaksPerRun[args.workspace->getPeak(i).getRunNumber()].emplace_back(args.workspace->getPeakPtr(i));
    if (peaksPerRun.size() == 1) {
      if (args.updateUB) {
        // Save the optimized UB to the workspace
        auto peaks = peaksPerRun.begin()->second;
        std::vector<V3D> qSample(peaks.size());
        std::generate(std::begin(qSample), std::end(qSample),
                      [&peaks, i = 0u]() mutable { return peaks[i++]->getQSampleFrame(); });
 
        DblMatrix optimizedUB = optimizeUBMatrix(sampleUB, qSample, args.mainTolerance);
        if (optimizedUB != sampleUB) {
          args.workspace->mutableSample().getOrientedLattice().setUB(optimizedUB);
          g_log.notice() << "Updating the UB matrix with an improved optimized version.\n";
        } else {
          g_log.notice() << "No improved UB matrix was found, so retaining the original UB.\n";
        }
      } else {
        g_log.warning(
            "Peaks from only one run exist but CommonUBForAll=False so peaks will be indexed with an optimised "
            "UB which will not be saved in the workspace. Use UpdateUB=True to save the optimized UB matrix.\n");
      }
    }
    const bool optimizeUB{true};
    for (const auto &runPeaks : peaksPerRun) {
      const auto &peaks = runPeaks.second;
      const auto indexedInRun =
          indexPeaks(peaks, sampleUB, args.mainTolerance, args.roundHKLs, optimizeUB, args.satellites);
      logIndexingResults(g_log.notice(), indexedInRun, runPeaks.first, peaks.size(), args);
      indexingInfo += indexedInRun;
    }
  }
 
  setProperty("NumIndexed", indexingInfo.totalNumIndexed());
  setProperty("MainNumIndexed", indexingInfo.main.numIndexed);
  setProperty("SateNumIndexed", indexingInfo.satellites.numIndexed);
  setProperty("AverageError", indexingInfo.averageError());
  setProperty("MainError", indexingInfo.mainError());
  setProperty("SatelliteError", indexingInfo.satelliteError());
 
  // Final results
  logIndexingResults(g_log.notice(), indexingInfo, -1, args.workspace->getNumberPeaks(), args);
  // Show the lattice parameters
  g_log.notice() << args.workspace->sample().getOrientedLattice() << "\n";
}
 
} // namespace Mantid::Crystal