Re #5571. First cut at algorithm to subtract matching peaks.

Processes the peaks in the rhs workspace, removing the first match
it finds (in Q, within the given tolerance) from the peaks in the
lhs workspace.

Code/Mantid/Framework/Crystal/CMakeLists.txt

@@ -4,6 +4,7 @@ set ( SRC_FILES
 	src/CalculateUMatrix.cpp
 	src/CentroidPeaks.cpp
 	src/CombinePeaksWorkspaces.cpp
+	src/DiffPeaksWorkspaces.cpp
 	src/FindSXPeaks.cpp
 	src/FindUBUsingFFT.cpp
 	src/FindUBUsingIndexedPeaks.cpp
@@ -43,6 +44,7 @@ set ( INC_FILES
 	inc/MantidCrystal/CalculateUMatrix.h
 	inc/MantidCrystal/CentroidPeaks.h
 	inc/MantidCrystal/CombinePeaksWorkspaces.h
+	inc/MantidCrystal/DiffPeaksWorkspaces.h
 	inc/MantidCrystal/FindSXPeaks.h
 	inc/MantidCrystal/FindUBUsingFFT.h
 	inc/MantidCrystal/FindUBUsingIndexedPeaks.h
@@ -82,6 +84,7 @@ set ( TEST_FILES
 	test/CalculateUMatrixTest.h
 	test/CentroidPeaksTest.h
 	test/CombinePeaksWorkspacesTest.h
+	test/DiffPeaksWorkspacesTest.h
 	test/FindSXPeaksTest.h
 	test/FindUBUsingFFTTest.h
 	test/FindUBUsingIndexedPeaksTest.h

Code/Mantid/Framework/Crystal/inc/MantidCrystal/DiffPeaksWorkspaces.h

@@ -0,0 +1,51 @@
+#ifndef MANTID_CRYSTAL_DIFFPEAKSWORKSPACES_H_
+#define MANTID_CRYSTAL_DIFFPEAKSWORKSPACES_H_
+
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid
+{
+namespace Crystal
+{
+  /** An algorithm that subtracts from a workspace (the LHSWorkspace) any peaks that match
+      entries in a second workspace (the RHSWorkspace). Such peaks are those that are within
+      the given tolerance in all components of Q. Note that a peak in the RHSWorkspace will
+      only be matched to the first in the LHSWorkspace that is within tolerance.
+
+      Copyright © 2012 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+      This file is part of Mantid.
+
+      Mantid is free software; you can redistribute it and/or modify
+      it under the terms of the GNU General Public License as published by
+      the Free Software Foundation; either version 3 of the License, or
+      (at your option) any later version.
+
+      Mantid is distributed in the hope that it will be useful,
+      but WITHOUT ANY WARRANTY; without even the implied warranty of
+      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+      GNU General Public License for more details.
+
+      You should have received a copy of the GNU General Public License
+      along with this program.  If not, see <http://www.gnu.org/licenses/>.
+  */
+  class DLLExport DiffPeaksWorkspaces  : public API::Algorithm
+  {
+  public:
+    DiffPeaksWorkspaces();
+    virtual ~DiffPeaksWorkspaces();
+
+    virtual const std::string name() const;
+    virtual int version() const;
+    virtual const std::string category() const;
+
+  private:
+    virtual void initDocs();
+    void init();
+    void exec();
+  };
+
+} // namespace Crystal
+} // namespace Mantid
+
+#endif  /* MANTID_CRYSTAL_DIFFPEAKSWORKSPACES_H_ */

Code/Mantid/Framework/Crystal/src/DiffPeaksWorkspaces.cpp

@@ -0,0 +1,118 @@
+/*WIKI*
+TODO: Enter a full wiki-markup description of your algorithm here. You can then use the Build/wiki_maker.py script to generate your full wiki page.
+*WIKI*/
+
+#include "MantidCrystal/DiffPeaksWorkspaces.h"
+#include "MantidKernel/BoundedValidator.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+
+namespace Mantid
+{
+namespace Crystal
+{
+  // Register the algorithm into the AlgorithmFactory
+  DECLARE_ALGORITHM(DiffPeaksWorkspaces)
+
+  using namespace Kernel;
+  using namespace API;
+  using DataObjects::PeaksWorkspace;
+  using DataObjects::PeaksWorkspace_const_sptr;
+  using DataObjects::PeaksWorkspace_sptr;
+  using DataObjects::Peak;
+
+  /** Constructor
+   */
+  DiffPeaksWorkspaces::DiffPeaksWorkspaces()
+  {
+  }
+
+  /** Destructor
+   */
+  DiffPeaksWorkspaces::~DiffPeaksWorkspaces()
+  {
+  }
+
+  /// Algorithm's name for identification. @see Algorithm::name
+  const std::string DiffPeaksWorkspaces::name() const { return "DiffPeaksWorkspaces"; };
+  /// Algorithm's version for identification. @see Algorithm::version
+  int DiffPeaksWorkspaces::version() const { return 1; };
+  /// Algorithm's category for identification. @see Algorithm::category
+  const std::string DiffPeaksWorkspaces::category() const { return "Crystal"; }
+
+  /// Sets documentation strings for this algorithm
+  void DiffPeaksWorkspaces::initDocs()
+  {
+    this->setWikiSummary("Removes from a workspace any peaks that match a peak in a second workspace.");
+    this->setOptionalMessage("Removes from a workspace any peaks that match a peak in a second workspace.");
+  }
+
+  /** Initialises the algorithm's properties.
+   */
+  void DiffPeaksWorkspaces::init()
+  {
+    declareProperty(new WorkspaceProperty<PeaksWorkspace>("LHSWorkspace","",Direction::Input),
+        "The first of peaks.");
+    declareProperty(new WorkspaceProperty<PeaksWorkspace>("RHSWorkspace","",Direction::Input),
+        "The second set of peaks.");
+    declareProperty(new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace","",Direction::Output),
+        "The set of peaks that are in the first, but not the second, workspace.");
+
+    auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
+    mustBePositive->setLower(0.0);
+    // N.B. Andrei reckons it should be delta_q/q
+    declareProperty("Tolerance", 0.0, mustBePositive,
+        "Maximum difference in each component of Q for which peaks are considered identical");
+  }
+
+  /** Executes the algorithm.
+   */
+  void DiffPeaksWorkspaces::exec()
+  {
+    PeaksWorkspace_const_sptr LHSWorkspace = getProperty("LHSWorkspace");
+    PeaksWorkspace_const_sptr RHSWorkspace = getProperty("RHSWorkspace");
+    const double Tolerance = getProperty("Tolerance");
+
+    // Warn if not the same instrument, sample
+    if ( LHSWorkspace->getInstrument()->getName() != RHSWorkspace->getInstrument()->getName() )
+    {
+      g_log.warning("The two input workspaces do not appear to come from data take on the same instrument");
+    }
+    if ( LHSWorkspace->sample().getName() != RHSWorkspace->sample().getName() )
+    {
+      g_log.warning("The two input workspaces do not appear to relate to the same sample");
+    }
+
+    // Copy the first workspace to our output workspace
+    PeaksWorkspace_sptr output(LHSWorkspace->clone());
+    // Get hold of the peaks in the second workspace
+    auto & rhsPeaks = RHSWorkspace->getPeaks();
+    // Get hold of the peaks in the first workspace as we'll need to examine them
+    auto & lhsPeaks = output->getPeaks();
+
+    Progress progress(this, 0, 1, rhsPeaks.size());
+
+    // Loop over the peaks in the second workspace, searching for a match in the first
+    for ( size_t i = 0; i < rhsPeaks.size(); ++i )
+    {
+      const Peak& currentPeak = rhsPeaks[i];
+      // Now have to go through the first workspace checking for matches
+      // Not doing anything clever as peaks workspace are typically not large - just a linear search
+      for ( int j = 0; j < output->getNumberPeaks(); ++j )
+      {
+        const V3D deltaQ = currentPeak.getQSampleFrame() - lhsPeaks[j].getQSampleFrame();
+        if ( deltaQ.nullVector(Tolerance) )  // Using a V3D method that does the job
+        {
+          // As soon as we find a match, remove it from the output and move onto the next rhs peak
+          output->removePeak(j);
+          break;
+        }
+      }
+
+      progress.report();
+    }
+
+    setProperty("OutputWorkspace", output);
+  }
+
+} // namespace Crystal
+} // namespace Mantid

Code/Mantid/Framework/Crystal/test/DiffPeaksWorkspacesTest.h

@@ -0,0 +1,195 @@
+#ifndef MANTID_CRYSTAL_DIFFPEAKSWORKSPACESTEST_H_
+#define MANTID_CRYSTAL_DIFFPEAKSWORKSPACESTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidCrystal/DiffPeaksWorkspaces.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+
+using Mantid::Crystal::DiffPeaksWorkspaces;
+
+class DiffPeaksWorkspacesTest : public CxxTest::TestSuite
+{
+public:
+  // This pair of boilerplate methods prevent the suite being created statically
+  // This means the constructor isn't called when running other tests
+  static DiffPeaksWorkspacesTest *createSuite() { return new DiffPeaksWorkspacesTest(); }
+  static void destroySuite( DiffPeaksWorkspacesTest *suite ) { delete suite; }
+
+
+  void test_init()
+  {
+    DiffPeaksWorkspaces alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT( alg.isInitialized() )
+  }
+
+  void test_invalid_input()
+  {
+    DiffPeaksWorkspaces alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    // Tolerance has to be positive.
+    TS_ASSERT_THROWS( alg.setProperty("Tolerance", -1.0), std::invalid_argument )
+  }
+
+  // It shouldn't be a problem to subtract identical workspaces. You just get an empty one out.
+  void test_diff_identical_peaksworkspaces()
+  {
+    using namespace Mantid::API;
+    using namespace Mantid::DataObjects;
+
+    PeaksWorkspace_sptr inWS = WorkspaceCreationHelper::createPeaksWorkspace();
+
+    // Name of the output workspace.
+    std::string outWSName("CombinePeaksWorkspacesTest_OutputWS");
+
+    DiffPeaksWorkspaces alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("LHSWorkspace", inWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("RHSWorkspace", inWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) )
+    TS_ASSERT( alg.execute() )
+
+    // Retrieve the workspace from data service.
+    IPeaksWorkspace_const_sptr ws;
+    TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<IPeaksWorkspace>(outWSName) );
+    TS_ASSERT(ws);
+    if (!ws) return;
+
+    TS_ASSERT_EQUALS( ws->getNumberPeaks(), 0 )
+
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(outWSName);
+  }
+
+  // In fact, it should be fine for the second one to be larger
+  void test_diff_larger_rhs()
+  {
+    using namespace Mantid::API;
+    using namespace Mantid::DataObjects;
+
+    PeaksWorkspace_sptr lhsWS = WorkspaceCreationHelper::createPeaksWorkspace();
+    PeaksWorkspace_sptr rhsWS = WorkspaceCreationHelper::createPeaksWorkspace(6);
+
+    // Name of the output workspace.
+    std::string outWSName("CombinePeaksWorkspacesTest_OutputWS");
+
+    DiffPeaksWorkspaces alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("LHSWorkspace", lhsWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("RHSWorkspace", rhsWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) )
+    TS_ASSERT( alg.execute() )
+
+    // Retrieve the workspace from data service.
+    IPeaksWorkspace_const_sptr ws;
+    TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<IPeaksWorkspace>(outWSName) );
+    TS_ASSERT(ws);
+    if (!ws) return;
+
+    TS_ASSERT_EQUALS( ws->getNumberPeaks(), 0 )
+
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(outWSName);
+  }
+
+  void test_no_matches()
+  {
+    using namespace Mantid::API;
+    using namespace Mantid::DataObjects;
+
+    PeaksWorkspace_sptr lhsWS = WorkspaceCreationHelper::createPeaksWorkspace();
+    PeaksWorkspace_sptr rhsWS = WorkspaceCreationHelper::createPeaksWorkspace();
+
+    // Shift the rhs peaks in Q
+    auto & rhsPeaks = rhsWS->getPeaks();
+    rhsPeaks[0].setWavelength(rhsPeaks[0].getWavelength()*1.01);
+    rhsPeaks[1].setDetectorID(50);
+
+    // Name of the output workspace.
+    std::string outWSName("CombinePeaksWorkspacesTest_OutputWS");
+
+    DiffPeaksWorkspaces alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("LHSWorkspace", lhsWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("RHSWorkspace", rhsWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("Tolerance", 0.1) )
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) )
+    TS_ASSERT( alg.execute() )
+
+    // Retrieve the workspace from data service.
+    IPeaksWorkspace_const_sptr ws;
+    TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<IPeaksWorkspace>(outWSName) );
+    TS_ASSERT(ws);
+    if (!ws) return;
+
+    TS_ASSERT_EQUALS( ws->getNumberPeaks(), 2 )
+    TS_ASSERT_EQUALS( ws->getPeak(0).getQLabFrame(), lhsWS->getPeak(0).getQLabFrame() )
+    TS_ASSERT_EQUALS( ws->getPeak(1).getQLabFrame(), lhsWS->getPeak(1).getQLabFrame() )
+
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(outWSName);
+  }
+
+  void test_match_peaks_within_tolerance()
+  {
+    using namespace Mantid::API;
+    using namespace Mantid::DataObjects;
+
+    PeaksWorkspace_sptr lhsWS = WorkspaceCreationHelper::createPeaksWorkspace(6);
+    PeaksWorkspace_sptr rhsWS = WorkspaceCreationHelper::createPeaksWorkspace(4);
+
+    // Slightly adjust the peaks in one of the workspaces
+    auto & rhsPeaks = rhsWS->getPeaks();
+    auto & lhsPeaks = lhsWS->getPeaks();
+
+    // Need to change a couple of detector IDs so that I can get peaks with larger |Q_z| than |Q_x|
+    lhsPeaks[2].setDetectorID(50);
+    lhsPeaks[3].setDetectorID(51);
+    rhsPeaks[2].setDetectorID(50);
+    rhsPeaks[3].setDetectorID(51);
+
+    // And need to shift some peaks in one workspace to test the delta checking
+    // This one will fails to match in x & z
+    rhsPeaks[0].setWavelength(rhsPeaks[0].getWavelength()*1.01);
+    // This one matches in z but not in x
+    rhsPeaks[1].setWavelength(rhsPeaks[1].getWavelength()*1.02);
+    // This one matches in x but not z
+    rhsPeaks[2].setWavelength(rhsPeaks[2].getWavelength()*1.0335);
+    // This one will be matched (to lhsPeaks[0]) and will not appear in the output
+    rhsPeaks[3].setWavelength(rhsPeaks[3].getWavelength()*1.04);
+
+    // Name of the output workspace.
+    std::string outWSName("CombinePeaksWorkspacesTest_OutputWS");
+
+    DiffPeaksWorkspaces alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("LHSWorkspace", lhsWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("RHSWorkspace", rhsWS) )
+    TS_ASSERT_THROWS_NOTHING( alg.setProperty("Tolerance", 0.08145) )
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) )
+    TS_ASSERT( alg.execute() )
+
+    // Retrieve the workspace from data service.
+    IPeaksWorkspace_const_sptr ws;
+    TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<IPeaksWorkspace>(outWSName) );
+    TS_ASSERT(ws);
+    if (!ws) return;
+
+    TS_ASSERT_EQUALS( ws->getNumberPeaks(), 5 )
+    TS_ASSERT_EQUALS( ws->getPeak(0).getQLabFrame(), lhsWS->getPeak(0).getQLabFrame() )
+    TS_ASSERT_EQUALS( ws->getPeak(1).getQLabFrame(), lhsWS->getPeak(1).getQLabFrame() )
+    TS_ASSERT_EQUALS( ws->getPeak(2).getQLabFrame(), lhsWS->getPeak(2).getQLabFrame() )
+    TS_ASSERT_EQUALS( ws->getPeak(3).getQLabFrame(), lhsWS->getPeak(4).getQLabFrame() )
+    TS_ASSERT_EQUALS( ws->getPeak(4).getQLabFrame(), lhsWS->getPeak(5).getQLabFrame() )
+    TS_ASSERT_EQUALS( ws->getInstrument()->baseInstrument(), lhsWS->getInstrument()->baseInstrument() )
+
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(outWSName);
+  }
+
+};
+
+
+#endif /* MANTID_CRYSTAL_DIFFPEAKSWORKSPACESTEST_H_ */