MatrixWorkspaceTest.h

#ifndef WORKSPACETEST_H_
#define WORKSPACETEST_H_

#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/VMD.h"
#include "MantidKernel/make_cow.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include "MantidTestHelpers/FakeObjects.h"
#include "MantidTestHelpers/InstrumentCreationHelper.h"
#include "MantidTestHelpers/NexusTestHelper.h"
#include "MantidTestHelpers/ParallelRunner.h"
#include "MantidTypes/SpectrumDefinition.h"
#include "PropertyManagerHelper.h"

#include <cxxtest/TestSuite.h>

#include <boost/make_shared.hpp>

#include <algorithm>
#include <atomic>
#include <cmath>
#include <functional>
#include <numeric>

using std::size_t;
using namespace Mantid;
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using Mantid::Indexing::IndexInfo;
using Mantid::Types::Core::DateAndTime;

// Declare into the factory.
DECLARE_WORKSPACE(WorkspaceTester)

namespace {
/** Create a workspace with numSpectra, with
 * each spectrum having one detector, at id = workspace index.
 * @param numSpectra
 * @return
 */
boost::shared_ptr<MatrixWorkspace> makeWorkspaceWithDetectors(size_t numSpectra,
                                                              size_t numBins) {
  boost::shared_ptr<MatrixWorkspace> ws2 =
      boost::make_shared<WorkspaceTester>();
  ws2->initialize(numSpectra, numBins, numBins);

  auto inst = boost::make_shared<Instrument>("TestInstrument");
  // We get a 1:1 map by default so the detector ID should match the spectrum
  // number
  for (size_t i = 0; i < ws2->getNumberHistograms(); ++i) {
    // Create a detector for each spectra
    Detector *det = new Detector("pixel", static_cast<detid_t>(i), inst.get());
    det->setShape(
        ComponentCreationHelper::createSphere(0.01, V3D(0, 0, 0), "1"));
    inst->add(det);
    inst->markAsDetector(det);
    ws2->getSpectrum(i).addDetectorID(static_cast<detid_t>(i));
  }
  ws2->setInstrument(inst);
  return ws2;
}

void run_legacy_setting_spectrum_numbers_with_MPI(
    const Parallel::Communicator &comm) {
  using namespace Parallel;
  for (const auto storageMode : {StorageMode::MasterOnly, StorageMode::Cloned,
                                 StorageMode::Distributed}) {
    WorkspaceTester ws;
    if (comm.rank() == 0 || storageMode != StorageMode::MasterOnly) {
      Indexing::IndexInfo indexInfo(1000, storageMode, comm);
      ws.initialize(indexInfo,
                    HistogramData::Histogram(HistogramData::Points(1)));
    }
    if (storageMode == StorageMode::Distributed && comm.size() > 1) {
      TS_ASSERT_THROWS_EQUALS(ws.getSpectrum(0).setSpectrumNo(42),
                              const std::logic_error &e, std::string(e.what()),
                              "Setting spectrum numbers in MatrixWorkspace via "
                              "ISpectrum::setSpectrumNo is not possible in MPI "
                              "runs for distributed workspaces. Use "
                              "IndexInfo.");
    } else {
      if (comm.rank() == 0 || storageMode != StorageMode::MasterOnly) {
        TS_ASSERT_THROWS_NOTHING(ws.getSpectrum(0).setSpectrumNo(42));
      }
    }
  }
}
} // namespace

class MatrixWorkspaceTest : 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 MatrixWorkspaceTest *createSuite() {
    return new MatrixWorkspaceTest();
  }
  static void destroySuite(MatrixWorkspaceTest *suite) { delete suite; }

  MatrixWorkspaceTest() : ws(boost::make_shared<WorkspaceTester>()) {
    ws->initialize(1, 1, 1);
  }

  void test_indexInfo() {
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    const auto &indexInfo = ws.indexInfo();
    // IndexInfo contains spectrum numbers generated by WorkspaceTester.
    TS_ASSERT_EQUALS(indexInfo.spectrumNumber(0), 1);
    TS_ASSERT_EQUALS(indexInfo.spectrumNumber(1), 2);
    TS_ASSERT_EQUALS(indexInfo.spectrumNumber(2), 3);
    // Workspace tester contains detector IDs...
    TS_ASSERT_EQUALS(ws.getSpectrum(0).getDetectorIDs().size(), 1);
    TS_ASSERT_EQUALS(ws.getSpectrum(1).getDetectorIDs().size(), 1);
    TS_ASSERT_EQUALS(ws.getSpectrum(2).getDetectorIDs().size(), 1);
    // ... but spectrum definitions in IndexInfo are empty...
    TS_ASSERT_EQUALS((*indexInfo.spectrumDefinitions())[0],
                     SpectrumDefinition{});
    TS_ASSERT_EQUALS((*indexInfo.spectrumDefinitions())[1],
                     SpectrumDefinition{});
    TS_ASSERT_EQUALS((*indexInfo.spectrumDefinitions())[2],
                     SpectrumDefinition{});
    // ... since there is no instrument, i.e., all detector IDs are invalid.
    TS_ASSERT_EQUALS(ws.detectorInfo().size(), 0);
  }

  void test_setIndexInfo_size_failure() {
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    IndexInfo bad(2);
    TS_ASSERT_THROWS_EQUALS(ws.setIndexInfo(std::move(bad)),
                            const std::invalid_argument &e,
                            std::string(e.what()),
                            "MatrixWorkspace::setIndexInfo: IndexInfo size "
                            "does not match number of histograms in workspace");
  }

  void test_setIndexInfo_bad_detector_index() {
    WorkspaceTester ws;
    ws.initialize(1, 1, 1);
    IndexInfo indices(1);
    indices.setSpectrumNumbers({2});
    std::vector<SpectrumDefinition> specDefs(1);
    specDefs[0].add(0);
    indices.setSpectrumDefinitions(specDefs);
    TS_ASSERT_THROWS_EQUALS(ws.setIndexInfo(std::move(indices)),
                            const std::invalid_argument &e,
                            std::string(e.what()),
                            "MatrixWorkspace: SpectrumDefinition contains an "
                            "out-of-range detector index, i.e., the spectrum "
                            "definition does not match the instrument in the "
                            "workspace.");
  }

  void test_setIndexInfo_bad_detector_time_index() {
    WorkspaceTester ws;
    ws.initialize(1, 1, 1);
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 1));
    IndexInfo indices(1);
    indices.setSpectrumNumbers({2});
    std::vector<SpectrumDefinition> specDefs(1);
    specDefs[0].add(0, 1);
    indices.setSpectrumDefinitions(specDefs);
    TS_ASSERT_THROWS_EQUALS(ws.setIndexInfo(std::move(indices)),
                            const std::invalid_argument &e,
                            std::string(e.what()),
                            "MatrixWorkspace: SpectrumDefinition contains an "
                            "out-of-range time index for a detector, i.e., the "
                            "spectrum definition does not match the instrument "
                            "in the workspace.");
  }

  void test_setIndexInfo_default_mapping_failure() {
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    // 2x2 = 4 pixels
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    IndexInfo indices(3);
    TS_ASSERT_THROWS_EQUALS(ws.setIndexInfo(std::move(indices)),
                            const std::invalid_argument &e,
                            std::string(e.what()),
                            "MatrixWorkspace: IndexInfo does not contain "
                            "spectrum definitions so building a 1:1 mapping "
                            "from spectra to detectors was attempted, but the "
                            "number of spectra in the workspace is not equal "
                            "to the number of detectors in the instrument.");
  }

  void test_setIndexInfo_default_mapping() {
    WorkspaceTester ws;
    ws.initialize(4, 1, 1);
    // 2x2 = 4 pixels
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    IndexInfo indices(4);

    TS_ASSERT_THROWS_NOTHING(ws.setIndexInfo(std::move(indices)));

    const auto &info = ws.indexInfo();
    TS_ASSERT(info.spectrumDefinitions());
    std::vector<SpectrumDefinition> specDefs(4);
    specDefs[0].add(0);
    specDefs[1].add(1);
    specDefs[2].add(2);
    specDefs[3].add(3);
    TS_ASSERT_EQUALS(*info.spectrumDefinitions(), specDefs);
    TS_ASSERT_EQUALS(ws.getSpectrum(0).getSpectrumNo(), 1);
    TS_ASSERT_EQUALS(ws.getSpectrum(1).getSpectrumNo(), 2);
    TS_ASSERT_EQUALS(ws.getSpectrum(2).getSpectrumNo(), 3);
    TS_ASSERT_EQUALS(ws.getSpectrum(3).getSpectrumNo(), 4);
    TS_ASSERT_EQUALS(ws.getSpectrum(0).getDetectorIDs(),
                     (std::set<detid_t>{4}));
    TS_ASSERT_EQUALS(ws.getSpectrum(1).getDetectorIDs(),
                     (std::set<detid_t>{5}));
    TS_ASSERT_EQUALS(ws.getSpectrum(2).getDetectorIDs(),
                     (std::set<detid_t>{6}));
    TS_ASSERT_EQUALS(ws.getSpectrum(3).getDetectorIDs(),
                     (std::set<detid_t>{7}));
  }

  void test_setIndexInfo_updates_ISpectrum() {
    // NOTE: This test checks if the IndexInfo set via
    // MatrixWorkspace::setIndexInfo() affects data stored in ISpectrum and
    // obtained via the legacy interface. THIS TEST SHOULD BE REMOVED ONCE THAT
    // INTERFACE IS BEING REMOVED.
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    TS_ASSERT_EQUALS(ws.detectorInfo().size(), 4);
    IndexInfo indices(3);
    indices.setSpectrumNumbers({2, 4, 6});
    std::vector<SpectrumDefinition> specDefs(3);
    specDefs[0].add(0);
    specDefs[1].add(1);
    specDefs[2].add(2);
    specDefs[2].add(3);
    indices.setSpectrumDefinitions(specDefs);
    TS_ASSERT_THROWS_NOTHING(ws.setIndexInfo(std::move(indices)));
    TS_ASSERT_EQUALS(ws.getSpectrum(0).getSpectrumNo(), 2);
    TS_ASSERT_EQUALS(ws.getSpectrum(1).getSpectrumNo(), 4);
    TS_ASSERT_EQUALS(ws.getSpectrum(2).getSpectrumNo(), 6);
    TS_ASSERT_EQUALS(ws.getSpectrum(0).getDetectorIDs(),
                     (std::set<detid_t>{4}));
    TS_ASSERT_EQUALS(ws.getSpectrum(1).getDetectorIDs(),
                     (std::set<detid_t>{5}));
    TS_ASSERT_EQUALS(ws.getSpectrum(2).getDetectorIDs(),
                     (std::set<detid_t>{6, 7}));
  }

  void test_indexInfo_legacy_compatibility() {
    // NOTE: This test checks if the IndexInfo reference returned by
    // MatrixWorkspace::indexInfo() reflects changes done via the legacy
    // interface of ISpectrum. THIS TEST SHOULD BE REMOVED ONCE THAT INTERFACE
    // IS BEING REMOVED.
    WorkspaceTester ws;
    ws.initialize(1, 1, 1);
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    const auto &indexInfo = ws.indexInfo();
    TS_ASSERT_EQUALS(indexInfo.spectrumNumber(0), 1);
    TS_ASSERT_EQUALS((*indexInfo.spectrumDefinitions())[0],
                     SpectrumDefinition{});
    ws.getSpectrum(0).setSpectrumNo(7);
    ws.getSpectrum(0).setDetectorID(7);
    // No changes -- old and new interface should not be mixed!
    TS_ASSERT_EQUALS(indexInfo.spectrumNumber(0), 1);
    TS_ASSERT_EQUALS((*indexInfo.spectrumDefinitions())[0],
                     SpectrumDefinition{});
    // After getting a new reference we should see the changes.
    const auto &indexInfo2 = ws.indexInfo();
    TS_ASSERT_EQUALS(indexInfo2.spectrumNumber(0), 7);
    SpectrumDefinition specDef;
    specDef.add(ws.detectorInfo().indexOf(7));
    TS_ASSERT_EQUALS((*indexInfo.spectrumDefinitions())[0], specDef);
  }

  void test_IndexInfo_copy() {
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    IndexInfo indices(3);
    indices.setSpectrumNumbers({2, 4, 6});
    std::vector<SpectrumDefinition> specDefs(3);
    specDefs[0].add(0);
    specDefs[1].add(1);
    specDefs[2].add(2);
    specDefs[2].add(3);
    indices.setSpectrumDefinitions(specDefs);
    ws.setIndexInfo(std::move(indices));

    // Internally this references data in ISpectrum
    const auto &indexInfo = ws.indexInfo();
    // This should create a copy, dropping any links to MatrixWorkspace or
    // ISpectrum
    const auto copy(indexInfo);

    TS_ASSERT_EQUALS(copy.spectrumNumber(0), 2);
    TS_ASSERT_EQUALS(copy.spectrumNumber(1), 4);
    TS_ASSERT_EQUALS(copy.spectrumNumber(2), 6);
    TS_ASSERT_EQUALS(*copy.spectrumDefinitions(), specDefs);
    // Changing data in workspace affects indexInfo, but not copy
    ws.getSpectrum(0).setSpectrumNo(7);
    ws.getSpectrum(0).addDetectorID(7);
    const auto &indexInfo2 = ws.indexInfo();
    TS_ASSERT_EQUALS(indexInfo2.spectrumNumber(0), 7);
    SpectrumDefinition specDef;
    specDef.add(ws.detectorInfo().indexOf(4));
    specDef.add(ws.detectorInfo().indexOf(7));
    TS_ASSERT_EQUALS((*indexInfo2.spectrumDefinitions())[0], specDef);
    TS_ASSERT_EQUALS(copy.spectrumNumber(0), 2);
    TS_ASSERT_EQUALS((*copy.spectrumDefinitions())[0], specDefs[0]);
  }

  void test_setIndexInfo_shares_spectrumDefinition() {
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    IndexInfo indices(3);
    indices.setSpectrumNumbers({2, 4, 6});

    auto defs = Kernel::make_cow<std::vector<SpectrumDefinition>>(3);
    TS_ASSERT_THROWS_NOTHING(indices.setSpectrumDefinitions(defs));
    TS_ASSERT_THROWS_NOTHING(ws.setIndexInfo(std::move(indices)));
    TS_ASSERT_EQUALS(ws.indexInfo().spectrumDefinitions().get(), defs.get());
  }

  void test_clone_shares_data_in_IndexInfo() {
    WorkspaceTester ws;
    ws.initialize(3, 1, 1);
    ws.setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    const auto clone = ws.clone();
    const auto &info1 = ws.indexInfo();
    const auto &info2 = clone->indexInfo();
    // Spectrum numbers should also be shared, but there is no access by
    // reference, so we cannot check.
    TS_ASSERT_EQUALS(info1.spectrumDefinitions(), info2.spectrumDefinitions());
    TS_ASSERT(info1.spectrumDefinitions()); // should not be nullptr
    TS_ASSERT_EQUALS(info1.spectrumDefinitions(), info2.spectrumDefinitions());
  }

  void test_WorkspaceFactory_shares_data_in_IndexInfo() {
    const auto ws =
        WorkspaceFactory::Instance().create("WorkspaceTester", 3, 1, 1);
    ws->setInstrument(
        ComponentCreationHelper::createTestInstrumentRectangular(1, 2));
    const auto copy = WorkspaceFactory::Instance().create(ws);
    const auto &info1 = ws->indexInfo();
    const auto &info2 = copy->indexInfo();
    // Spectrum numbers should also be shared, but there is no access by
    // reference, so we cannot check.
    TS_ASSERT_EQUALS(info1.spectrumDefinitions(), info2.spectrumDefinitions());
    TS_ASSERT(info1.spectrumDefinitions()); // should not be nullptr
    TS_ASSERT_EQUALS(info1.spectrumDefinitions(), info2.spectrumDefinitions());
  }

  void test_toString_Produces_Expected_Contents() {
    auto testWS = boost::make_shared<WorkspaceTester>();
    testWS->initialize(1, 2, 1);
    testWS->setTitle("A test run");
    testWS->getAxis(0)->setUnit("TOF");
    testWS->setYUnitLabel("Counts");

    std::string expected = "WorkspaceTester\n"
                           "Title: A test run\n"
                           "Histograms: 1\n"
                           "Bins: 1\n"
                           "Histogram\n"
                           "X axis: Time-of-flight / microsecond\n"
                           "Y axis: Counts\n"
                           "Distribution: False\n"
                           "Instrument: None\n"
                           "Run start: not available\n"
                           "Run end:  not available\n";

    TS_ASSERT_EQUALS(expected, testWS->toString());
  }

  void test_initialize_with_IndexInfo_does_not_set_default_detectorIDs() {
    WorkspaceTester ws;
    Indexing::IndexInfo indexInfo(1);
    ws.initialize(indexInfo,
                  HistogramData::Histogram(HistogramData::Points(1)));
    TS_ASSERT_EQUALS(ws.getSpectrum(0).getDetectorIDs().size(), 0);
  }

  void testCloneClearsWorkspaceName() {
    auto ws = boost::make_shared<WorkspaceTester>();
    ws->initialize(1, 1, 1);
    const std::string name{"MatrixWorkspace_testCloneClearsWorkspaceName"};
    AnalysisDataService::Instance().add(name, ws);
    TS_ASSERT_EQUALS(ws->getName(), name)
    auto cloned = ws->clone();
    TS_ASSERT(cloned->getName().empty())
    AnalysisDataService::Instance().clear();
  }

  void testGetSetTitle() {
    TS_ASSERT_EQUALS(ws->getTitle(), "");
    ws->setTitle("something");
    TS_ASSERT_EQUALS(ws->getTitle(), "something");
    ws->setTitle("");
  }

  void testGetSetComment() {
    TS_ASSERT_EQUALS(ws->getComment(), "");
    ws->setComment("commenting");
    TS_ASSERT_EQUALS(ws->getComment(), "commenting");
    ws->setComment("");
  }

  void test_getIndicesFromDetectorIDs() {
    WorkspaceTester ws;
    ws.initialize(10, 1, 1);
    for (size_t i = 0; i < 10; i++)
      ws.getSpectrum(i).setDetectorID(detid_t(i * 10));
    std::vector<detid_t> dets;
    dets.push_back(60);
    dets.push_back(20);
    dets.push_back(90);
    std::vector<size_t> indices = ws.getIndicesFromDetectorIDs(dets);
    TS_ASSERT_EQUALS(indices.size(), 3);
    TS_ASSERT_EQUALS(indices[0], 6);
    TS_ASSERT_EQUALS(indices[1], 2);
    TS_ASSERT_EQUALS(indices[2], 9);
  }

  void
  test_That_A_Workspace_Gets_SpectraMap_When_Initialized_With_NVector_Elements() {
    WorkspaceTester testWS;
    const size_t nhist(10);
    testWS.initialize(nhist, 1, 1);
    for (size_t i = 0; i < testWS.getNumberHistograms(); i++) {
      TS_ASSERT_EQUALS(testWS.getSpectrum(i).getSpectrumNo(), specnum_t(i + 1));
      TS_ASSERT(testWS.getSpectrum(i).hasDetectorID(detid_t(i)));
    }
  }

  void testEmptyWorkspace() {
    WorkspaceTester ws;
    TS_ASSERT(ws.isCommonBins());
    TS_ASSERT_EQUALS(ws.blocksize(), 0);
    TS_ASSERT_EQUALS(ws.size(), 0);
  }

  void test_updateSpectraUsing() {
    WorkspaceTester testWS;
    testWS.initialize(3, 1, 1);

    specnum_t specs[] = {1, 2, 2, 3};
    detid_t detids[] = {10, 99, 20, 30};
    TS_ASSERT_THROWS_NOTHING(
        testWS.updateSpectraUsing(SpectrumDetectorMapping(specs, detids, 4)));

    TS_ASSERT(testWS.getSpectrum(0).hasDetectorID(10));
    TS_ASSERT(testWS.getSpectrum(1).hasDetectorID(20));
    TS_ASSERT(testWS.getSpectrum(1).hasDetectorID(99));
    TS_ASSERT(testWS.getSpectrum(2).hasDetectorID(30));
  }

  void testDetectorMappingCopiedWhenAWorkspaceIsCopied() {
    boost::shared_ptr<MatrixWorkspace> parent =
        boost::make_shared<WorkspaceTester>();
    parent->initialize(1, 1, 1);
    parent->getSpectrum(0).setSpectrumNo(99);
    parent->getSpectrum(0).setDetectorID(999);

    MatrixWorkspace_sptr copied = WorkspaceFactory::Instance().create(parent);
    // Has it been copied?
    TS_ASSERT_EQUALS(copied->getSpectrum(0).getSpectrumNo(), 99);
    TS_ASSERT(copied->getSpectrum(0).hasDetectorID(999));
  }

  void testGetMemorySize() { TS_ASSERT_THROWS_NOTHING(ws->getMemorySize()); }

  void testHistory() { TS_ASSERT_THROWS_NOTHING(ws->history()); }

  void testAxes() { TS_ASSERT_EQUALS(ws->axes(), 2); }

  void testGetAxis() {
    TS_ASSERT_THROWS(ws->getAxis(-1), Exception::IndexError);
    TS_ASSERT_THROWS_NOTHING(ws->getAxis(0));
    TS_ASSERT(ws->getAxis(0));
    TS_ASSERT(ws->getAxis(0)->isNumeric());
    TS_ASSERT_THROWS(ws->getAxis(2), Exception::IndexError);
  }

  void testReplaceAxis() {
    Axis *ax = new SpectraAxis(ws.get());
    TS_ASSERT_THROWS(ws->replaceAxis(2, ax), Exception::IndexError);
    TS_ASSERT_THROWS_NOTHING(ws->replaceAxis(0, ax));
    TS_ASSERT(ws->getAxis(0)->isSpectra());
  }

  void testIsDistribution() {
    TS_ASSERT(!ws->isDistribution());
    ws->setDistribution(true);
    TS_ASSERT(ws->isDistribution());
  }

  void testGetSetYUnit() {
    TS_ASSERT_EQUALS(ws->YUnit(), "");
    TS_ASSERT_THROWS_NOTHING(ws->setYUnit("something"));
    TS_ASSERT_EQUALS(ws->YUnit(), "something");
  }

  void testGetSpectrum() {
    WorkspaceTester ws;
    ws.initialize(4, 1, 1);
    TS_ASSERT_THROWS_NOTHING(ws.getSpectrum(0));
    TS_ASSERT_THROWS_NOTHING(ws.getSpectrum(3));
  }

  /** Get a detector sptr for each spectrum */
  void testGetDetector() {
    // Workspace has 3 spectra, each 1 in length
    const int numHist(3);
    boost::shared_ptr<MatrixWorkspace> workspace(
        makeWorkspaceWithDetectors(3, 1));

    // Initially un masked
    for (int i = 0; i < numHist; ++i) {
      IDetector_const_sptr det;
      TS_ASSERT_THROWS_NOTHING(det = workspace->getDetector(i));
      if (det) {
        TS_ASSERT_EQUALS(det->getID(), i);
      } else {
        TS_FAIL("No detector defined");
      }
    }

    // Now a detector group
    auto &spec = workspace->getSpectrum(0);
    spec.addDetectorID(1);
    spec.addDetectorID(2);
    IDetector_const_sptr det;
    TS_ASSERT_THROWS_NOTHING(det = workspace->getDetector(0));
    TS_ASSERT(det);

    // Now an empty (no detector) pixel
    auto &spec2 = workspace->getSpectrum(1);
    spec2.clearDetectorIDs();
    IDetector_const_sptr det2;
    TS_ASSERT_THROWS_ANYTHING(det2 = workspace->getDetector(1));
    TS_ASSERT(!det2);
  }

  void testWholeSpectraMasking() {
    // Workspace has 3 spectra, each 1 in length
    const int numHist(3);
    boost::shared_ptr<MatrixWorkspace> workspace(
        makeWorkspaceWithDetectors(3, 1));

    // Initially un masked
    const auto &spectrumInfo = workspace->spectrumInfo();
    for (int i = 0; i < numHist; ++i) {
      TS_ASSERT_EQUALS(workspace->readY(i)[0], 1.0);
      TS_ASSERT_EQUALS(workspace->readE(i)[0], 1.0);
      TS_ASSERT(spectrumInfo.hasDetectors(i));
      TS_ASSERT_EQUALS(spectrumInfo.isMasked(i), false);
    }

    // Mask a spectra
    workspace->getSpectrum(1).clearData();
    workspace->getSpectrum(2).clearData();
    workspace->mutableSpectrumInfo().setMasked(1, true);
    workspace->mutableSpectrumInfo().setMasked(2, true);

    const auto &spectrumInfo2 = workspace->spectrumInfo();
    for (int i = 0; i < numHist; ++i) {
      double expectedValue(0.0);
      bool expectedMasked(false);
      if (i == 0) {
        expectedValue = 1.0;
        expectedMasked = false;
      } else {
        expectedMasked = true;
      }
      TS_ASSERT_EQUALS(workspace->readY(i)[0], expectedValue);
      TS_ASSERT_EQUALS(workspace->readE(i)[0], expectedValue);
      TS_ASSERT(spectrumInfo2.hasDetectors(i));
      TS_ASSERT_EQUALS(spectrumInfo2.isMasked(i), expectedMasked);
    }
  }

  void testWholeSpectraMasking_SpectrumInfo() {
    // Workspace has 3 spectra, each 1 in length
    const int numHist(3);
    auto workspace = makeWorkspaceWithDetectors(numHist, 1);
    workspace->getSpectrum(1).clearData();
    workspace->getSpectrum(2).clearData();
    workspace->mutableSpectrumInfo().setMasked(1, true);
    workspace->mutableSpectrumInfo().setMasked(2, true);

    const auto &spectrumInfo = workspace->spectrumInfo();
    for (int i = 0; i < numHist; ++i) {
      bool expectedMasked(false);
      if (i == 0) {
        expectedMasked = false;
      } else {
        expectedMasked = true;
      }
      TS_ASSERT_EQUALS(spectrumInfo.isMasked(i), expectedMasked);
    }
  }

  void test_spectrumInfo_works_unthreaded() {
    const int numHist(3);
    auto workspace = makeWorkspaceWithDetectors(numHist, 1);
    std::atomic<bool> parallelException{false};
    for (int i = 0; i < numHist; ++i) {
      try {
        static_cast<void>(workspace->spectrumInfo());
      } catch (...) {
        parallelException = true;
      }
    }
    TS_ASSERT(!parallelException);
  }

  void test_spectrumInfo_works_threaded() {
    const int numHist(3);
    auto workspace = makeWorkspaceWithDetectors(numHist, 1);
    std::vector<const SpectrumInfo *> spectrumInfos(numHist);
    std::atomic<bool> parallelException{false};
    PARALLEL_FOR_IF(Kernel::threadSafe(*workspace))
    for (int i = 0; i < numHist; ++i) {
      try {
        spectrumInfos[i] = &(workspace->spectrumInfo());
      } catch (...) {
        parallelException = true;
      }
    }
    TS_ASSERT(!parallelException);
    for (int i = 0; i < numHist; ++i)
      TS_ASSERT_EQUALS(spectrumInfos[0], spectrumInfos[i]);
  }

  void testFlagMasked() {
    auto ws = makeWorkspaceWithDetectors(2, 2);
    // Now do a valid masking
    TS_ASSERT_THROWS_NOTHING(ws->flagMasked(0, 1, 0.75));
    TS_ASSERT(ws->hasMaskedBins(0));
    TS_ASSERT_EQUALS(ws->maskedBins(0).size(), 1);
    TS_ASSERT_EQUALS(ws->maskedBins(0).begin()->first, 1);
    TS_ASSERT_EQUALS(ws->maskedBins(0).begin()->second, 0.75);
    // flagMasked() shouldn't change the y-value maskBins() tested below does
    // that
    TS_ASSERT_EQUALS(ws->dataY(0)[1], 1.0);

    // Now mask a bin earlier than above and check it's sorting properly
    TS_ASSERT_THROWS_NOTHING(ws->flagMasked(1, 1))
    TS_ASSERT_EQUALS(ws->maskedBins(1).size(), 1)
    TS_ASSERT_EQUALS(ws->maskedBins(1).begin()->first, 1)
    TS_ASSERT_EQUALS(ws->maskedBins(1).begin()->second, 1.0)
    // Check the previous masking is still OK
    TS_ASSERT_EQUALS(ws->maskedBins(0).rbegin()->first, 1)
    TS_ASSERT_EQUALS(ws->maskedBins(0).rbegin()->second, 0.75)
  }

  void testMasking() {
    auto ws2 = makeWorkspaceWithDetectors(1, 2);

    TS_ASSERT(!ws2->hasMaskedBins(0));
    // Doesn't throw on invalid spectrum number, just returns false
    TS_ASSERT(!ws2->hasMaskedBins(1));
    TS_ASSERT(!ws2->hasMaskedBins(-1));

    // Will throw if nothing masked for spectrum
    TS_ASSERT_THROWS(ws2->maskedBins(0), Mantid::Kernel::Exception::IndexError);
    // Will throw if attempting to mask invalid spectrum
    TS_ASSERT_THROWS(ws2->maskBin(-1, 1),
                     Mantid::Kernel::Exception::IndexError);
    TS_ASSERT_THROWS(ws2->maskBin(1, 1), Mantid::Kernel::Exception::IndexError);
    // ...or an invalid bin
    TS_ASSERT_THROWS(ws2->maskBin(0, -1),
                     Mantid::Kernel::Exception::IndexError);
    TS_ASSERT_THROWS(ws2->maskBin(0, 2), Mantid::Kernel::Exception::IndexError);

    // Now do a valid masking
    TS_ASSERT_THROWS_NOTHING(ws2->maskBin(0, 1, 0.5));
    TS_ASSERT(ws2->hasMaskedBins(0));
    TS_ASSERT_EQUALS(ws2->maskedBins(0).size(), 1);
    TS_ASSERT_EQUALS(ws2->maskedBins(0).begin()->first, 1);
    TS_ASSERT_EQUALS(ws2->maskedBins(0).begin()->second, 0.5);
    TS_ASSERT_EQUALS(ws2->dataY(0)[1], 0.5);

    // Now mask a bin earlier than above and check it's sorting properly
    TS_ASSERT_THROWS_NOTHING(ws2->maskBin(0, 0));
    TS_ASSERT_EQUALS(ws2->maskedBins(0).begin()->first, 0);
    TS_ASSERT_EQUALS(ws2->maskedBins(0).begin()->second, 1.0);
    TS_ASSERT_EQUALS(ws2->dataY(0)[0], 0.0);
    // Check the previous masking is still OK
    TS_ASSERT_EQUALS(ws2->maskedBins(0).rbegin()->first, 1);
    TS_ASSERT_EQUALS(ws2->maskedBins(0).rbegin()->second, 0.5);
    TS_ASSERT_EQUALS(ws2->dataY(0)[1], 0.5);
  }

  void testMaskingNaNInf() {
    const size_t s = 4;
    const double y[s] = {NAN, INFINITY, -INFINITY, 2.};
    WorkspaceTester ws;
    ws.initialize(1, s + 1, s);

    // initialize and mask first with 0 weights
    // masking with 0 weight should be equiavalent to flagMasked
    // i.e. values should not change, even Inf and NaN
    for (size_t i = 0; i < s; ++i) {
      ws.mutableY(0)[i] = y[i];
      ws.maskBin(0, i, 0);
    }

    TS_ASSERT(std::isnan(ws.y(0)[0]));
    TS_ASSERT(std::isinf(ws.y(0)[1]));
    TS_ASSERT(std::isinf(ws.y(0)[2]));
    TS_ASSERT_EQUALS(ws.y(0)[3], 2.);

    // now mask w/o specifying weight (e.g. 1 by default)
    // in this case everything should be 0, even NaN and Inf
    for (size_t i = 0; i < s; ++i) {
      ws.maskBin(0, i);
      TS_ASSERT_EQUALS(ws.y(0)[i], 0.);
    }
  }

  void testSetMaskedBins() {
    auto ws = makeWorkspaceWithDetectors(2, 2);
    ws->flagMasked(0, 1);
    ws->flagMasked(1, 0);
    ws->setMaskedBins(1, ws->maskedBins(0));
    TS_ASSERT(ws->hasMaskedBins(1));
    TS_ASSERT_EQUALS(ws->maskedBins(1).size(), 1);
    TS_ASSERT_EQUALS(ws->maskedBins(0).begin()->first, 1);
  }

  void testSize() {
    WorkspaceTester wkspace;
    wkspace.initialize(1, 4, 3);
    TS_ASSERT_EQUALS(wkspace.blocksize(), 3);
    TS_ASSERT_EQUALS(wkspace.size(), 3);
  }

  void testBinIndexOf() {
    WorkspaceTester wkspace;
    wkspace.initialize(1, 4, 3);
    // Data is all 1.0s
    wkspace.dataX(0)[1] = 2.0;
    wkspace.dataX(0)[2] = 3.0;
    wkspace.dataX(0)[3] = 4.0;

    TS_ASSERT_EQUALS(wkspace.getNumberHistograms(), 1);

    // First bin
    TS_ASSERT_EQUALS(wkspace.binIndexOf(1.3), 0);
    // Bin boundary
    TS_ASSERT_EQUALS(wkspace.binIndexOf(2.0), 0);
    // Mid range
    TS_ASSERT_EQUALS(wkspace.binIndexOf(2.5), 1);
    // Still second bin
    TS_ASSERT_EQUALS(wkspace.binIndexOf(2.001), 1);
    // Last bin
    TS_ASSERT_EQUALS(wkspace.binIndexOf(3.1), 2);
    // Last value
    TS_ASSERT_EQUALS(wkspace.binIndexOf(4.0), 2);

    // Error handling

    // Bad index value
    TS_ASSERT_THROWS(wkspace.binIndexOf(2.5, 1), std::out_of_range);
    TS_ASSERT_THROWS(wkspace.binIndexOf(2.5, -1), std::out_of_range);

    // Bad X values
    TS_ASSERT_THROWS(wkspace.binIndexOf(5.), std::out_of_range);
    TS_ASSERT_THROWS(wkspace.binIndexOf(0.), std::out_of_range);
  }

  void testBinIndexOfDescendingBinning() {
    WorkspaceTester wkspace;
    wkspace.initialize(1, 4, 3);

    wkspace.dataX(0)[0] = 5.3;
    wkspace.dataX(0)[1] = 4.3;
    wkspace.dataX(0)[2] = 3.3;
    wkspace.dataX(0)[3] = 2.3;

    TS_ASSERT_EQUALS(wkspace.getNumberHistograms(), 1);

    // First boundary
    TS_ASSERT_EQUALS(wkspace.binIndexOf(5.3), 0)
    // First bin
    TS_ASSERT_EQUALS(wkspace.binIndexOf(5.2), 0);
    // Bin boundary
    TS_ASSERT_EQUALS(wkspace.binIndexOf(4.3), 0);
    // Mid range
    TS_ASSERT_EQUALS(wkspace.binIndexOf(3.8), 1);
    // Still second bin
    TS_ASSERT_EQUALS(wkspace.binIndexOf(std::nextafter(3.3, 10.0)), 1);
    // Last bin
    TS_ASSERT_EQUALS(wkspace.binIndexOf(3.1), 2);
    // Last value
    TS_ASSERT_EQUALS(wkspace.binIndexOf(2.3), 2);

    // Error handling

    // Bad index value
    TS_ASSERT_THROWS(wkspace.binIndexOf(2.5, 1), std::out_of_range);
    TS_ASSERT_THROWS(wkspace.binIndexOf(2.5, -1), std::out_of_range);

    // Bad X values
    TS_ASSERT_THROWS(wkspace.binIndexOf(std::nextafter(5.3, 10.0)),
                     std::out_of_range);
    TS_ASSERT_THROWS(wkspace.binIndexOf(5.4), std::out_of_range);
    TS_ASSERT_THROWS(wkspace.binIndexOf(std::nextafter(2.3, 0.0)),
                     std::out_of_range);
    TS_ASSERT_THROWS(wkspace.binIndexOf(0.), std::out_of_range);
  }

  void test_hasGroupedDetectors() {
    auto ws = makeWorkspaceWithDetectors(5, 1);
    TS_ASSERT_EQUALS(ws->hasGroupedDetectors(), false);

    ws->getSpectrum(0).addDetectorID(3);
    TS_ASSERT_EQUALS(ws->hasGroupedDetectors(), true);
  }

  void test_getSpectrumToWorkspaceIndexMap() {
    WorkspaceTester ws;
    ws.initialize(2, 1, 1);
    auto map = ws.getSpectrumToWorkspaceIndexMap();
    TS_ASSERT_EQUALS(0, map[1]);
    TS_ASSERT_EQUALS(1, map[2]);
    TS_ASSERT_EQUALS(map.size(), 2);

    // Check it throws for non-spectra axis
    ws.replaceAxis(1, new NumericAxis(1));
    TS_ASSERT_THROWS(ws.getSpectrumToWorkspaceIndexMap(), std::runtime_error);
  }

  void test_getDetectorIDToWorkspaceIndexMap() {
    auto ws = makeWorkspaceWithDetectors(5, 1);
    detid2index_map idmap = ws->getDetectorIDToWorkspaceIndexMap(true);

    TS_ASSERT_EQUALS(idmap.size(), 5);
    int i = 0;
    for (auto it = idmap.begin(); it != idmap.end(); ++it, ++i) {
      TS_ASSERT_EQUALS(idmap.count(i), 1);
      TS_ASSERT_EQUALS(idmap[i], i);
    }

    ws->getSpectrum(2).addDetectorID(99); // Set a second ID on one spectrum
    TS_ASSERT_THROWS(ws->getDetectorIDToWorkspaceIndexMap(true),
                     std::runtime_error);
    detid2index_map idmap2 = ws->getDetectorIDToWorkspaceIndexMap();
    TS_ASSERT_EQUALS(idmap2.size(), 6);
  }

  void test_getDetectorIDToWorkspaceIndexVector() {
    auto ws = makeWorkspaceWithDetectors(100, 10);
    std::vector<size_t> out;
    detid_t offset = -1234;
    TS_ASSERT_THROWS_NOTHING(
        out = ws->getDetectorIDToWorkspaceIndexVector(offset));
    TS_ASSERT_EQUALS(offset, 0);
    TS_ASSERT_EQUALS(out.size(), 100);
    TS_ASSERT_EQUALS(out[0], 0);
    TS_ASSERT_EQUALS(out[1], 1);
    TS_ASSERT_EQUALS(out[99], 99);

    // Create some discontinuities and check that the default value is there
    // Have to create a whole new instrument to keep things consistent, since
    // the detector ID
    // is stored in at least 3 places
    auto inst = boost::make_shared<Instrument>("TestInstrument");
    // We get a 1:1 map by default so the detector ID should match the spectrum
    // number
    for (size_t i = 0; i < ws->getNumberHistograms(); ++i) {
      detid_t detid = static_cast<detid_t>(i);
      // Create a detector for each spectra
      if (i == 0)
        detid = -1;
      if (i == 99)
        detid = 110;
      Detector *det = new Detector("pixel", detid, inst.get());
      inst->add(det);
      inst->markAsDetector(det);
      ws->getSpectrum(i).addDetectorID(detid);
    }
    ws->setInstrument(inst);
    ws->getSpectrum(66).clearDetectorIDs();

    TS_ASSERT_THROWS_NOTHING(
        out = ws->getDetectorIDToWorkspaceIndexVector(offset));
    TS_ASSERT_EQUALS(offset, 1);
    TS_ASSERT_EQUALS(out.size(), 112);
    TS_ASSERT_EQUALS(out[66 + offset], std::numeric_limits<size_t>::max());
    TS_ASSERT_EQUALS(out[99 + offset], 99);
    TS_ASSERT_EQUALS(out[105 + offset], std::numeric_limits<size_t>::max());
    TS_ASSERT_EQUALS(out[110 + offset], 99);
  }

  void test_getSpectrumToWorkspaceIndexVector() {
    auto ws = makeWorkspaceWithDetectors(100, 10);
    std::vector<size_t> out;
    detid_t offset = -1234;
    TS_ASSERT_THROWS_NOTHING(out =
                                 ws->getSpectrumToWorkspaceIndexVector(offset));
    TS_ASSERT_EQUALS(offset, -1);
    TS_ASSERT_EQUALS(out.size(), 100);
    TS_ASSERT_EQUALS(out[0], 0);
    TS_ASSERT_EQUALS(out[1], 1);
    TS_ASSERT_EQUALS(out[99], 99);
  }

  void test_getSignalAtCoord_histoData() {
    // Create a test workspace
    const auto ws = createTestWorkspace(4, 6, 5);

    // Get signal at coordinates
    std::vector<coord_t> coords = {0.5, 1.0};
    TS_ASSERT_DELTA(
        ws.getSignalAtCoord(coords.data(), Mantid::API::NoNormalization), 0.0,
        1e-5);
    coords[0] = 1.5;
    TS_ASSERT_DELTA(
        ws.getSignalAtCoord(coords.data(), Mantid::API::NoNormalization), 1.0,
        1e-5);
  }

  void test_getSignalAtCoord_pointData() {
    // Create a test workspace
    const auto ws = createTestWorkspace(4, 5, 5);
    auto normType = Mantid::API::NoNormalization;

    // Get signal at coordinates
    std::vector<coord_t> coords = {-1.0, 1.0};
    coords[0] = -0.75;
    TS_ASSERT(std::isnan(ws.getSignalAtCoord(coords.data(), normType)));
    coords[0] = -0.25;
    TS_ASSERT_DELTA(ws.getSignalAtCoord(coords.data(), normType), 0.0, 1e-5);
    coords[0] = 0.0;
    TS_ASSERT_DELTA(ws.getSignalAtCoord(coords.data(), normType), 0.0, 1e-5);
    coords[0] = 0.25;
    TS_ASSERT_DELTA(ws.getSignalAtCoord(coords.data(), normType), 0.0, 1e-5);
    coords[0] = 0.75;
    TS_ASSERT_DELTA(ws.getSignalAtCoord(coords.data(), normType), 1.0, 1e-5);
    coords[0] = 1.0;
    TS_ASSERT_DELTA(ws.getSignalAtCoord(coords.data(), normType), 1.0, 1e-5);
    coords[0] = 4.25;
    TS_ASSERT_DELTA(ws.getSignalAtCoord(coords.data(), normType), 4.0, 1e-5);
    coords[0] = 4.75;
    TS_ASSERT(std::isnan(ws.getSignalAtCoord(coords.data(), normType)));
  }

  void test_getCoordAtSignal_regression() {
    /*
    Having more spectrum numbers (acutally vertical axis increments) than x bins
    in VolumeNormalisation mode
    should not cause any issues.
    */
    WorkspaceTester ws;
    const int nVertical = 4;

    const int nBins = 2;
    const int nYValues = 1;
    ws.initialize(nVertical, nBins, nYValues);
    NumericAxis *verticalAxis = new NumericAxis(nVertical);
    for (int i = 0; i < nVertical; ++i) {
      for (int j = 0; j < nBins; ++j) {
        if (j < nYValues) {
          ws.dataY(i)[j] = 1.0; // All y values are 1.
          ws.dataE(i)[j] = j;
        }
        ws.dataX(i)[j] = j; // x increments by 1
      }
      verticalAxis->setValue(i, double(i)); // Vertical axis increments by 1.
    }
    ws.replaceAxis(1, verticalAxis);
    // Signal is always 1 and volume of each box is 1. Therefore normalized