Refs #21752 make weightedZScore input option and test

Framework/Crystal/inc/MantidCrystal/LoadIsawPeaks.h

@@ -52,7 +52,7 @@ private:
   /// Read a single peak from peaks file
   DataObjects::Peak readPeak(DataObjects::PeaksWorkspace_sptr outWS,
                              std::string &lastStr, std::ifstream &in,
-                             int &seqNum, std::string bankName, double qSign);
+                             size_t &seqNum, std::string bankName, double qSign);
 
   int findPixelID(Geometry::Instrument_const_sptr inst, std::string bankName,
                   int col, int row);

Framework/Crystal/inc/MantidCrystal/PeakStatisticsTools.h

@@ -38,7 +38,7 @@ public:
   std::vector<double> getIntensities() const;
   std::vector<double> getSigmas() const;
 
-  UniqueReflection removeOutliers(double sigmaCritical = 3.0) const;
+  UniqueReflection removeOutliers(double sigmaCritical = 3.0, bool weightedZ = false) const;
   void setPeaksIntensityAndSigma(double intensity, double sigma);
 
 private:
@@ -105,15 +105,26 @@ private:
  */
 class DLLExport PeaksStatistics {
 public:
+  explicit PeaksStatistics(const UniqueReflectionCollection &reflections)
+      : m_measuredReflections(0), m_uniqueReflections(0), m_completeness(0.0),
+        m_redundancy(0.0), m_rMerge(0.0), m_rPim(0.0), m_meanIOverSigma(0.0),
+        m_dspacingMin(0.0), m_dspacingMax(0.0), m_chiSquared(0.0), m_peaks() {
+    m_peaks.reserve(reflections.getObservedReflectionCount());
+    std::string equivalentIntensities = "Mean";
+    double sigmaCritical = 3.0;
+    bool weightedZ = false;
+    calculatePeaksStatistics(reflections.getReflections(),
+                             equivalentIntensities, sigmaCritical, weightedZ);
+  }
   explicit PeaksStatistics(const UniqueReflectionCollection &reflections,
                            std::string &equivalentIntensities,
-                           double &sigmaCritical)
+                           double &sigmaCritical, bool &weightedZ)
       : m_measuredReflections(0), m_uniqueReflections(0), m_completeness(0.0),
         m_redundancy(0.0), m_rMerge(0.0), m_rPim(0.0), m_meanIOverSigma(0.0),
         m_dspacingMin(0.0), m_dspacingMax(0.0), m_chiSquared(0.0), m_peaks() {
     m_peaks.reserve(reflections.getObservedReflectionCount());
     calculatePeaksStatistics(reflections.getReflections(),
-                             equivalentIntensities, sigmaCritical);
+                             equivalentIntensities, sigmaCritical, weightedZ);
   }
 
   /// Total number of observed reflections - no symmetry is taken into
@@ -157,7 +168,7 @@ public:
 private:
   void calculatePeaksStatistics(
       const std::map<Kernel::V3D, UniqueReflection> &uniqueReflections,
-      std::string &equivalentIntensities, double &sigmaCritical);
+      std::string &equivalentIntensities, double &sigmaCritical, bool& weightedZ);
 
   double getIOverSigmaSum(const std::vector<double> &sigmas,
                           const std::vector<double> &intensities) const;

Framework/Crystal/inc/MantidCrystal/PredictPeaks.h

@@ -64,7 +64,7 @@ private:
 
   void calculateQAndAddToOutput(const Kernel::V3D &hkl,
                                 const Kernel::DblMatrix &orientedUB,
-                                const Kernel::DblMatrix &goniometerMatrix);
+                                const Kernel::DblMatrix &goniometerMatrix, size_t &seqNum);
 
 private:
   /// Get the predicted detector direction from Q

Framework/Crystal/src/LoadHKL.cpp

@@ -84,7 +84,8 @@ void LoadHKL::exec() {
     double SigI = std::stod(line.substr(20, 8));
     double wl = std::stod(line.substr(32, 8));
     double tbar, trans, scattering;
-    int run, bank, seqNum;
+    int run, bank;
+    size_t seqNum;
     if (cosines) {
       tbar = std::stod(line.substr(40, 8)); // tbar
       run = std::stoi(line.substr(102, 6));

Framework/Crystal/src/LoadIsawPeaks.cpp

@@ -256,7 +256,7 @@ std::string LoadIsawPeaks::readHeader(PeaksWorkspace_sptr outWS,
  */
 DataObjects::Peak LoadIsawPeaks::readPeak(PeaksWorkspace_sptr outWS,
                                           std::string &lastStr,
-                                          std::ifstream &in, int &seqNum,
+                                          std::ifstream &in, size_t &seqNum,
                                           std::string bankName, double qSign) {
   double h;
   double k;
@@ -268,8 +268,6 @@ DataObjects::Peak LoadIsawPeaks::readPeak(PeaksWorkspace_sptr outWS,
   double Inti;
   double SigI;
 
-  seqNum = -1;
-
   std::string s = lastStr;
 
   if (s.length() < 1 && in.good()) // blank line
@@ -489,7 +487,7 @@ void LoadIsawPeaks::appendFile(PeaksWorkspace_sptr outWS,
     oss << bankString << bankNum;
     std::string bankName = oss.str();
 
-    int seqNum = -1;
+    size_t seqNum;
 
     try {
       // Read the peak

Framework/Crystal/src/PeakStatisticsTools.cpp

@@ -57,7 +57,7 @@ std::vector<double> UniqueReflection::getSigmas() const {
 
 /// Removes peaks whose intensity deviates more than sigmaCritical from the
 /// intensities' mean.
-UniqueReflection UniqueReflection::removeOutliers(double sigmaCritical) const {
+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.");
@@ -67,8 +67,13 @@ UniqueReflection UniqueReflection::removeOutliers(double sigmaCritical) const {
 
   if (m_peaks.size() > 2) {
     auto intensities = getIntensities();
-    auto sigmas = getSigmas();
-    auto zScores = Kernel::getWeightedZscore(intensities, sigmas);
+    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) {
@@ -213,10 +218,11 @@ UniqueReflectionCollection::getReflections() const {
  * @param equivalentIntensities :: Mean or median for statistics of equivalent
  *peaks.
  * @param sigmaCritical :: Number of standard deviations for outliers.
+ * @param weightedZ :: True for weighted Zscore
  */
 void PeaksStatistics::calculatePeaksStatistics(
     const std::map<V3D, UniqueReflection> &uniqueReflections,
-    std::string &equivalentIntensities, double &sigmaCritical) {
+    std::string &equivalentIntensities, double &sigmaCritical, bool &weightedZ) {
   double rMergeNumerator = 0.0;
   double rPimNumerator = 0.0;
   double intensitySumRValues = 0.0;
@@ -230,7 +236,7 @@ void PeaksStatistics::calculatePeaksStatistics(
       ++m_uniqueReflections;
 
       // Possibly remove outliers.
-      auto outliersRemoved = unique.second.removeOutliers(sigmaCritical);
+      auto outliersRemoved = unique.second.removeOutliers(sigmaCritical, weightedZ);
 
       // I/sigma is calculated for all reflections, even if there is only one
       // observation.

Framework/Crystal/src/PredictPeaks.cpp

@@ -274,6 +274,7 @@ void PredictPeaks::exec() {
     HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
 
     size_t allowedPeakCount = 0;
+    size_t seqNum = 1;
 
     bool useExtendedDetectorSpace = getProperty("PredictPeaksOutsideDetectors");
     if (useExtendedDetectorSpace &&
@@ -284,8 +285,8 @@ void PredictPeaks::exec() {
 
     for (auto &possibleHKL : possibleHKLs) {
       if (lambdaFilter.isAllowed(possibleHKL)) {
+        calculateQAndAddToOutput(possibleHKL, orientedUB, goniometerMatrix, seqNum);
         ++allowedPeakCount;
-        calculateQAndAddToOutput(possibleHKL, orientedUB, goniometerMatrix);
       }
       prog.report();
     }
@@ -474,7 +475,7 @@ void PredictPeaks::setStructureFactorCalculatorFromSample(
  */
 void PredictPeaks::calculateQAndAddToOutput(const V3D &hkl,
                                             const DblMatrix &orientedUB,
-                                            const DblMatrix &goniometerMatrix) {
+                                            const DblMatrix &goniometerMatrix, size_t &seqNum) {
   // The q-vector direction of the peak is = goniometer * ub * hkl_vector
   // This is in inelastic convention: momentum transfer of the LATTICE!
   // Also, q does have a 2pi factor = it is equal to 2pi/wavelength.
@@ -535,6 +536,8 @@ void PredictPeaks::calculateQAndAddToOutput(const V3D &hkl,
   // Save the run number found before.
   peak->setRunNumber(m_runNumber);
   peak->setHKL(hkl * m_qConventionFactor);
+  peak->setPeakNumber(seqNum);
+  seqNum++;
 
   if (m_sfCalculator) {
     peak->setIntensity(m_sfCalculator->getFSquared(hkl));

Framework/Crystal/src/SortHKL.cpp

@@ -97,6 +97,9 @@ void SortHKL::init() {
       make_unique<WorkspaceProperty<MatrixWorkspace>>(
           "EquivalentsWorkspace", "EquivalentIntensities", Direction::Output),
       "Output Equivalent Intensities");
+  declareProperty("WeightedZScore", false,
+                  "Use weighted ZScore if true.\n"
+                  "If false, standard ZScore (default).");
 }
 
 void SortHKL::exec() {
@@ -116,6 +119,7 @@ void SortHKL::exec() {
       getUniqueReflections(peaks, cell);
   std::string equivalentIntensities = getPropertyValue("EquivalentIntensities");
   double sigmaCritical = getProperty("SigmaCritical");
+  bool weightedZ = getProperty("WeightedZScore");
 
   MatrixWorkspace_sptr UniqWksp =
       Mantid::API::WorkspaceFactory::Instance().create(
@@ -137,17 +141,24 @@ void SortHKL::exec() {
       counter++;
       auto wavelengths = unique.second.getWavelengths();
       auto intensities = unique.second.getIntensities();
-      auto sigmas = unique.second.getSigmas();
       g_log.debug() << "HKL " << unique.second.getHKL() << "\n";
       g_log.debug() << "Intensities ";
       for (const auto &e : intensities)
         g_log.debug() << e << "  ";
       g_log.debug() << "\n";
-      auto zScores = Kernel::getWeightedZscore(intensities, sigmas);
+      std::vector<double> zScores;
+      bool weighted = false;
+      if (!weighted) {
+        zScores = Kernel::getZscore(intensities);
+      } else {
+        auto sigmas = unique.second.getSigmas();
+        zScores = Kernel::getWeightedZscore(intensities, sigmas);
+      } 
+
       if (zScores.size() > maxPeaks)
         maxPeaks = zScores.size();
       // Possibly remove outliers.
-      auto outliersRemoved = unique.second.removeOutliers(sigmaCritical);
+      auto outliersRemoved = unique.second.removeOutliers(sigmaCritical, weightedZ);
 
       auto intensityStatistics =
           Kernel::getStatistics(outliersRemoved.getIntensities(),
@@ -211,7 +222,7 @@ void SortHKL::exec() {
   }
 
   PeaksStatistics peaksStatistics(uniqueReflections, equivalentIntensities,
-                                  sigmaCritical);
+                                  sigmaCritical, weightedZ);
 
   // Store the statistics for output.
   const std::string tableName = getProperty("StatisticsTable");

Framework/Crystal/src/StatisticsOfPeaksWorkspace.cpp

@@ -89,6 +89,9 @@ void StatisticsOfPeaksWorkspace::init() {
       make_unique<WorkspaceProperty<MatrixWorkspace>>(
           "EquivalentsWorkspace", "EquivalentIntensities", Direction::Output),
       "Output Equivalent Intensities");
+  declareProperty("WeightedZScore", false,
+                  "Use weighted ZScore if true.\n"
+                  "If false, standard ZScore (default).");
 }
 
 //----------------------------------------------------------------------------------------------
@@ -193,6 +196,7 @@ void StatisticsOfPeaksWorkspace::doSortHKL(Mantid::API::Workspace_sptr ws,
   std::string tableName = getPropertyValue("StatisticsTable");
   std::string equivalentIntensities = getPropertyValue("EquivalentIntensities");
   double sigmaCritical = getProperty("SigmaCritical");
+  bool weightedZ = getProperty("WeightedZScore");
   API::IAlgorithm_sptr statsAlg = createChildAlgorithm("SortHKL");
   statsAlg->setProperty("InputWorkspace", ws);
   statsAlg->setPropertyValue("OutputWorkspace", wkspName);
@@ -204,6 +208,7 @@ void StatisticsOfPeaksWorkspace::doSortHKL(Mantid::API::Workspace_sptr ws,
     statsAlg->setProperty("Append", true);
   statsAlg->setPropertyValue("EquivalentIntensities", equivalentIntensities);
   statsAlg->setProperty("SigmaCritical", sigmaCritical);
+  statsAlg->setProperty("WeightedZScore", weightedZ);
   statsAlg->executeAsChildAlg();
   PeaksWorkspace_sptr statsWksp = statsAlg->getProperty("OutputWorkspace");
   ITableWorkspace_sptr tablews = statsAlg->getProperty("StatisticsTable");

Framework/Crystal/test/PeakStatisticsToolsTest.h

@@ -161,14 +161,39 @@ public:
     auto cleanReflection = reflection.removeOutliers();
     TSM_ASSERT_EQUALS(
         "UniqueReflection removed outlier although it should not.",
-        cleanReflection.count(), 3);
+        cleanReflection.count(), 4);
 
     cleanReflection = reflection.removeOutliers(2.0);
     TSM_ASSERT_EQUALS(
         "UniqueReflection removed outlier although it should not.",
-        cleanReflection.count(), 2);
+        cleanReflection.count(), 4);
 
     cleanReflection = reflection.removeOutliers(1.0);
+    TSM_ASSERT_EQUALS(
+        "UniqueReflection did not remove outliers although it should have.",
+        cleanReflection.count(), 2);
+
+    std::vector<double> cleanIntensities = cleanReflection.getIntensities();
+    TS_ASSERT_EQUALS(cleanIntensities[0], 32.0);
+    TS_ASSERT_EQUALS(cleanIntensities[1], 31.0);
+  }
+
+  void test_UniqueReflectionRemoveOutliersWeighted() {
+    UniqueReflection reflection =
+        getReflectionWithPeaks({30.0, 34.0, 32.0, 31.0}, {4.5, 6.5, 10.0, 2.3});
+
+    // standard deviation is 1.70782512765993
+    auto cleanReflection = reflection.removeOutliers(3.0, true);
+    TSM_ASSERT_EQUALS(
+        "UniqueReflection removed outlier although it should not.",
+        cleanReflection.count(), 3);
+
+    cleanReflection = reflection.removeOutliers(2.0, true);
+    TSM_ASSERT_EQUALS(
+        "UniqueReflection removed outlier although it should not.",
+        cleanReflection.count(), 2);
+
+    cleanReflection = reflection.removeOutliers(1.0, true);
     TSM_ASSERT_EQUALS(
         "UniqueReflection did not remove outliers although it should have.",
         cleanReflection.count(), 1);
@@ -262,9 +287,7 @@ public:
         std::make_pair(V3D(1, 1, 1), UniqueReflection(V3D(1, 1, 1))));
     MockUniqueReflectionCollection reflections(uniques);
 
-    std::string statType = "Mean";
-    double sigmas = 3.0;
-    PeaksStatistics statistics(reflections, statType, sigmas);
+    PeaksStatistics statistics(reflections);
     TS_ASSERT_EQUALS(statistics.m_peaks.size(), 0);
     TS_ASSERT_EQUALS(statistics.m_uniqueReflections, 0);
     TS_ASSERT_EQUALS(statistics.m_redundancy, 0.0);
@@ -279,9 +302,7 @@ public:
         {{1, 1, 1}, getReflectionWithPeaks({56.0}, {4.5}, 1.0)}};
     MockUniqueReflectionCollection reflections(uniques);
 
-    std::string statType = "Mean";
-    double sigmas = 3.0;
-    PeaksStatistics statistics(reflections, statType, sigmas);
+    PeaksStatistics statistics(reflections);
     TS_ASSERT_EQUALS(statistics.m_peaks.size(), 1);
     TS_ASSERT_EQUALS(statistics.m_uniqueReflections, 1);
     TS_ASSERT_EQUALS(statistics.m_redundancy, 1.0);
@@ -297,9 +318,7 @@ public:
         {{1, 1, 2}, UniqueReflection(V3D(1, 1, 2))}};
     MockUniqueReflectionCollection reflections(uniques);
 
-    std::string statType = "Mean";
-    double sigmas = 3.0;
-    PeaksStatistics statistics(reflections, statType, sigmas);
+    PeaksStatistics statistics(reflections);
     TS_ASSERT_EQUALS(statistics.m_peaks.size(), 1);
     TS_ASSERT_EQUALS(statistics.m_uniqueReflections, 1);
     TS_ASSERT_EQUALS(statistics.m_redundancy, 1.0);
@@ -316,9 +335,7 @@ public:
         {{1, 1, 2}, getReflectionWithPeaks({20.0}, {1.0}, 2.0)}};
     MockUniqueReflectionCollection reflections(uniques);
 
-    std::string statType = "Mean";
-    double sigmas = 3.0;
-    PeaksStatistics statistics(reflections, statType, sigmas);
+    PeaksStatistics statistics(reflections);
     TS_ASSERT_EQUALS(statistics.m_peaks.size(), 2);
     TS_ASSERT_EQUALS(statistics.m_uniqueReflections, 2);
     TS_ASSERT_EQUALS(statistics.m_redundancy, 1.0);
@@ -333,9 +350,7 @@ public:
         {{1, 1, 1}, getReflectionWithPeaks({10.0, 20.0}, {0.1, 0.1}, 1.0)}};
     MockUniqueReflectionCollection reflections(uniques);
 
-    std::string statType = "Mean";
-    double sigmas = 3.0;
-    PeaksStatistics statistics(reflections, statType, sigmas);
+    PeaksStatistics statistics(reflections);
     TS_ASSERT_EQUALS(statistics.m_peaks.size(), 2);
     TS_ASSERT_EQUALS(statistics.m_uniqueReflections, 1);
     TS_ASSERT_EQUALS(statistics.m_redundancy, 2.0);
@@ -353,9 +368,7 @@ public:
          getReflectionWithPeaks({10.0, 20.0, 15.0}, {0.1, 0.1, 0.1}, 1.0)}};
     MockUniqueReflectionCollection reflections(uniques);
 
-    std::string statType = "Mean";
-    double sigmas = 3.0;
-    PeaksStatistics statistics(reflections, statType, sigmas);
+    PeaksStatistics statistics(reflections);
     TS_ASSERT_EQUALS(statistics.m_peaks.size(), 3);
     TS_ASSERT_EQUALS(statistics.m_uniqueReflections, 1);
     TS_ASSERT_EQUALS(statistics.m_redundancy, 3.0);

Framework/Crystal/test/SortHKLTest.h

@@ -89,7 +89,7 @@ public:
     TS_ASSERT(wsout);
     if (!wsout)
       return;
-    TS_ASSERT_EQUALS(wsout->getNumberPeaks(), 16);
+    TS_ASSERT_EQUALS(wsout->getNumberPeaks(), 24);
 
     Peak p = wsout->getPeaks()[0];
     TS_ASSERT_EQUALS(p.getH(), 1);

Framework/Crystal/test/StatisticsOfPeaksWorkspaceTest.h

@@ -81,9 +81,9 @@ public:
     TS_ASSERT_EQUALS(tableOut->String(0, 0), "Overall");
     TS_ASSERT_EQUALS(tableOut->Int(0, 1), 3);
     TS_ASSERT_DELTA(tableOut->Double(0, 2), 3.6, .1);
-    TS_ASSERT_DELTA(tableOut->Double(0, 3), 14.03, .1);
-    TS_ASSERT_DELTA(tableOut->Double(0, 4), 5.3, .1);
-    TS_ASSERT_DELTA(tableOut->Double(0, 5), 1.25, .1);
+    TS_ASSERT_DELTA(tableOut->Double(0, 3), 14.3, .1);
+    TS_ASSERT_DELTA(tableOut->Double(0, 4), 8., .1);
+    TS_ASSERT_DELTA(tableOut->Double(0, 5), 1.4195, .1);
     TS_ASSERT_DELTA(tableOut->Double(0, 6), 0., .1);
     TS_ASSERT_DELTA(tableOut->Double(0, 7), 0., .1);
 
@@ -94,7 +94,7 @@ public:
     TS_ASSERT(wsout);
     if (!wsout)
       return;
-    TS_ASSERT_EQUALS(wsout->getNumberPeaks(), 16);
+    TS_ASSERT_EQUALS(wsout->getNumberPeaks(), 24);
 
     Peak p = wsout->getPeaks()[0];
     TS_ASSERT_EQUALS(p.getH(), 1);

Framework/DataObjects/inc/MantidDataObjects/Peak.h

@@ -150,7 +150,7 @@ public:
   int getCol() const override;
   void setRow(int m_row);
   void setCol(int m_col);
-  void setPeakNumber(int m_peakNumber) override;
+  void setPeakNumber(size_t m_peakNumber) override;
   int getPeakNumber() const override;
 
   virtual Mantid::Kernel::V3D getDetPos() const override;

Framework/DataObjects/src/Peak.cpp

@@ -900,7 +900,7 @@ void Peak::setCol(int m_col) { this->m_col = m_col; }
 // -------------------------------------------------------------------------------------
 /** Sets the unique peak number
  * @param m_peakNumber :: unique peak number value   */
-void Peak::setPeakNumber(int m_peakNumber) {
+void Peak::setPeakNumber(size_t m_peakNumber) {
   this->m_peakNumber = m_peakNumber;
 }
 

Framework/Geometry/inc/MantidGeometry/Crystal/IPeak.h

@@ -78,7 +78,7 @@ public:
   virtual void setBinCount(double m_BinCount) = 0;
 
   virtual int getPeakNumber() const = 0;
-  virtual void setPeakNumber(int m_PeakNumber) = 0;
+  virtual void setPeakNumber(size_t m_PeakNumber) = 0;
 
   virtual Mantid::Kernel::Matrix<double> getGoniometerMatrix() const = 0;
   virtual void setGoniometerMatrix(

Framework/Geometry/test/MockObjects.h

@@ -67,7 +67,7 @@ public:
   MOCK_CONST_METHOD0(getRunNumber, int());
   MOCK_CONST_METHOD0(getPeakNumber, int());
   MOCK_METHOD1(setRunNumber, void(int m_RunNumber));
-  MOCK_METHOD1(setPeakNumber, void(int m_PeakNumber));
+  MOCK_METHOD1(setPeakNumber, void(size_t m_PeakNumber));
   MOCK_CONST_METHOD0(getMonitorCount, double());
   MOCK_METHOD1(setMonitorCount, void(double m_MonitorCount));
   MOCK_CONST_METHOD0(getH, double());

Framework/MDAlgorithms/src/FindPeaksMD.cpp

@@ -748,7 +748,7 @@ void FindPeaksMD::exec() {
 
   for (auto i = 0; i != peakWS->getNumberPeaks(); ++i) {
     Peak &p = peakWS->getPeak(i);
-    p.setPeakNumber(static_cast<int>(i + 1));
+    p.setPeakNumber(i + 1);
   }
   // Save the output
   setProperty("OutputWorkspace", peakWS);

Testing/SystemTests/tests/analysis/SortHKLTest.py

@@ -116,7 +116,7 @@ class SortHKLTest(HKLStatisticsTestMixin, stresstesting.MantidStressTest):
         centering_name = self._centering_map[space_group[0]]
 
         # pylint: disable=unused-variable
-        sorted_hkls, chi2, statistics = SortHKL(InputWorkspace=reflections,
+        sorted_hkls, chi2, statistics, equivInten = SortHKL(InputWorkspace=reflections,
                                                 PointGroup=point_group_name,
                                                 LatticeCentering=centering_name)