MatrixWorkspace.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 "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/MatrixWorkspaceMDIterator.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/MDGeometry/GeneralFrame.h"
#include "MantidGeometry/MDGeometry/MDFrame.h"
#include "MantidIndexing/GlobalSpectrumIndex.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/MDUnit.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/VectorHelper.h"
 
#include "MantidParallel/Collectives.h"
#include "MantidParallel/Communicator.h"
#include "MantidTypes/SpectrumDefinition.h"
 
#include <boost/algorithm/string/join.hpp>
#include <boost/algorithm/string/regex.hpp>
 
#include <cmath>
#include <functional>
#include <numeric>
#include <utility>
 
using Mantid::Kernel::TimeSeriesProperty;
using Mantid::Types::Core::DateAndTime;
 
namespace {
auto accumulate_if_finite = [](const double accumulator, const double newValue) {
  if (std::isfinite(newValue)) {
    return accumulator + newValue;
  } else {
    return accumulator;
  }
};
} // namespace
 
namespace Mantid::API {
using std::size_t;
using namespace Geometry;
using Kernel::V3D;
 
namespace {
Kernel::Logger g_log("MatrixWorkspace");
constexpr const double EPSILON{1.0e-9};
 
std::string appendUnitDenominatorUsingPer(std::string yLabel, const MatrixWorkspace &workspace, bool useLatex) {
  if (useLatex) {
    std::string xLabel = workspace.getAxis(0)->unit()->label().latex();
    if (!xLabel.empty())
      yLabel += " per $" + workspace.getAxis(0)->unit()->label().latex() + "$";
  } else {
    std::string xLabel = workspace.getAxis(0)->unit()->label().ascii();
    if (!xLabel.empty())
      yLabel += " per " + workspace.getAxis(0)->unit()->label().ascii();
  }
  return yLabel;
}
 
std::string replacePerWithLatex(std::string yLabel) {
  std::vector<std::string> splitVec;
  boost::split_regex(splitVec, yLabel, boost::regex(" per "));
  if (splitVec.size() > 1) {
    yLabel = splitVec[0];
    splitVec.erase(splitVec.begin());
    std::string unitString = boost::algorithm::join(splitVec, " ");
    if (!yLabel.empty())
      yLabel += " ";
    yLabel += "(" + unitString + ")$^{-1}$";
  }
  return yLabel;
}
 
} // namespace
const std::string MatrixWorkspace::xDimensionId = "xDimension";
const std::string MatrixWorkspace::yDimensionId = "yDimension";
 
MatrixWorkspace::MatrixWorkspace(const Parallel::StorageMode storageMode)
    : IMDWorkspace(storageMode), ExperimentInfo(), m_axes(), m_isInitialized(false), m_YUnit(), m_YUnitLabel(),
      m_masks() {}
 
MatrixWorkspace::MatrixWorkspace(const MatrixWorkspace &other)
    : IMDWorkspace(other), ExperimentInfo(other), m_indexInfo(std::make_unique<Indexing::IndexInfo>(other.indexInfo())),
      m_isInitialized(other.m_isInitialized), m_YUnit(other.m_YUnit), m_YUnitLabel(other.m_YUnitLabel),
      m_isCommonBinsFlag(other.m_isCommonBinsFlag), m_masks(other.m_masks), m_indexInfoNeedsUpdate(false) {
  m_axes.resize(other.m_axes.size());
  for (size_t i = 0; i < m_axes.size(); ++i)
    m_axes[i] = std::unique_ptr<Axis>(other.m_axes[i]->clone(this));
  m_isCommonBinsFlagValid.store(other.m_isCommonBinsFlagValid.load());
  // TODO: Do we need to init m_monitorWorkspace?
}
 
// RJT, 3/10/07: The Analysis Data Service needs to be able to delete
// workspaces, so I moved this from protected to public.
MatrixWorkspace::~MatrixWorkspace() = default;
 
const Indexing::IndexInfo &MatrixWorkspace::indexInfo() const {
  std::lock_guard<std::mutex> lock(m_indexInfoMutex);
  // Individual SpectrumDefinitions in SpectrumInfo may have changed. Due to a
  // copy-on-write mechanism the definitions stored in IndexInfo may then be out
  // of sync (definitions in SpectrumInfo have been updated).
  m_indexInfo->setSpectrumDefinitions(spectrumInfo().sharedSpectrumDefinitions());
  // If spectrum numbers are set in ISpectrum this flag will be true. However,
  // for MPI builds we will forbid changing spectrum numbers in this way, so we
  // should never enter this branch. Thus it is sufficient to set only the local
  // spectrum numbers here.
  if (m_indexInfoNeedsUpdate) {
    std::vector<Indexing::SpectrumNumber> spectrumNumbers;
    for (size_t i = 0; i < getNumberHistograms(); ++i)
      spectrumNumbers.emplace_back(getSpectrum(i).getSpectrumNo());
    m_indexInfo->setSpectrumNumbers(std::move(spectrumNumbers));
    m_indexInfoNeedsUpdate = false;
  }
  return *m_indexInfo;
}
 
void MatrixWorkspace::setIndexInfo(const Indexing::IndexInfo &indexInfo) {
  if (indexInfo.storageMode() != storageMode())
    throw std::invalid_argument("MatrixWorkspace::setIndexInfo: "
                                "Parallel::StorageMode in IndexInfo does not "
                                "match storage mode in workspace");
 
  // Comparing the *local* size of the indexInfo.
  if (indexInfo.size() != getNumberHistograms())
    throw std::invalid_argument("MatrixWorkspace::setIndexInfo: IndexInfo size "
                                "does not match number of histograms in "
                                "workspace");
 
  m_indexInfo = std::make_unique<Indexing::IndexInfo>(indexInfo);
  m_indexInfoNeedsUpdate = false;
  if (!m_indexInfo->spectrumDefinitions())
    buildDefaultSpectrumDefinitions();
  // Fails if spectrum definitions contain invalid indices.
  rebuildDetectorIDGroupings();
  // This sets the SpectrumDefinitions for the SpectrumInfo, which may seem
  // counterintuitive at first -- why would setting IndexInfo modify internals
  // of SpectrumInfo? However, logically it would not make sense to assign
  // SpectrumDefinitions in an assignment of SpectrumInfo: Changing
  // SpectrumDefinitions requires also changes at a higher level of a workspace
  // (in particular the histograms, which would need to be regrouped as well).
  // Thus, assignment of SpectrumInfo should just check for compatible
  // SpectrumDefinitions and assign other data (such as per-spectrum masking
  // flags, which do not exist yet). Furthermore, since currently detector
  // groupings are still stored in ISpectrum (in addition to the
  // SpectrumDefinitions in SpectrumInfo), an assigment of SpectrumDefinitions
  // in SpectrumInfo would lead to inconsistent workspaces. SpectrumDefinitions
  // are thus assigned by IndexInfo, which acts at a highler level and is
  // typically used at construction time of a workspace, i.e., there is no data
  // in histograms yet which would need to be regrouped.
  setSpectrumDefinitions(m_indexInfo->spectrumDefinitions());
}
 
void MatrixWorkspace::setIndexInfoWithoutISpectrumUpdate(const Indexing::IndexInfo &indexInfo) {
  // Workspace is already initialized (m_isInitialized == true), but this is
  // called by initializedFromParent which is some sort of second-stage
  // initialization, so there is no check for storage mode compatibility here,
  // in contrast to what setIndexInfo() does.
  setStorageMode(indexInfo.storageMode());
  // Comparing the *local* size of the indexInfo.
  if (indexInfo.size() != getNumberHistograms())
    throw std::invalid_argument("MatrixWorkspace::setIndexInfo: IndexInfo size "
                                "does not match number of histograms in "
                                "workspace");
  *m_indexInfo = indexInfo;
  m_indexInfoNeedsUpdate = false;
  setSpectrumDefinitions(m_indexInfo->spectrumDefinitions());
}
 
const std::string MatrixWorkspace::toString() const {
  std::ostringstream os;
  os << id() << "\n"
     << "Title: " << getTitle() << "\n"
     << "Histograms: " << getNumberHistograms() << "\n"
     << "Bins: ";
 
  try {
    os << blocksize() << "\n";
  } catch (std::length_error &) {
    os << "variable\n"; // TODO shouldn't use try/catch
  }
 
  if (isHistogramData())
    os << "Histogram\n";
  else
    os << "Data points\n";
 
  os << "X axis: ";
  if (axes() > 0) {
    Axis *ax = getAxis(0);
    if (ax && ax->unit())
      os << ax->unit()->caption() << " / " << ax->unit()->label().ascii();
    else
      os << "Not set";
  } else {
    os << "N/A";
  }
  os << "\n"
     << "Y axis: " << YUnitLabel() << "\n";
 
  os << "Distribution: " << (isDistribution() ? "True" : "False") << "\n";
 
  os << ExperimentInfo::toString();
  return os.str();
}
 
void MatrixWorkspace::initialize(const std::size_t &NVectors, const std::size_t &XLength, const std::size_t &YLength) {
  // Check validity of arguments
  if (NVectors == 0 || XLength == 0 || YLength == 0) {
    throw std::out_of_range("All arguments to init must be positive and non-zero");
  }
 
  // Bypass the initialization if the workspace has already been initialized.
  if (m_isInitialized)
    return;
 
  setNumberOfDetectorGroups(NVectors);
  m_indexInfo = std::make_unique<Indexing::IndexInfo>(NVectors);
 
  // Invoke init() method of the derived class inside a try/catch clause
  try {
    this->init(NVectors, XLength, YLength);
  } catch (std::runtime_error &) {
    throw;
  }
 
  // Indicate that this workspace has been initialized to prevent duplicate
  // attempts.
  m_isInitialized = true;
}
 
void MatrixWorkspace::initialize(const std::size_t &NVectors, const HistogramData::Histogram &histogram) {
  Indexing::IndexInfo indices(NVectors);
  // Empty SpectrumDefinitions to indicate no default mapping to detectors.
  indices.setSpectrumDefinitions(std::vector<SpectrumDefinition>(NVectors));
  return initialize(indices, histogram);
}
 
void MatrixWorkspace::initialize(const Indexing::IndexInfo &indexInfo, const HistogramData::Histogram &histogram) {
  // Check validity of arguments
  if (indexInfo.size() == 0 || histogram.x().empty()) {
    throw std::out_of_range("All arguments to init must be positive and non-zero");
  }
 
  // Bypass the initialization if the workspace has already been initialized.
  if (m_isInitialized)
    return;
  setStorageMode(indexInfo.storageMode());
  setNumberOfDetectorGroups(indexInfo.size());
  init(histogram);
  setIndexInfo(indexInfo);
 
  // Indicate that this workspace has been initialized to prevent duplicate
  // attempts.
  m_isInitialized = true;
}
 
//---------------------------------------------------------------------------------------
void MatrixWorkspace::setTitle(const std::string &t) {
  Workspace::setTitle(t);
 
  // A MatrixWorkspace contains uniquely one Run object, hence for this
  // workspace
  // keep the Run object run_title property the same as the workspace title
  Run &run = mutableRun();
  run.addProperty("run_title", t, true);
}
 
const std::string MatrixWorkspace::getTitle() const {
  if (run().hasProperty("run_title")) {
    std::string title = run().getProperty("run_title")->value();
    return title;
  } else
    return Workspace::getTitle();
}
 
void MatrixWorkspace::updateSpectraUsing(const SpectrumDetectorMapping &map) {
  for (size_t j = 0; j < getNumberHistograms(); ++j) {
    auto &spec = getSpectrum(j);
    try {
      if (map.indexIsSpecNumber())
        spec.setDetectorIDs(map.getDetectorIDsForSpectrumNo(spec.getSpectrumNo()));
      else
        spec.setDetectorIDs(map.getDetectorIDsForSpectrumIndex(j));
    } catch (std::out_of_range &e) {
      // Get here if the spectrum number is not in the map.
      spec.clearDetectorIDs();
      g_log.debug(e.what());
      g_log.debug() << "Spectrum number " << spec.getSpectrumNo() << " not in map.\n";
    }
  }
}
 
void MatrixWorkspace::rebuildSpectraMapping(const bool includeMonitors, const specnum_t specNumOffset) {
  if (sptr_instrument->nelements() == 0) {
    return;
  }
 
  std::vector<detid_t> pixelIDs = this->getInstrument()->getDetectorIDs(!includeMonitors);
 
  try {
    size_t index = 0;
    std::vector<detid_t>::const_iterator iend = pixelIDs.end();
    for (std::vector<detid_t>::const_iterator it = pixelIDs.begin(); it != iend; ++it) {
      // The detector ID
      const detid_t detId = *it;
      const auto specNo = specnum_t(index + specNumOffset);
 
      if (index < this->getNumberHistograms()) {
        auto &spec = getSpectrum(index);
        spec.setSpectrumNo(specNo);
        spec.setDetectorID(detId);
      }
 
      index++;
    }
 
  } catch (std::runtime_error &) {
    throw;
  }
}
 
spec2index_map MatrixWorkspace::getSpectrumToWorkspaceIndexMap() const {
  auto *ax = dynamic_cast<SpectraAxis *>(this->m_axes[1].get());
  if (!ax)
    throw std::runtime_error("MatrixWorkspace::getSpectrumToWorkspaceIndexMap: "
                             "axis[1] is not a SpectraAxis, so I cannot "
                             "generate a map.");
  try {
    return ax->getSpectraIndexMap();
  } catch (std::runtime_error &) {
    g_log.error() << "MatrixWorkspace::getSpectrumToWorkspaceIndexMap: no elements!";
    throw;
  }
}
 
std::vector<size_t> MatrixWorkspace::getSpectrumToWorkspaceIndexVector(specnum_t &offset) const {
  auto *ax = dynamic_cast<SpectraAxis *>(this->m_axes[1].get());
  if (!ax)
    throw std::runtime_error("MatrixWorkspace::getSpectrumToWorkspaceIndexMap: "
                             "axis[1] is not a SpectraAxis, so I cannot "
                             "generate a map.");
 
  // Find the min/max spectra IDs
  specnum_t min = std::numeric_limits<specnum_t>::max(); // So that any number
  // will be less than this
  specnum_t max = -std::numeric_limits<specnum_t>::max(); // So that any number
  // will be greater than
  // this
  size_t length = ax->length();
  for (size_t i = 0; i < length; i++) {
    specnum_t spec = ax->spectraNo(i);
    if (spec < min)
      min = spec;
    if (spec > max)
      max = spec;
  }
 
  // Offset so that the "min" value goes to index 0
  offset = -min;
 
  // Size correctly
  std::vector<size_t> out(max - min + 1, 0);
 
  // Make the vector
  for (size_t i = 0; i < length; i++) {
    specnum_t spec = ax->spectraNo(i);
    out[spec + offset] = i;
  }
 
  return out;
}
 
bool MatrixWorkspace::hasGroupedDetectors() const {
  bool retVal = false;
 
  // Loop through the workspace index
  for (size_t workspaceIndex = 0; workspaceIndex < this->getNumberHistograms(); workspaceIndex++) {
    auto detList = getSpectrum(workspaceIndex).getDetectorIDs();
    if (detList.size() > 1) {
      retVal = true;
      break;
    }
  }
  return retVal;
}
 
detid2index_map MatrixWorkspace::getDetectorIDToWorkspaceIndexMap(bool throwIfMultipleDets,
                                                                  bool ignoreIfNoValidDets) const {
  detid2index_map map;
  const auto &specInfo = spectrumInfo();
 
  // Loop through the workspace index
  for (size_t workspaceIndex = 0; workspaceIndex < this->getNumberHistograms(); ++workspaceIndex) {
    // Workspaces can contain invalid detector IDs. hasDetectors will silently ignore them until this is fixed.
    if (ignoreIfNoValidDets && !specInfo.hasDetectors(workspaceIndex)) {
      continue;
    }
    auto detList = getSpectrum(workspaceIndex).getDetectorIDs();
 
    if (throwIfMultipleDets) {
      if (detList.size() > 1) {
        throw std::runtime_error("MatrixWorkspace::getDetectorIDToWorkspaceIndexMap(): more than 1 "
                                 "detector for one histogram! I cannot generate a map of detector "
                                 "ID to workspace index.");
      }
 
      // Set the KEY to the detector ID and the VALUE to the workspace index.
      if (detList.size() == 1)
        map[*detList.begin()] = workspaceIndex;
    } else {
      // Allow multiple detectors per workspace index
      for (auto det : detList)
        map[det] = workspaceIndex;
    }
 
    // Ignore if the detector list is empty.
  }
 
  return map;
}
 
std::vector<size_t> MatrixWorkspace::getDetectorIDToWorkspaceIndexVector(detid_t &offset,
                                                                         bool throwIfMultipleDets) const {
  // Make a correct initial size
  std::vector<size_t> out;
  detid_t minId = 0;
  detid_t maxId = 0;
  this->getInstrument()->getMinMaxDetectorIDs(minId, maxId);
  offset = -minId;
  const int outSize = maxId - minId + 1;
  // Allocate at once
  out.resize(outSize, std::numeric_limits<size_t>::max());
 
  for (size_t workspaceIndex = 0; workspaceIndex < getNumberHistograms(); ++workspaceIndex) {
    // Get the list of detectors from the WS index
    const auto &detList = this->getSpectrum(workspaceIndex).getDetectorIDs();
 
    if (throwIfMultipleDets && (detList.size() > 1))
      throw std::runtime_error("MatrixWorkspace::getDetectorIDToWorkspaceIndexVector(): more than 1 "
                               "detector for one histogram! I cannot generate a map of detector ID "
                               "to workspace index.");
 
    // Allow multiple detectors per workspace index, or,
    // If only one is allowed, then this has thrown already
    for (auto det : detList) {
      int index = det + offset;
      if (index < 0 || index >= outSize) {
        g_log.debug() << "MatrixWorkspace::getDetectorIDToWorkspaceIndexVector("
                         "): detector ID found ("
                      << det << " at workspace index " << workspaceIndex << ") is invalid.\n";
      } else
        // Save it at that point.
        out[index] = workspaceIndex;
    }
 
  } // (for each workspace index)
  return out;
}
 
std::vector<size_t> MatrixWorkspace::getIndicesFromSpectra(const std::vector<specnum_t> &spectraList) const {
  // Clear the output index list
  std::vector<size_t> indexList;
  indexList.reserve(this->getNumberHistograms());
 
  auto iter = spectraList.cbegin();
  while (iter != spectraList.cend()) {
    for (size_t i = 0; i < this->getNumberHistograms(); ++i) {
      if (this->getSpectrum(i).getSpectrumNo() == *iter) {
        indexList.emplace_back(i);
        break;
      }
    }
    ++iter;
  }
  return indexList;
}
 
size_t MatrixWorkspace::getIndexFromSpectrumNumber(const specnum_t specNo) const {
  for (size_t i = 0; i < this->getNumberHistograms(); ++i) {
    if (this->getSpectrum(i).getSpectrumNo() == specNo)
      return i;
  }
  throw std::runtime_error("Could not find spectrum number in any spectrum.");
}
 
std::vector<size_t> MatrixWorkspace::getIndicesFromDetectorIDs(const std::vector<detid_t> &detIdList) const {
  if (m_indexInfo->size() != m_indexInfo->globalSize())
    throw std::runtime_error("MatrixWorkspace: Using getIndicesFromDetectorIDs "
                             "in a parallel run is most likely incorrect. "
                             "Aborting.");
 
  std::map<detid_t, std::set<size_t>> detectorIDtoWSIndices;
  for (size_t i = 0; i < getNumberHistograms(); ++i) {
    auto detIDs = getSpectrum(i).getDetectorIDs();
    for (auto detID : detIDs) {
      detectorIDtoWSIndices[detID].insert(i);
    }
  }
 
  std::vector<size_t> indexList;
  indexList.reserve(detIdList.size());
  for (const auto detId : detIdList) {
    auto wsIndices = detectorIDtoWSIndices.find(detId);
    if (wsIndices != detectorIDtoWSIndices.end()) {
      std::copy(wsIndices->second.cbegin(), wsIndices->second.cend(), std::back_inserter(indexList));
    }
  }
  return indexList;
}
 
std::vector<specnum_t> MatrixWorkspace::getSpectraFromDetectorIDs(const std::vector<detid_t> &detIdList) const {
  std::vector<specnum_t> spectraList;
 
  // Try every detector in the list
  for (auto detId : detIdList) {
    bool foundDet = false;
    specnum_t foundSpecNum = 0;
 
    // Go through every histogram
    for (size_t i = 0; i < this->getNumberHistograms(); i++) {
      if (this->getSpectrum(i).hasDetectorID(detId)) {
        foundDet = true;
        foundSpecNum = this->getSpectrum(i).getSpectrumNo();
        break;
      }
    }
 
    if (foundDet)
      spectraList.emplace_back(foundSpecNum);
  } // for each detector ID in the list
  return spectraList;
}
 
double MatrixWorkspace::getXMin() const {
  double xmin;
  double xmax;
  this->getXMinMax(xmin, xmax); // delegate to the proper code
  return xmin;
}
 
double MatrixWorkspace::getXMax() const {
  double xmin;
  double xmax;
  this->getXMinMax(xmin, xmax); // delegate to the proper code
  return xmax;
}
 
void MatrixWorkspace::getXMinMax(double &xmin, double &xmax) const {
  // set to crazy values to start
  xmin = std::numeric_limits<double>::max();
  xmax = -1.0 * xmin;
  size_t numberOfSpectra = this->getNumberHistograms();
 
  // determine the data range
  for (size_t workspaceIndex = 0; workspaceIndex < numberOfSpectra; workspaceIndex++) {
    const auto &dataX = this->x(workspaceIndex);
    const double xfront = dataX.front();
    const double xback = dataX.back();
    if (std::isfinite(xfront) && std::isfinite(xback)) {
      if (xfront < xmin)
        xmin = xfront;
      if (xback > xmax)
        xmax = xback;
    }
  }
  if (m_indexInfo->size() != m_indexInfo->globalSize()) {
    auto &comm = m_indexInfo->communicator();
    std::vector<double> extrema(comm.size());
    Parallel::all_gather(comm, xmin, extrema);
    xmin = *std::min_element(extrema.begin(), extrema.end());
    Parallel::all_gather(comm, xmax, extrema);
    xmax = *std::max_element(extrema.begin(), extrema.end());
  }
}
 
void MatrixWorkspace::getIntegratedSpectra(std::vector<double> &out, const double minX, const double maxX,
                                           const bool entireRange) const {
  out.resize(this->getNumberHistograms(), 0.0);
 
  // offset for histogram data, because the x axis is not the same size for histogram and point data.
  const size_t histogramOffset = this->isHistogramData() ? 1 : 0;
 
  // Run in parallel if the implementation is threadsafe
  PARALLEL_FOR_IF(this->threadSafe())
  for (int wksp_index = 0; wksp_index < static_cast<int>(this->getNumberHistograms()); wksp_index++) {
    // Get Handle to data
    const Mantid::MantidVec &x = this->readX(wksp_index);
    const auto &y = this->y(wksp_index);
    // If it is a 1D workspace, no need to integrate
    if ((x.size() <= 1 + histogramOffset) && (!y.empty())) {
      out[wksp_index] = y[0];
    } else {
      // Iterators for limits - whole range by default
      Mantid::MantidVec::const_iterator lowit, highit;
      lowit = x.begin();
      highit = x.end() - histogramOffset;
 
      // But maybe we don't want the entire range?
      if (!entireRange) {
        // If the first element is lower that the xmin then search for new lowit
        if ((*lowit) < minX)
          lowit = std::lower_bound(x.begin(), x.end(), minX);
        // If the last element is higher that the xmax then search for new highit
        if (*(highit - 1 + histogramOffset) > maxX)
          highit = std::upper_bound(lowit, x.end(), maxX);
      }
 
      // Get the range for the y vector
      Mantid::MantidVec::difference_type distmin = std::distance(x.begin(), lowit);
      Mantid::MantidVec::difference_type distmax = std::distance(x.begin(), highit);
      double sum(0.0);
      if (distmin <= distmax) {
        // Integrate
        sum = std::accumulate(y.begin() + distmin, y.begin() + distmax, 0.0, accumulate_if_finite);
      }
      // Save it in the vector
      out[wksp_index] = sum;
    }
  }
}
 
Geometry::IDetector_const_sptr MatrixWorkspace::getDetector(const size_t workspaceIndex) const {
  const auto &dets = getSpectrum(workspaceIndex).getDetectorIDs();
  Instrument_const_sptr localInstrument = getInstrument();
  if (!localInstrument) {
    g_log.debug() << "No instrument defined.\n";
    throw Kernel::Exception::NotFoundError("Instrument not found", "");
  }
 
  const size_t ndets = dets.size();
  if (ndets == 1) {
    // If only 1 detector for the spectrum number, just return it
    return localInstrument->getDetector(*dets.begin());
  } else if (ndets == 0) {
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): No "
                                           "detectors for this workspace "
                                           "index.",
                                           "");
  }
  // Else need to construct a DetectorGroup and return that
  auto dets_ptr = localInstrument->getDetectors(dets);
  return std::make_shared<Geometry::DetectorGroup>(dets_ptr);
}
 
double MatrixWorkspace::detectorSignedTwoTheta(const Geometry::IDetector &det) const {
 
  Instrument_const_sptr instrument = getInstrument();
  Geometry::IComponent_const_sptr source = instrument->getSource();
  Geometry::IComponent_const_sptr sample = instrument->getSample();
  if (source == nullptr || sample == nullptr) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Instrument not sufficiently defined: failed to get source and/or "
        "sample");
  }
 
  const Kernel::V3D samplePos = sample->getPos();
  const Kernel::V3D beamLine = samplePos - source->getPos();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
  // Get the instrument up axis.
  const V3D &thetaSignAxis = instrument->getReferenceFrame()->vecThetaSign();
  return det.getSignedTwoTheta(samplePos, beamLine, thetaSignAxis);
}
 
double MatrixWorkspace::detectorTwoTheta(const Geometry::IDetector &det) const {
  Instrument_const_sptr instrument = this->getInstrument();
  Geometry::IComponent_const_sptr source = instrument->getSource();
  Geometry::IComponent_const_sptr sample = instrument->getSample();
  if (source == nullptr || sample == nullptr) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Instrument not sufficiently defined: failed to get source and/or "
        "sample");
  }
 
  const Kernel::V3D samplePos = sample->getPos();
  const Kernel::V3D beamLine = samplePos - source->getPos();
 
  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError("Source and sample are at same position!");
  }
  return det.getTwoTheta(samplePos, beamLine);
}
 
int MatrixWorkspace::axes() const { return static_cast<int>(m_axes.size()); }
 
Axis *MatrixWorkspace::getAxis(const std::size_t &axisIndex) const {
  if (axisIndex >= m_axes.size()) {
    throw Kernel::Exception::IndexError(axisIndex, m_axes.size(), "Argument to getAxis is invalid for this workspace");
  }
 
  return m_axes[axisIndex].get();
}
 
void MatrixWorkspace::replaceAxis(const std::size_t &axisIndex, std::unique_ptr<Axis> newAxis) {
  // First check that axisIndex is in range
  if (axisIndex >= m_axes.size()) {
    throw Kernel::Exception::IndexError(axisIndex, m_axes.size(), "Value of axisIndex is invalid for this workspace");
  }
  // If we're OK, then delete the old axis and set the pointer to the new one
  m_axes[axisIndex] = std::move(newAxis);
}
 
bool MatrixWorkspace::isCommonLogBins() const {
  if (!this->isCommonBins()) {
    return false;
  }
 
  if (this->getNumberHistograms() == 0) {
    return false;
  }
 
  const auto &x0 = this->x(0);
  if (x0.size() < 2) {
    return false;
  }
 
  // guard against all axis elements being equal
  if (x0[1] == x0[0]) {
    return false;
  }
 
  double diff = x0[1] / x0[0];
  if (!std::isfinite(diff)) {
    return false;
  }
  // ignore final bin, since it may be a different size
  for (size_t i = 1; i < x0.size() - 2; ++i) {
    if (std::isfinite(x0[i + 1]) && std::isfinite(x0[i])) {
      if (std::abs(x0[i + 1] / x0[i] - diff) > EPSILON) {
        return false;
      }
    } else {
      return false;
    }
  }
  return true;
}
 
size_t MatrixWorkspace::numberOfAxis() const { return m_axes.size(); }
 
std::string MatrixWorkspace::YUnit() const { return m_YUnit; }
 
void MatrixWorkspace::setYUnit(const std::string &newUnit) { m_YUnit = newUnit; }
 
std::string MatrixWorkspace::YUnitLabel(bool useLatex /* = false */, bool plotAsDistribution /* = false */) const {
  std::string retVal;
  if (!m_YUnitLabel.empty()) {
    retVal = m_YUnitLabel;
    // If a custom label has been set and we are dividing by bin width when
    // plotting (i.e. plotAsDistribution = true and the workspace is not a
    // distribution), we must append the x-unit as a divisor. We assume the
    // custom label contains the correct units for the data.
    if (plotAsDistribution && !this->isDistribution())
      retVal = appendUnitDenominatorUsingPer(retVal, *this, useLatex);
  } else {
    retVal = m_YUnit;
    // If no custom label is set and the workspace is a distribution we need to
    // append the divisor's unit to the label. If the workspace is not a
    // distribution, but we are plotting it as a distribution, we must append
    // the divisor's unit.
    if (plotAsDistribution || this->isDistribution())
      retVal = appendUnitDenominatorUsingPer(retVal, *this, useLatex);
  }
  if (useLatex)
    retVal = replacePerWithLatex(retVal);
  return retVal;
}
 
void MatrixWorkspace::setYUnitLabel(const std::string &newLabel) { m_YUnitLabel = newLabel; }
 
bool MatrixWorkspace::isDistribution() const {
  return getSpectrum(0).yMode() == HistogramData::Histogram::YMode::Frequencies;
}
 
void MatrixWorkspace::setDistribution(bool newValue) {
  if (isDistribution() == newValue)
    return;
  HistogramData::Histogram::YMode ymode =
      newValue ? HistogramData::Histogram::YMode::Frequencies : HistogramData::Histogram::YMode::Counts;
  for (size_t i = 0; i < getNumberHistograms(); ++i)
    getSpectrum(i).setYMode(ymode);
}
 
bool MatrixWorkspace::isHistogramData() const {
  // all spectra *should* have the same behavior
  return isHistogramDataByIndex(0);
}
 
bool MatrixWorkspace::isHistogramDataByIndex(const std::size_t index) const {
  bool isHist = (x(index).size() != y(index).size());
 
  // TODOHIST temporary sanity check
  if (isHist) {
    if (getSpectrum(index).histogram().xMode() != HistogramData::Histogram::XMode::BinEdges) {
      throw std::logic_error("In MatrixWorkspace::isHistogramData(): "
                             "Histogram::Xmode is not BinEdges");
    }
  } else {
    if (getSpectrum(index).histogram().xMode() != HistogramData::Histogram::XMode::Points) {
      throw std::logic_error("In MatrixWorkspace::isHistogramData(): "
                             "Histogram::Xmode is not Points");
    }
  }
  return isHist;
}
 
bool MatrixWorkspace::isCommonBins() const {
  std::lock_guard<std::mutex> lock{m_isCommonBinsMutex};
  const bool isFlagValid{m_isCommonBinsFlagValid.exchange(true)};
  if (isFlagValid) {
    return m_isCommonBinsFlag;
  }
  m_isCommonBinsFlag = true;
  const size_t numHist = this->getNumberHistograms();
  // there being only one or zero histograms is accepted as not being an error
  if (numHist <= 1) {
    return m_isCommonBinsFlag;
  }
 
  // First check if the x-axis shares a common ptr.
  const HistogramData::HistogramX *first = &x(0);
  for (size_t i = 1; i < numHist; ++i) {
    if (&x(i) != first) {
      m_isCommonBinsFlag = false;
      break;
    }
  }
 
  // If true, we may return here.
  if (m_isCommonBinsFlag) {
    return m_isCommonBinsFlag;
  }
 
  m_isCommonBinsFlag = true;
  // Check that that size of each histogram is identical.
  const size_t numBins = x(0).size();
  for (size_t i = 1; i < numHist; ++i) {
    if (x(i).size() != numBins) {
      m_isCommonBinsFlag = false;
      break;
    }
  }
 
  // Check that the values of each histogram are identical.
  if (m_isCommonBinsFlag) {
    const size_t lastSpec = numHist - 1;
    for (size_t i = 0; i < lastSpec; ++i) {
      const auto &xi = x(i);
      const auto &xip1 = x(i + 1);
      for (size_t j = 0; j < numBins; ++j) {
        const double a = xi[j];
        const double b = xip1[j];
        // Check for NaN and infinity before comparing for equality
        if (std::isfinite(a) && std::isfinite(b)) {
          if (std::abs(a - b) > EPSILON) {
            m_isCommonBinsFlag = false;
            break;
          }
          // Otherwise we check that both are NaN or both are infinity
        } else if ((std::isnan(a) != std::isnan(b)) || (std::isinf(a) != std::isinf(b))) {
          m_isCommonBinsFlag = false;
          break;
        }
      }
    }
  }
  return m_isCommonBinsFlag;
}
 
bool MatrixWorkspace::isIntegerBins() const {
  if (!this->isCommonBins())
    return false;
 
  const HistogramData::HistogramX bins = x(0);
 
  for (size_t i = 0; i < bins.size(); ++i) {
    if (std::trunc(bins[i]) != bins[i])
      return false;
  }
  return true;
}
 
void MatrixWorkspace::maskBin(const size_t &workspaceIndex, const size_t &binIndex, const double &weight) {
  // First check the workspaceIndex is valid
  if (workspaceIndex >= this->getNumberHistograms())
    throw Kernel::Exception::IndexError(workspaceIndex, this->getNumberHistograms(),
                                        "MatrixWorkspace::maskBin,workspaceIndex");
  // Then check the bin index
  if (binIndex >= y(workspaceIndex).size())
    throw Kernel::Exception::IndexError(binIndex, y(workspaceIndex).size(), "MatrixWorkspace::maskBin,binIndex");
 
  // this function is marked parallel critical
  flagMasked(workspaceIndex, binIndex, weight);
 
  // this is the actual result of the masking that most algorithms and plotting
  // implementations will see, the bin mask flags defined above are used by only
  // some algorithms
  // If the weight is 0, nothing more needs to be done after flagMasked above
  // (i.e. NaN and Inf will also stay intact)
  // If the weight is not 0, NaN and Inf values are set to 0,
  // whereas other values are scaled by (1 - weight)
  if (weight != 0.) {
    double &y = this->mutableY(workspaceIndex)[binIndex];
    (std::isnan(y) || std::isinf(y)) ? y = 0. : y *= (1 - weight);
    double &e = this->mutableE(workspaceIndex)[binIndex];
    (std::isnan(e) || std::isinf(e)) ? e = 0. : e *= (1 - weight);
  }
}
 
void MatrixWorkspace::flagMasked(const size_t &index, const size_t &binIndex, const double &weight) {
  // Writing to m_masks is not thread-safe, so put in some protection
  PARALLEL_CRITICAL(maskBin) {
    // First get a reference to the list for this spectrum (or create a new
    // list)
    MaskList &binList = m_masks[index];
    binList[binIndex] = weight;
  }
}
 
bool MatrixWorkspace::hasMaskedBins(const size_t &workspaceIndex) const {
  // First check the workspaceIndex is valid. Return false if it isn't (decided
  // against throwing here).
  if (workspaceIndex >= this->getNumberHistograms())
    return false;
  return m_masks.find(workspaceIndex) != m_masks.end();
}
 
bool MatrixWorkspace::hasAnyMaskedBins() const { return !m_masks.empty(); }
 
const MatrixWorkspace::MaskList &MatrixWorkspace::maskedBins(const size_t &workspaceIndex) const {
  auto it = m_masks.find(workspaceIndex);
  // Throw if there are no masked bins for this spectrum. The caller should
  // check first using hasMaskedBins!
  if (it == m_masks.end()) {
    throw Kernel::Exception::IndexError(workspaceIndex, 0, "MatrixWorkspace::maskedBins");
  }
 
  return it->second;
}
 
std::vector<size_t> MatrixWorkspace::maskedBinsIndices(const size_t &workspaceIndex) const {
  auto it = m_masks.find(workspaceIndex);
  // Throw if there are no masked bins for this spectrum. The caller should
  // check first using hasMaskedBins!
  if (it == m_masks.end()) {
    throw Kernel::Exception::IndexError(workspaceIndex, 0, "MatrixWorkspace::maskedBins");
  }
 
  auto maskedBins = it->second;
  std::vector<size_t> maskedIds;
  maskedIds.reserve(maskedBins.size());
 
  std::transform(maskedBins.begin(), maskedBins.end(), std::back_inserter(maskedIds),
                 [](const auto &mb) { return mb.first; });
  return maskedIds;
}
 
void MatrixWorkspace::setMaskedBins(const size_t workspaceIndex, const MaskList &maskedBins) {
  m_masks[workspaceIndex] = maskedBins;
}
 
void MatrixWorkspace::setUnmaskedBins(const size_t workspaceIndex) { m_masks.erase(workspaceIndex); }
 
void MatrixWorkspace::setMonitorWorkspace(const std::shared_ptr<MatrixWorkspace> &monitorWS) {
  if (monitorWS.get() == this) {
    throw std::runtime_error("To avoid memory leak, monitor workspace"
                             " can not be the same workspace as the host workspace");
  }
  m_monitorWorkspace = monitorWS;
}
 
std::shared_ptr<MatrixWorkspace> MatrixWorkspace::monitorWorkspace() const { return m_monitorWorkspace; }
 
size_t MatrixWorkspace::getMemorySize() const {
  // 3 doubles per histogram bin.
  return 3 * size() * sizeof(double) + run().getMemorySize();
}
 
size_t MatrixWorkspace::getMemorySizeForXAxes() const {
  size_t total = 0;
  auto lastX = this->refX(0);
  for (size_t wi = 0; wi < getNumberHistograms(); wi++) {
    auto X = this->refX(wi);
    // If the pointers are the same
    if (!(X == lastX) || wi == 0)
      total += (*X).size() * sizeof(double);
  }
  return total;
}
 
Types::Core::DateAndTime MatrixWorkspace::getFirstPulseTime() const {
  TimeSeriesProperty<double> *log = this->run().getTimeSeriesProperty<double>("proton_charge");
 
  DateAndTime startDate = log->firstTime();
  DateAndTime reference("1991-01-01T00:00:00");
 
  int i = 0;
  // Find the first pulse after 1991
  while (startDate < reference && i < 100) {
    i++;
    startDate = log->nthTime(i);
  }
 
  // Return as DateAndTime.
  return startDate;
}
 
Types::Core::DateAndTime MatrixWorkspace::getLastPulseTime() const {
  TimeSeriesProperty<double> *log = this->run().getTimeSeriesProperty<double>("proton_charge");
  return log->lastTime();
}
 
std::size_t MatrixWorkspace::yIndexOfX(const double xValue, const std::size_t &index,
                                       [[maybe_unused]] const double tolerance) const {
  if (index >= getNumberHistograms())
    throw std::out_of_range("MatrixWorkspace::yIndexOfX - Index out of range.");
 
  const auto &xValues = this->x(index);
  const bool ascendingOrder = xValues.front() < xValues.back();
  const auto minX = ascendingOrder ? xValues.front() : xValues.back();
  const auto maxX = ascendingOrder ? xValues.back() : xValues.front();
 
  if (isHistogramDataByIndex(index)) {
    if (xValue < minX || xValue > maxX)
      throw std::out_of_range("MatrixWorkspace::yIndexOfX - X value is out of "
                              "the range of the min and max bin edges.");
 
    return binIndexOfValue(xValues, xValue, ascendingOrder);
  } else {
    if (xValue < minX - tolerance || xValue > maxX + tolerance)
      throw std::out_of_range("MatrixWorkspace::yIndexOfX - X value is out of "
                              "range for this point data.");
 
    return xIndexOfValue(xValues, xValue, tolerance);
  }
}
 
std::size_t MatrixWorkspace::binIndexOfValue(HistogramData::HistogramX const &xValues, const double xValue,
                                             const bool ascendingOrder) const {
  std::size_t hops;
  if (ascendingOrder) {
    auto lowerIter = std::lower_bound(xValues.cbegin(), xValues.cend(), xValue);
 
    // If we are pointing at the first value then we want to be in the first bin
    if (lowerIter == xValues.cbegin())
      ++lowerIter;
 
    hops = std::distance(xValues.cbegin(), lowerIter);
  } else {
    auto lowerIter = std::lower_bound(xValues.crbegin(), xValues.crend(), xValue);
 
    if (lowerIter == xValues.crbegin())
      ++lowerIter;
 
    hops = xValues.size() - std::distance(xValues.crbegin(), lowerIter);
  }
  // The bin index is offset by one from the number of hops between iterators as
  // they start at zero (for a histogram workspace)
  return hops - 1;
}
 
std::size_t MatrixWorkspace::xIndexOfValue(const HistogramData::HistogramX &xValues, const double xValue,
                                           const double tolerance) const {
  auto const iter = std::find_if(xValues.cbegin(), xValues.cend(), [&xValue, &tolerance](double const &value) {
    return std::abs(xValue - value) <= tolerance;
  });
  if (iter != xValues.cend())
    return std::distance(xValues.cbegin(), iter);
  else
    throw std::invalid_argument("MatrixWorkspace::yIndexOfX - the X value provided could not be found "
                                "in the workspace containing point data.");
}
 
uint64_t MatrixWorkspace::getNPoints() const { return static_cast<uint64_t>(this->size()); }
 
//================================= FOR MDGEOMETRY
//====================================================
 
size_t MatrixWorkspace::getNumDims() const { return 2; }
 
std::string MatrixWorkspace::getDimensionIdFromAxis(const int &axisIndex) const {
  std::string id;
  if (0 == axisIndex) {
    id = xDimensionId;
  } else if (1 == axisIndex) {
    id = yDimensionId;
  } else {
    throw std::invalid_argument("Cannot have an index for a MatrixWorkspace "
                                "axis that is not == 0 or == 1");
  }
  return id;
}
 
//===============================================================================
class MWDimension : public Mantid::Geometry::IMDDimension {
public:
  MWDimension(const Axis *axis, std::string dimensionId)
      : m_axis(*axis), m_dimensionId(std::move(dimensionId)),
        m_haveEdges(dynamic_cast<const BinEdgeAxis *>(&m_axis) != nullptr),
        m_frame(std::make_unique<Geometry::GeneralFrame>(m_axis.unit()->label(), m_axis.unit()->label())) {}
 
  std::string getName() const override {
    const auto &unit = m_axis.unit();
    if (unit && unit->unitID() != "Empty")
      return unit->caption();
    else
      return m_axis.title();
  }
 
  const Kernel::UnitLabel getUnits() const override { return m_axis.unit()->label(); }
 
  const std::string &getDimensionId() const override { return m_dimensionId; }
 
  bool getIsIntegrated() const override { return m_axis.length() == 1; }
 
  coord_t getMinimum() const override { return coord_t(m_axis.getMin()); }
 
  coord_t getMaximum() const override { return coord_t(m_axis.getMax()); }
 
  size_t getNBins() const override {
    if (m_haveEdges)
      return m_axis.length() - 1;
    else
      return m_axis.length();
  }
 
  size_t getNBoundaries() const override { return m_axis.length(); }
 
  void setRange(size_t /*nBins*/, coord_t /*min*/, coord_t /*max*/) override {
    throw std::runtime_error("Not implemented");
  }
 
  coord_t getX(size_t ind) const override { return coord_t(m_axis(ind)); }
 
  coord_t getBinWidth() const override {
    size_t nsteps = (m_haveEdges) ? this->getNBins() : this->getNBins() - 1;
    return (getMaximum() - getMinimum()) / static_cast<coord_t>(nsteps);
  }
 
  // Dimensions must be xml serializable.
  std::string toXMLString() const override { throw std::runtime_error("Not implemented"); }
 
  const Kernel::MDUnit &getMDUnits() const override { return m_frame->getMDUnit(); }
  const Geometry::MDFrame &getMDFrame() const override { return *m_frame; }
 
private:
  const Axis &m_axis;
  const std::string m_dimensionId;
  const bool m_haveEdges;
  const Geometry::MDFrame_const_uptr m_frame;
};
 
//===============================================================================
class MWXDimension : public Mantid::Geometry::IMDDimension {
public:
  MWXDimension(const MatrixWorkspace *ws, std::string dimensionId)
      : m_ws(ws), m_X(ws->readX(0)), m_dimensionId(std::move(dimensionId)),
        m_frame(std::make_unique<Geometry::GeneralFrame>(m_ws->getAxis(0)->unit()->label(),
                                                         m_ws->getAxis(0)->unit()->label())) {}
 
  std::string getName() const override {
    const auto *axis = m_ws->getAxis(0);
    const auto &unit = axis->unit();
    if (unit && unit->unitID() != "Empty")
      return unit->caption();
    else
      return axis->title();
  }
 
  const Kernel::UnitLabel getUnits() const override { return m_ws->getAxis(0)->unit()->label(); }
 
  const std::string &getDimensionId() const override { return m_dimensionId; }
 
  bool getIsIntegrated() const override { return m_X.size() == 1; }
 
  coord_t getMinimum() const override { return coord_t(m_X.front()); }
 
  coord_t getMaximum() const override { return coord_t(m_X.back()); }
 
  size_t getNBins() const override { return (m_ws->isHistogramData()) ? m_X.size() - 1 : m_X.size(); }
 
  size_t getNBoundaries() const override { return m_X.size(); }
 
  void setRange(size_t /*nBins*/, coord_t /*min*/, coord_t /*max*/) override {
    throw std::runtime_error("Not implemented");
  }
 
  coord_t getX(size_t ind) const override { return coord_t(m_X[ind]); }
 
  // Dimensions must be xml serializable.
  std::string toXMLString() const override { throw std::runtime_error("Not implemented"); }
  const Kernel::MDUnit &getMDUnits() const override { return m_frame->getMDUnit(); }
  const Geometry::MDFrame &getMDFrame() const override { return *m_frame; }
 
private:
  const MatrixWorkspace *m_ws;
  MantidVec m_X;
  const std::string m_dimensionId;
  const Geometry::MDFrame_const_uptr m_frame;
};
 
std::shared_ptr<const Mantid::Geometry::IMDDimension> MatrixWorkspace::getDimension(size_t index) const {
  if (index == 0) {
    return std::make_shared<MWXDimension>(this, xDimensionId);
  } else if (index == 1) {
    Axis *yAxis = this->getAxis(1);
    return std::make_shared<MWDimension>(yAxis, yDimensionId);
  } else
    throw std::invalid_argument("MatrixWorkspace only has 2 dimensions.");
}
 
std::shared_ptr<const Mantid::Geometry::IMDDimension> MatrixWorkspace::getDimensionWithId(std::string id) const {
  int nAxes = this->axes();
  std::shared_ptr<IMDDimension> dim;
  for (int i = 0; i < nAxes; i++) {
    const std::string knownId = getDimensionIdFromAxis(i);
    if (knownId == id) {
      dim = std::make_shared<MWDimension>(this->getAxis(i), id);
      break;
    }
  }
 
  if (nullptr == dim) {
    std::string message = "Cannot find id : " + id;
    throw std::overflow_error(message);
  }
  return dim;
}
 
std::vector<std::unique_ptr<IMDIterator>>
MatrixWorkspace::createIterators(size_t suggestedNumCores, Mantid::Geometry::MDImplicitFunction *function) const {
  // Find the right number of cores to use
  size_t numCores = suggestedNumCores;
  if (!this->threadSafe())
    numCores = 1;
  size_t numElements = this->getNumberHistograms();
  if (numCores > numElements)
    numCores = numElements;
  if (numCores < 1)
    numCores = 1;
 
  // Create one iterator per core, splitting evenly amongst spectra
  std::vector<std::unique_ptr<IMDIterator>> out;
  for (size_t i = 0; i < numCores; i++) {
    size_t begin = (i * numElements) / numCores;
    size_t end = ((i + 1) * numElements) / numCores;
    if (end > numElements)
      end = numElements;
    out.emplace_back(std::make_unique<MatrixWorkspaceMDIterator>(this, function, begin, end));
  }
  return out;
}
 
IMDWorkspace::LinePlot MatrixWorkspace::getLinePlot(const Mantid::Kernel::VMD &start, const Mantid::Kernel::VMD &end,
                                                    Mantid::API::MDNormalization normalize) const {
  return IMDWorkspace::getLinePlot(start, end, normalize);
}
 
signal_t MatrixWorkspace::getSignalAtCoord(const coord_t *coords,
                                           const Mantid::API::MDNormalization &normalization) const {
  if (this->axes() != 2)
    throw std::invalid_argument("MatrixWorkspace::getSignalAtCoord() - "
                                "Workspace can only have 2 axes, found " +
                                std::to_string(this->axes()));
 
  coord_t xCoord = coords[0];
  coord_t yCoord = coords[1];
  // First, find the workspace index
  Axis *ax1 = this->getAxis(1);
  size_t wi(-1);
  try {
    wi = ax1->indexOfValue(yCoord);
  } catch (std::out_of_range &) {
    return std::numeric_limits<double>::quiet_NaN();
  }
 
  const size_t nhist = this->getNumberHistograms();
  const auto &yVals = this->y(wi);
  double yBinSize(1.0); // only applies for volume normalization & numeric axis
  if (normalization == VolumeNormalization && ax1->isNumeric()) {
    size_t uVI = 0; // unused vertical index.
    double currentVertical = ax1->operator()(wi, uVI);
    if (wi + 1 == nhist && nhist > 1) // On the boundary, look back to get diff
    {
      yBinSize = currentVertical - ax1->operator()(wi - 1, uVI);
    } else {
      yBinSize = ax1->operator()(wi + 1, uVI) - currentVertical;
    }
  }
 
  if (wi < nhist) {
    const auto &xVals = x(wi);
    size_t i;
    try {
      if (isHistogramData())
        i = Kernel::VectorHelper::indexOfValueFromEdges(xVals.rawData(), xCoord);
      else
        i = Kernel::VectorHelper::indexOfValueFromCenters(xVals.rawData(), xCoord);
    } catch (std::out_of_range &) {
      return std::numeric_limits<double>::quiet_NaN();
    }
 
    double y = yVals[i];
    // What is our normalization factor?
    switch (normalization) {
    case NoNormalization:
      return y;
    case VolumeNormalization: {
      // Divide the signal by the area
      auto volume = yBinSize * (xVals[i + 1] - xVals[i]);
      if (volume == 0.0) {
        return std::numeric_limits<double>::quiet_NaN();
      }
      return y / volume;
    }
    case NumEventsNormalization:
      // Not yet implemented, may not make sense
      return y;
    }
    // This won't happen
    return y;
  } else {
    return std::numeric_limits<double>::quiet_NaN();
  }
}
 
signal_t MatrixWorkspace::getSignalWithMaskAtCoord(const coord_t *coords,
                                                   const Mantid::API::MDNormalization &normalization) const {
  return getSignalAtCoord(coords, normalization);
}
 
/*
MDMasking for a Matrix Workspace has not been implemented.
@param :
*/
void MatrixWorkspace::setMDMasking(std::unique_ptr<Mantid::Geometry::MDImplicitFunction> /*maskingRegion*/) {
  throw std::runtime_error("MatrixWorkspace::setMDMasking has no implementation");
}
 
/*
Clear MDMasking for a Matrix Workspace has not been implemented.
*/
void MatrixWorkspace::clearMDMasking() {
  throw std::runtime_error("MatrixWorkspace::clearMDMasking has no implementation");
}
 
Mantid::Kernel::SpecialCoordinateSystem MatrixWorkspace::getSpecialCoordinateSystem() const {
  return Mantid::Kernel::None;
}
 
// Check if this class has an oriented lattice on a sample object
bool MatrixWorkspace::hasOrientedLattice() const { return Mantid::API::ExperimentInfo::sample().hasOrientedLattice(); }
 
MantidImage_sptr MatrixWorkspace::getImage(const MantidVec &(MatrixWorkspace::*read)(std::size_t const) const,
                                           size_t start, size_t stop, size_t width, size_t indexStart,
                                           size_t indexEnd) const {
  // width must be provided (for now)
  if (width == 0) {
    throw std::runtime_error("Cannot create image with width 0");
  }
 
  size_t nHist = getNumberHistograms();
  // use all spectra by default
  if (stop == 0) {
    stop = nHist;
  }
 
  // check start and stop
  if (stop < start) {
    throw std::runtime_error("Cannot create image for an empty data set.");
  }
 
  if (start >= nHist) {
    throw std::runtime_error("Cannot create image: start index is out of range");
  }
 
  if (stop >= nHist) {
    throw std::runtime_error("Cannot create image: stop index is out of range");
  }
 
  // calculate image geometry
  size_t dataSize = stop - start + 1;
  size_t height = dataSize / width;
 
  // and check that the data fits exactly into this geometry
  if (height * width != dataSize) {
    throw std::runtime_error("Cannot create image: the data set cannot form a rectangle.");
  }
 
  size_t nBins = blocksize();
  bool isHisto = isHistogramData();
 
  // default indexEnd is the last index of the X vector
  if (indexEnd == 0) {
    indexEnd = nBins;
    if (!isHisto && indexEnd > 0)
      --indexEnd;
  }
 
  // check the x-range indices
  if (indexEnd < indexStart) {
    throw std::runtime_error("Cannot create image for an empty data set.");
  }
 
  if (indexStart >= nBins || indexEnd > nBins || (!isHisto && indexEnd == nBins)) {
    throw std::runtime_error("Cannot create image: integration interval is out of range.");
  }
 
  // initialize the image
  auto image = std::make_shared<MantidImage>(height);
  if (!isHisto)
    ++indexEnd;
 
  // deal separately with single-binned workspaces: no integration is required
  if (isHisto && indexEnd == indexStart + 1) {
    PARALLEL_FOR_NO_WSP_CHECK()
    for (int i = 0; i < static_cast<int>(height); ++i) {
      auto &row = (*image)[i];
      row.resize(width);
      size_t spec = start + static_cast<size_t>(i) * width;
      for (size_t j = 0; j < width; ++j, ++spec) {
        row[j] = (this->*read)(spec)[indexStart];
      }
    }
  } else {
    // each image pixel is integrated over the x-range [indexStart,indexEnd)
    PARALLEL_FOR_NO_WSP_CHECK()
    for (int i = 0; i < static_cast<int>(height); ++i) {
      auto &row = (*image)[i];
      row.resize(width);
      size_t spec = start + static_cast<size_t>(i) * width;
      for (size_t j = 0; j < width; ++j, ++spec) {
        auto &V = (this->*read)(spec);
        row[j] = std::accumulate(V.begin() + indexStart, V.begin() + indexEnd, 0.0);
      }
    }
  }
 
  return image;
}
 
std::pair<int64_t, int64_t> MatrixWorkspace::findY(double value, const std::pair<int64_t, int64_t> &idx) const {
  std::pair<int64_t, int64_t> out(-1, -1);
  const int64_t numHists = static_cast<int64_t>(this->getNumberHistograms());
  if (std::isnan(value)) {
    for (int64_t i = idx.first; i < numHists; ++i) {
      const auto &Y = this->y(i);
      if (auto it = std::find_if(std::next(Y.begin(), idx.second), Y.end(), [](double v) { return std::isnan(v); });
          it != Y.end()) {
        out = {i, std::distance(Y.begin(), it)};
        break;
      }
    }
  } else {
    for (int64_t i = idx.first; i < numHists; ++i) {
      const auto &Y = this->y(i);
      if (auto it = std::find(std::next(Y.begin(), idx.second), Y.end(), value); it != Y.end()) {
        out = {i, std::distance(Y.begin(), it)};
        break;
      }
    }
  }
  return out;
}
 
std::pair<size_t, size_t> MatrixWorkspace::getImageStartEndXIndices(size_t i, double startX, double endX) const {
  if (startX == EMPTY_DBL())
    startX = x(i).front();
  auto pStart = getXIndex(i, startX, true);
  if (pStart.second != 0.0) {
    throw std::runtime_error("Start X value is required to be on bin boundary.");
  }
  if (endX == EMPTY_DBL())
    endX = x(i).back();
  auto pEnd = getXIndex(i, endX, false, pStart.first);
  if (pEnd.second != 0.0) {
    throw std::runtime_error("End X value is required to be on bin boundary.");
  }
  return std::make_pair(pStart.first, pEnd.first);
}
 
MantidImage_sptr MatrixWorkspace::getImageY(size_t start, size_t stop, size_t width, double startX, double endX) const {
  auto p = getImageStartEndXIndices(0, startX, endX);
  return getImage(&MatrixWorkspace::readY, start, stop, width, p.first, p.second);
}
 
MantidImage_sptr MatrixWorkspace::getImageE(size_t start, size_t stop, size_t width, double startX, double endX) const {
  auto p = getImageStartEndXIndices(0, startX, endX);
  return getImage(&MatrixWorkspace::readE, start, stop, width, p.first, p.second);
}
 
std::pair<size_t, double> MatrixWorkspace::getXIndex(size_t i, double x, bool isLeft, size_t start) const {
  auto &X = this->x(i);
  auto nx = X.size();
 
  // if start out of range - search failed
  if (start >= nx)
    return std::make_pair(nx, 0.0);
  if (start > 0 && start == nx - 1) {
    // starting with the last index is allowed for right boundary search
    if (!isLeft)
      return std::make_pair(start, 0.0);
    return std::make_pair(nx, 0.0);
  }
 
  // consider point data with single value
  if (nx == 1) {
    assert(start == 0);
    if (isLeft)
      return x <= X[start] ? std::make_pair(start, 0.0) : std::make_pair(nx, 0.0);
    return x >= X[start] ? std::make_pair(start, 0.0) : std::make_pair(nx, 0.0);
  }
 
  // left boundaries below start value map to the start value
  if (x <= X[start]) {
    return isLeft ? std::make_pair(start, 0.0) : std::make_pair(nx, 0.0);
  }
  // right boundary search returns last x value for all values above it
  if (x >= X.back()) {
    return !isLeft ? std::make_pair(nx - 1, 0.0) : std::make_pair(nx, 0.0);
  }
 
  // general case: find the boundary index and bin fraction
  auto end = X.end();
  for (auto ix = X.begin() + start + 1; ix != end; ++ix) {
    if (*ix >= x) {
      auto index = static_cast<size_t>(std::distance(X.begin(), ix));
      if (isLeft)
        --index;
      return std::make_pair(index, fabs((X[index] - x) / (*ix - *(ix - 1))));
    }
  }
  // I don't think we can ever get here
  return std::make_pair(nx, 0.0);
}
 
void MatrixWorkspace::setImage(MantidVec &(MatrixWorkspace::*dataVec)(const std::size_t), const MantidImage &image,
                               size_t start, [[maybe_unused]] bool parallelExecution) {
 
  if (image.empty())
    return;
  if (image[0].empty())
    return;
 
  if (blocksize() != 1) {
    throw std::runtime_error("Cannot set image: a single bin workspace is expected.");
  }
 
  size_t height = image.size();
  size_t width = image.front().size();
  size_t dataSize = width * height;
 
  if (start + dataSize > getNumberHistograms()) {
    throw std::runtime_error("Cannot set image: image is bigger than workspace.");
  }
 
  PARALLEL_FOR_IF(parallelExecution)
  for (int i = 0; i < static_cast<int>(height); ++i) {
    auto &row = image[i];
    if (row.size() != width) {
      throw std::runtime_error("Canot set image: image is corrupted.");
    }
    size_t spec = start + static_cast<size_t>(i) * width;
    auto rowEnd = row.end();
    for (auto pixel = row.begin(); pixel != rowEnd; ++pixel, ++spec) {
      (this->*dataVec)(spec)[0] = *pixel;
    }
  }
}
 
void MatrixWorkspace::setImageY(const MantidImage &image, size_t start, bool parallelExecution) {
  setImage(&MatrixWorkspace::dataY, image, start, parallelExecution);
}
 
void MatrixWorkspace::setImageE(const MantidImage &image, size_t start, bool parallelExecution) {
  setImage(&MatrixWorkspace::dataE, image, start, parallelExecution);
}
 
void MatrixWorkspace::invalidateCachedSpectrumNumbers() {
  if (m_isInitialized && storageMode() == Parallel::StorageMode::Distributed && m_indexInfo->communicator().size() > 1)
    throw std::logic_error("Setting spectrum numbers in MatrixWorkspace via "
                           "ISpectrum::setSpectrumNo is not possible in MPI "
                           "runs for distributed workspaces. Use IndexInfo.");
  m_indexInfoNeedsUpdate = true;
}
 
void MatrixWorkspace::updateCachedDetectorGrouping(const size_t index) const {
  setDetectorGrouping(index, getSpectrum(index).getDetectorIDs());
}
 
void MatrixWorkspace::buildDefaultSpectrumDefinitions() {
  const auto &detInfo = detectorInfo();
  size_t numberOfDetectors{detInfo.size()};
  if (numberOfDetectors == 0) {
    // Default to empty spectrum definitions if there is no instrument.
    m_indexInfo->setSpectrumDefinitions(std::vector<SpectrumDefinition>(m_indexInfo->size()));
    return;
  }
  size_t numberOfSpectra = numberOfDetectors * detInfo.scanCount();
  if (numberOfSpectra != m_indexInfo->globalSize())
    throw std::invalid_argument("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.");
  std::vector<SpectrumDefinition> specDefs(m_indexInfo->size());
  if (!detInfo.isScanning() && (numberOfSpectra == m_indexInfo->size())) {
    for (size_t i = 0; i < numberOfSpectra; ++i)
      specDefs[i].add(i);
  } else {
    size_t specIndex = 0;
    size_t globalSpecIndex = 0;
    for (size_t detIndex = 0; detIndex < detInfo.size(); ++detIndex) {
      for (size_t time = 0; time < detInfo.scanCount(); ++time) {
        if (m_indexInfo->isOnThisPartition(Indexing::GlobalSpectrumIndex(globalSpecIndex++)))
          specDefs[specIndex++].add(detIndex, time);
      }
    }
  }
  m_indexInfo->setSpectrumDefinitions(std::move(specDefs));
}
 
void MatrixWorkspace::rebuildDetectorIDGroupings() {
  const auto &detInfo = detectorInfo();
  const auto &allDetIDs = detInfo.detectorIDs();
  const auto &specDefs = m_indexInfo->spectrumDefinitions();
  const auto size = static_cast<int64_t>(m_indexInfo->size());
  enum class ErrorCode { None, InvalidDetIndex, InvalidTimeIndex };
  std::atomic<ErrorCode> errorValue(ErrorCode::None);
#pragma omp parallel for
  for (int64_t i = 0; i < size; ++i) {
    auto &spec = getSpectrum(i);
    // Prevent setting flags that require spectrum definition updates
    spec.setMatrixWorkspace(nullptr, i);
    spec.setSpectrumNo(static_cast<specnum_t>(m_indexInfo->spectrumNumber(i)));
    std::set<detid_t> detIDs;
    for (const auto &index : (*specDefs)[i]) {
      const size_t detIndex = index.first;
      const size_t timeIndex = index.second;
      if (detIndex >= allDetIDs.size()) {
        errorValue = ErrorCode::InvalidDetIndex;
      } else if (timeIndex >= detInfo.scanCount()) {
        errorValue = ErrorCode::InvalidTimeIndex;
      } else {
        detIDs.insert(allDetIDs[detIndex]);
      }
    }
    spec.setDetectorIDs(std::move(detIDs));
  }
  switch (errorValue) {
  case ErrorCode::InvalidDetIndex:
    throw std::invalid_argument("MatrixWorkspace: SpectrumDefinition contains an out-of-range "
                                "detector index, i.e., the spectrum definition does not match "
                                "the instrument in the workspace.");
  case ErrorCode::InvalidTimeIndex:
    throw std::invalid_argument("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.");
  case ErrorCode::None:; // nothing to do
  }
}
 
} // namespace Mantid::API
 
namespace Mantid::Kernel {
 
template <>
MANTID_API_DLL Mantid::API::MatrixWorkspace_sptr
IPropertyManager::getValue<Mantid::API::MatrixWorkspace_sptr>(const std::string &name) const {
  auto *prop = dynamic_cast<PropertyWithValue<Mantid::API::MatrixWorkspace_sptr> *>(getPointerToProperty(name));
  if (prop) {
    return *prop;
  } else {
    std::string message =
        "Attempt to assign property " + name + " to incorrect type. Expected shared_ptr<MatrixWorkspace>.";
    throw std::runtime_error(message);
  }
}
 
template <>
MANTID_API_DLL Mantid::API::MatrixWorkspace_const_sptr
IPropertyManager::getValue<Mantid::API::MatrixWorkspace_const_sptr>(const std::string &name) const {
  auto *prop = dynamic_cast<PropertyWithValue<Mantid::API::MatrixWorkspace_sptr> *>(getPointerToProperty(name));
  if (prop) {
    return prop->operator()();
  } else {
    std::string message =
        "Attempt to assign property " + name + " to incorrect type. Expected const shared_ptr<MatrixWorkspace>.";
    throw std::runtime_error(message);
  }
}
 
} // namespace Mantid::Kernel