PeakStatisticsTools.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/PeakStatisticsTools.h"
 
#include "MantidGeometry/Crystal/BasicHKLFilters.h"
#include "MantidGeometry/Crystal/HKLGenerator.h"
 
#include "MantidKernel/Statistics.h"
 
#include <memory>
#include <numeric>
#include <utility>
 
namespace Mantid::Crystal::PeakStatisticsTools {
 
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
 
std::vector<double> UniqueReflection::getWavelengths() const {
  std::vector<double> wavelengths;
  wavelengths.reserve(m_peaks.size());
 
  std::transform(m_peaks.begin(), m_peaks.end(), std::back_inserter(wavelengths),
                 [](const DataObjects::Peak &peak) { return peak.getWavelength(); });
 
  return wavelengths;
}
 
std::vector<double> UniqueReflection::getIntensities() const {
  std::vector<double> intensities;
  intensities.reserve(m_peaks.size());
 
  std::transform(m_peaks.begin(), m_peaks.end(), std::back_inserter(intensities),
                 [](const DataObjects::Peak &peak) { return peak.getIntensity(); });
 
  return intensities;
}
 
std::vector<double> UniqueReflection::getSigmas() const {
  std::vector<double> sigmas;
  sigmas.reserve(m_peaks.size());
 
  std::transform(m_peaks.begin(), m_peaks.end(), std::back_inserter(sigmas),
                 [](const DataObjects::Peak &peak) { return peak.getSigmaIntensity(); });
 
  return sigmas;
}
 
UniqueReflection UniqueReflection::removeOutliers(double sigmaCritical, bool weightedZ) const {
  if (sigmaCritical <= 0.0) {
    throw std::invalid_argument("Critical sigma value has to be greater than 0.");
  }
 
  UniqueReflection newReflection(m_hkl);
 
  if (m_peaks.size() > 2) {
    auto intensities = getIntensities();
    std::vector<double> zScores;
    if (!weightedZ) {
      zScores = Kernel::getZscore(intensities);
    } else {
      auto sigmas = getSigmas();
      zScores = Kernel::getWeightedZscore(intensities, sigmas);
    }
 
    for (size_t i = 0; i < zScores.size(); ++i) {
      if (zScores[i] <= sigmaCritical) {
        newReflection.addPeak(m_peaks[i]);
      }
    }
  } else {
    for (const auto &peak : m_peaks) {
      newReflection.addPeak(peak);
    }
  }
 
  return newReflection;
}
 
void UniqueReflection::setPeaksIntensityAndSigma(double intensity, double sigma) {
  for (auto &peak : m_peaks) {
    peak.setIntensity(intensity);
    peak.setSigmaIntensity(sigma);
  }
}
 
UniqueReflectionCollection::UniqueReflectionCollection(const UnitCell &cell, const std::pair<double, double> &dLimits,
                                                       PointGroup_sptr pointGroup,
                                                       const ReflectionCondition_sptr &centering)
    : m_reflections(), m_pointgroup(std::move(pointGroup)) {
  HKLGenerator generator(cell, dLimits.first);
  auto dFilter = std::make_shared<const HKLFilterDRange>(cell, dLimits.first, dLimits.second);
  auto centeringFilter = std::make_shared<const HKLFilterCentering>(centering);
  auto filter = dFilter & centeringFilter;
 
  // Generate map of UniqueReflection-objects with reflection family as key.
  for (const auto &hkl : generator) {
    if (filter->isAllowed(hkl)) {
      V3D hklFamily = m_pointgroup->getReflectionFamily(hkl);
      m_reflections.emplace(hklFamily, UniqueReflection(hklFamily));
    }
  }
}
 
void UniqueReflectionCollection::addObservations(const std::vector<Peak> &peaks) {
  for (auto const &peak : peaks) {
    V3D hkl = peak.getHKL();
    hkl.round();
 
    auto reflection = m_reflections.find(m_pointgroup->getReflectionFamily(hkl));
 
    if (reflection != m_reflections.end()) {
      (*reflection).second.addPeak(peak);
    }
  }
}
 
UniqueReflection UniqueReflectionCollection::getReflection(const V3D &hkl) const {
  return m_reflections.at(m_pointgroup->getReflectionFamily(hkl));
}
 
size_t UniqueReflectionCollection::getUniqueReflectionCount() const { return m_reflections.size(); }
 
size_t UniqueReflectionCollection::getObservedUniqueReflectionCount(size_t moreThan) const {
  return std::count_if(
      m_reflections.cbegin(), m_reflections.cend(),
      [=](const std::pair<Kernel::V3D, UniqueReflection> &item) { return item.second.count() > moreThan; });
}
 
std::vector<V3D> UniqueReflectionCollection::getUnobservedUniqueReflections() const {
  std::vector<V3D> reflections;
  reflections.reserve(m_reflections.size());
 
  for (const auto &reflection : m_reflections) {
    if (reflection.second.count() == 0) {
      reflections.emplace_back(reflection.first);
    }
  }
 
  return reflections;
}
 
size_t UniqueReflectionCollection::getObservedReflectionCount() const {
  return std::accumulate(m_reflections.cbegin(), m_reflections.cend(), size_t(0),
                         [](size_t totalReflections, const std::pair<Kernel::V3D, UniqueReflection> &item) {
                           return totalReflections + item.second.count();
                         });
}
 
const std::map<V3D, UniqueReflection> &UniqueReflectionCollection::getReflections() const { return m_reflections; }
 
void PeaksStatistics::calculatePeaksStatistics(const std::map<V3D, UniqueReflection> &uniqueReflections,
                                               std::string &equivalentIntensities, double &sigmaCritical,
                                               bool &weightedZ) {
  double rMergeNumerator = 0.0;
  double rPimNumerator = 0.0;
  double intensitySumRValues = 0.0;
  double iOverSigmaSum = 0.0;
 
  for (const auto &unique : uniqueReflections) {
    /* Since all possible unique reflections are explored
     * there may be 0 observations for some of them.
     * In that case, nothing can be done.*/
    if (unique.second.count() > 0) {
      ++m_uniqueReflections;
 
      // Possibly remove outliers.
      auto outliersRemoved = unique.second.removeOutliers(sigmaCritical, weightedZ);
 
      // I/sigma is calculated for all reflections, even if there is only one
      // observation.
      auto intensities = outliersRemoved.getIntensities();
      auto sigmas = outliersRemoved.getSigmas();
 
      // Accumulate the I/sigma's for current reflection into sum
      iOverSigmaSum += getIOverSigmaSum(sigmas, intensities);
 
      if (outliersRemoved.count() > 1) {
        // Get mean, standard deviation for intensities
        auto intensityStatistics = Kernel::getStatistics(
            intensities, StatOptions::Mean | StatOptions::UncorrectedStdDev | StatOptions::Median);
 
        double meanIntensity = intensityStatistics.mean;
        if (equivalentIntensities == "Median")
          meanIntensity = intensityStatistics.median;
 
        /* This was in the original algorithm, not entirely sure where it is
         * used. It's basically the sum of all relative standard deviations.
         * In a perfect data set with all equivalent reflections exactly
         * equivalent that would be 0. */
        m_chiSquared += intensityStatistics.standard_deviation / meanIntensity;
 
        // For both RMerge and RPim sum(|I - <I>|) is required
        double sumOfDeviationsFromMean =
            std::accumulate(intensities.begin(), intensities.end(), 0.0, [meanIntensity](double sum, double intensity) {
              return sum + fabs(intensity - meanIntensity);
            });
 
        // Accumulate into total sum for numerator of RMerge
        rMergeNumerator += sumOfDeviationsFromMean;
 
        // For Rpim, the sum is weighted by a factor depending on N
        double rPimFactor = sqrt(1.0 / (static_cast<double>(outliersRemoved.count()) - 1.0));
        rPimNumerator += (rPimFactor * sumOfDeviationsFromMean);
 
        // Collect sum of intensities for R-value calculation
        intensitySumRValues += std::accumulate(intensities.begin(), intensities.end(), 0.0);
      }
 
      const std::vector<Peak> &reflectionPeaks = outliersRemoved.getPeaks();
      m_peaks.insert(m_peaks.end(), reflectionPeaks.begin(), reflectionPeaks.end());
    }
  }
 
  m_measuredReflections = static_cast<int>(m_peaks.size());
 
  if (m_uniqueReflections > 0) {
    m_redundancy = static_cast<double>(m_measuredReflections) / static_cast<double>(m_uniqueReflections);
  }
 
  m_completeness = static_cast<double>(m_uniqueReflections) / static_cast<double>(uniqueReflections.size());
 
  if (intensitySumRValues > 0.0) {
    m_rMerge = rMergeNumerator / intensitySumRValues;
    m_rPim = rPimNumerator / intensitySumRValues;
  }
 
  if (m_measuredReflections > 0) {
    m_meanIOverSigma = iOverSigmaSum / static_cast<double>(m_measuredReflections);
 
    auto dspacingLimits = getDSpacingLimits(m_peaks);
    m_dspacingMin = dspacingLimits.first;
    m_dspacingMax = dspacingLimits.second;
  }
}
 
double PeaksStatistics::getIOverSigmaSum(const std::vector<double> &sigmas,
                                         const std::vector<double> &intensities) const {
  return std::inner_product(intensities.begin(), intensities.end(), sigmas.begin(), 0.0, std::plus<double>(),
                            std::divides<double>());
}
 
double PeaksStatistics::getRMS(const std::vector<double> &data) const {
  double sumOfSquares = std::inner_product(data.begin(), data.end(), data.begin(), 0.0);
 
  return sqrt(sumOfSquares / static_cast<double>(data.size()));
}
 
std::pair<double, double> PeaksStatistics::getDSpacingLimits(const std::vector<Peak> &peaks) const {
  if (peaks.empty()) {
    return std::make_pair(0.0, 0.0);
  }
 
  auto dspacingLimitIterators = std::minmax_element(peaks.begin(), peaks.end(), [](const Peak &lhs, const Peak &rhs) {
    return lhs.getDSpacing() < rhs.getDSpacing();
  });
 
  return std::make_pair((*(dspacingLimitIterators.first)).getDSpacing(),
                        (*(dspacingLimitIterators.second)).getDSpacing());
}
 
} // namespace Mantid::Crystal::PeakStatisticsTools