MatrixWorkspace.cpp

#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/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/Strings.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/make_unique.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidParallel/Communicator.h"
#include "MantidTypes/SpectrumDefinition.h"

#include <cmath>

#include <functional>
#include <numeric>

using Mantid::Types::Core::DateAndTime;
using Mantid::Kernel::TimeSeriesProperty;

namespace Mantid {
namespace API {
using std::size_t;
using namespace Geometry;
using Kernel::V3D;

namespace {
/// static logger
Kernel::Logger g_log("MatrixWorkspace");
} // namespace
const std::string MatrixWorkspace::xDimensionId = "xDimension";
const std::string MatrixWorkspace::yDimensionId = "yDimension";

/// Default constructor
MatrixWorkspace::MatrixWorkspace(const Parallel::StorageMode storageMode)
    : IMDWorkspace(storageMode), ExperimentInfo(), m_axes(),
      m_isInitialized(false), m_YUnit(), m_YUnitLabel(),
      m_isCommonBinsFlagSet(false), m_isCommonBinsFlag(false), m_masks() {}

MatrixWorkspace::MatrixWorkspace(const MatrixWorkspace &other)
    : IMDWorkspace(other), ExperimentInfo(other) {
  m_indexInfo = Kernel::make_unique<Indexing::IndexInfo>(other.indexInfo());
  m_indexInfoNeedsUpdate = false;
  m_axes.resize(other.m_axes.size());
  for (size_t i = 0; i < m_axes.size(); ++i)
    m_axes[i] = other.m_axes[i]->clone(this);

  m_isInitialized = other.m_isInitialized;
  m_YUnit = other.m_YUnit;
  m_YUnitLabel = other.m_YUnitLabel;
  m_isCommonBinsFlagSet = other.m_isCommonBinsFlagSet;
  m_isCommonBinsFlag = other.m_isCommonBinsFlag;
  m_masks = other.m_masks;
  // TODO: Do we need to init m_monitorWorkspace?
}

/// Destructor
// 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() {
  for (auto &axis : m_axes) {
    delete axis;
  }
}

/** Returns a const reference to the IndexInfo object of the workspace.
 *
 * Used for access to spectrum number and detector ID information of spectra.
 * Writing spectrum number or detector ID groupings of any spectrum in the
 * workspace will invalidate this reference. */
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.push_back(getSpectrum(i).getSpectrumNo());
    m_indexInfo->setSpectrumNumbers(std::move(spectrumNumbers));
    m_indexInfoNeedsUpdate = false;
  }
  return *m_indexInfo;
}

/** Sets the IndexInfo object of the workspace.
 *
 * Used for setting spectrum number and detector ID information of spectra */
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 = Kernel::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());
}

/// Variant of setIndexInfo, used by WorkspaceFactoryImpl.
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());
}

/// @returns A human-readable string of the current state
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();
}

/** Initialize the workspace. Calls the protected init() method, which is
 * implemented in each type of
 *  workspace. Returns immediately if the workspace is already initialized.
 *  @param NVectors :: The number of spectra in the workspace (only relevant for
 * a 2D workspace
 *  @param XLength :: The number of X data points/bin boundaries in each vector
 * (must all be the same)
 *  @param YLength :: The number of data/error points in each vector (must all
 * be the same)
 */
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 = Kernel::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;
}

//---------------------------------------------------------------------------------------
/** Set the title of the workspace
 *
 *  @param t :: The title
 */
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);
}

/** Get the workspace title
 *
 *  @return The title
 */
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";
    }
  }
}

/**
 * Rebuild the default spectra mapping for a workspace. If a non-empty
 * instrument is set then the default maps each detector to a spectra with
 * the same ID. If an empty instrument is set then a 1:1 map from 1->NHistograms
 * is created.
 * @param includeMonitors :: If false the monitors are not included
 */
void MatrixWorkspace::rebuildSpectraMapping(const bool includeMonitors) {
  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;
      // By default: Spectrum number = index +  1
      const specnum_t specNo = specnum_t(index + 1);

      if (index < this->getNumberHistograms()) {
        auto &spec = getSpectrum(index);
        spec.setSpectrumNo(specNo);
        spec.setDetectorID(detId);
      }

      index++;
    }

  } catch (std::runtime_error &) {
    throw;
  }
}

/** Return a map where:
 *    KEY is the Spectrum #
 *    VALUE is the Workspace Index
 */
spec2index_map MatrixWorkspace::getSpectrumToWorkspaceIndexMap() const {
  SpectraAxis *ax = dynamic_cast<SpectraAxis *>(this->m_axes[1]);
  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;
  }
}

/** Return a vector where:
 *    The index into the vector = spectrum number + offset
 *    The value at that index = the corresponding Workspace Index
 *
 *  @returns :: vector set to above definition
 *  @param offset :: add this to the detector ID to get the index into the
 *vector.
 */
std::vector<size_t>
MatrixWorkspace::getSpectrumToWorkspaceIndexVector(specnum_t &offset) const {
  SpectraAxis *ax = dynamic_cast<SpectraAxis *>(this->m_axes[1]);
  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;
}

/** Does the workspace has any grouped detectors?
 *  @return true if the workspace has any grouped detectors, otherwise false
 */
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;
}

/** Return a map where:
 *    KEY is the DetectorID (pixel ID)
 *    VALUE is the Workspace Index
 *  @param throwIfMultipleDets :: set to true to make the algorithm throw an
 * error
 *         if there is more than one detector for a specific workspace index.
 *  @throw runtime_error if there is more than one detector per spectrum (if
 * throwIfMultipleDets is true)
 *  @return Index to Index Map object. THE CALLER TAKES OWNERSHIP OF THE MAP AND
 * IS RESPONSIBLE FOR ITS DELETION.
 */
detid2index_map MatrixWorkspace::getDetectorIDToWorkspaceIndexMap(
    bool throwIfMultipleDets) const {
  detid2index_map map;

  // Loop through the workspace index
  for (size_t workspaceIndex = 0; workspaceIndex < this->getNumberHistograms();
       ++workspaceIndex) {
    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;
}

/** Return a vector where:
 *    The index into the vector = DetectorID (pixel ID) + offset
 *    The value at that index = the corresponding Workspace Index
 *
 *  @param offset :: add this to the detector ID to get the index into the
 *vector.
 *  @param throwIfMultipleDets :: set to true to make the algorithm throw an
 *error if there is more than one detector for a specific workspace index.
 *  @throw runtime_error if there is more than one detector per spectrum (if
 *throwIfMultipleDets is true)
 *  @returns :: vector set to above definition
 */
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;
}

/** Converts a list of spectrum numbers to the corresponding workspace indices.
 *  Not a very efficient operation, but unfortunately it's sometimes required.
 *
 *  @param spectraList :: The list of spectrum numbers required
 *  @returns :: the vector of indices (empty if not a Workspace2D)
 */
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.push_back(i);
        break;
      }
    }
    ++iter;
  }
  return indexList;
}

/** Given a spectrum number, find the corresponding workspace index
 *
 * @param specNo :: spectrum number wanted
 * @return the workspace index
 * @throw runtime_error if not found.
 */
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.");
}

/** Converts a list of detector IDs to the corresponding workspace indices.
 *
 *  Note that only known detector IDs are converted (so an empty vector will
 *be returned
 *  if none of the IDs are recognised), and that the returned workspace
 *indices are
 *  effectively a set (i.e. there are no duplicates).
 *
 *  @param detIdList :: The list of detector IDs required
 *  @returns :: a vector of indices
 */
std::vector<size_t> MatrixWorkspace::getIndicesFromDetectorIDs(
    const std::vector<detid_t> &detIdList) const {
  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()) {
      for (auto index : wsIndices->second) {
        indexList.push_back(index);
      }
    }
  }
  return indexList;
}

/** Converts a list of detector IDs to the corresponding spectrum numbers. Might
 *be slow!
 *
 * @param detIdList :: The list of detector IDs required
 * @returns :: a reference to the vector of spectrum numbers.
 *                       0 for not-found detectors
 */
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.push_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 {
  if (m_indexInfo->size() != m_indexInfo->globalSize())
    throw std::runtime_error(
        "MatrixWorkspace: Parallel support for XMin and XMax not implemented.");

  // 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;
    }
  }
}

/** Integrate all the spectra in the matrix workspace within the range given.
 * Default implementation, can be overridden by base classes if they know
 *something smarter!
 *
 * @param out :: returns the vector where there is one entry per spectrum in the
 *workspace. Same
 *            order as the workspace indices.
 * @param minX :: minimum X bin to use in integrating.
 * @param maxX :: maximum X bin to use in integrating.
 * @param entireRange :: set to true to use the entire range. minX and maxX are
 *then ignored!
 */
void MatrixWorkspace::getIntegratedSpectra(std::vector<double> &out,
                                           const double minX, const double maxX,
                                           const bool entireRange) const {
  out.resize(this->getNumberHistograms(), 0.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() <= 2) && (!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() - 1;

      // 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 lowit
        if ((*highit) > 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);
      }
      // Save it in the vector
      out[wksp_index] = sum;
    }
  }
}

/** Get the effective detector for the given spectrum
*  @param  workspaceIndex The workspace index for which the detector is required
*  @return A single detector object representing the detector(s) contributing
*          to the given spectrum number. If more than one detector contributes
then
*          the returned object's concrete type will be DetectorGroup.
*  @throw  Kernel::Exception::NotFoundError If the Instrument is missing or the
requested workspace index does not have any associated detectors
*/
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 boost::make_shared<Geometry::DetectorGroup>(dets_ptr);
}

/** Returns the signed 2Theta scattering angle for a detector
 *  @param det :: A pointer to the detector object (N.B. might be a
 * DetectorGroup)
 *  @return The scattering angle (0 < theta < pi)
 *  @throws InstrumentDefinitionError if source or sample is missing, or they
 * are in the same place
 */
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 &instrumentUpAxis =
      instrument->getReferenceFrame()->vecPointingUp();
  return det.getSignedTwoTheta(samplePos, beamLine, instrumentUpAxis);
}

/** Returns the 2Theta scattering angle for a detector
 *  @param det :: A pointer to the detector object (N.B. might be a
 * DetectorGroup)
 *  @return The scattering angle (0 < theta < pi)
 *  @throws InstrumentDefinitionError if source or sample is missing, or they
 * are in the same place
 */
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);
}

/// @return The number of axes which this workspace has
int MatrixWorkspace::axes() const { return static_cast<int>(m_axes.size()); }

/** Get a pointer to a workspace axis
 *  @param axisIndex :: The index of the axis required
 *  @throw IndexError If the argument given is outside the range of axes held by
 * this workspace
 *  @return Pointer to Axis object
 */
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];
}

/** Replaces one of the workspace's axes with the new one provided.
 *  @param axisIndex :: The index of the axis to replace
 *  @param newAxis :: A pointer to the new axis. The class will take ownership.
 *  @throw IndexError If the axisIndex given is outside the range of axes held
 * by this workspace
 *  @throw std::runtime_error If the new axis is not of the correct length
 * (within one of the old one)
 */
void MatrixWorkspace::replaceAxis(const std::size_t &axisIndex,
                                  Axis *const 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
  delete m_axes[axisIndex];
  m_axes[axisIndex] = newAxis;
}

/// Returns the units of the data in the workspace
std::string MatrixWorkspace::YUnit() const { return m_YUnit; }

/// Sets a new unit for the data (Y axis) in the workspace
void MatrixWorkspace::setYUnit(const std::string &newUnit) {
  m_YUnit = newUnit;
}

/// Returns a caption for the units of the data in the workspace
std::string MatrixWorkspace::YUnitLabel() const {
  std::string retVal;
  if (!m_YUnitLabel.empty())
    retVal = m_YUnitLabel;
  else {
    retVal = m_YUnit;
    // If this workspace a distribution & has at least one axis & this axis has
    // its unit set
    // then append that unit to the string to be returned
    if (!retVal.empty() && this->isDistribution() && this->axes() &&
        this->getAxis(0)->unit()) {
      retVal = retVal + " per " + this->getAxis(0)->unit()->label().ascii();
    }
  }

  return retVal;
}

/// Sets a new caption for the data (Y axis) in the workspace
void MatrixWorkspace::setYUnitLabel(const std::string &newLabel) {
  m_YUnitLabel = newLabel;
}

/** Are the Y-values in this workspace dimensioned?
 * TODO: For example: ????
 * @return whether workspace is a distribution or not
 */
bool MatrixWorkspace::isDistribution() const {
  return getSpectrum(0).yMode() == HistogramData::Histogram::YMode::Frequencies;
}

/** Set the flag for whether the Y-values are dimensioned
 *  @return whether workspace is now a distribution
 */
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);
}

/**
 *  Whether the workspace contains histogram data
 *  @return whether the workspace contains histogram data
 */
bool MatrixWorkspace::isHistogramData() const {
  // all spectra *should* have the same behavior
  bool isHist = (x(0).size() != y(0).size());
  // TODOHIST temporary sanity check
  if (isHist) {
    if (getSpectrum(0).histogram().xMode() !=
        HistogramData::Histogram::XMode::BinEdges) {
      throw std::logic_error("In MatrixWorkspace::isHistogramData(): "
                             "Histogram::Xmode is not BinEdges");
    }
  } else {
    if (getSpectrum(0).histogram().xMode() !=
        HistogramData::Histogram::XMode::Points) {
      throw std::logic_error("In MatrixWorkspace::isHistogramData(): "
                             "Histogram::Xmode is not Points");
    }
  }
  return isHist;
}

/**
 *  Whether the workspace contains common X bins
 *  @return whether the workspace contains common X bins
 */
bool MatrixWorkspace::isCommonBins() const {
  if (!m_isCommonBinsFlagSet) {
    m_isCommonBinsFlag = true;

    const size_t numHist = getNumberHistograms();
    // there being only one or zero histograms is accepted as not being an error
    if (numHist > 1) {
      const size_t numBins = x(0).size();
      for (size_t i = 1; i < numHist; ++i) {
        if (x(i).size() != numBins) {
          m_isCommonBinsFlag = false;
          break;
        }
      }

      // there being only one or zero histograms is accepted as not being an
      // error
      if (m_isCommonBinsFlag) {
        // otherwise will compare some of the data, to save time just check two
        // the first and the last
        const size_t lastSpec = numHist - 1;
        // Quickest check is to see if they are actually the same vector
        if (&(x(0)[0]) != &(x(lastSpec)[0])) {
          // Now check numerically
          const double first = std::accumulate(x(0).begin(), x(0).end(), 0.);
          const double last =
              std::accumulate(x(lastSpec).begin(), x(lastSpec).end(), 0.);
          if (std::abs(first - last) / std::abs(first + last) > 1.0E-9) {
            m_isCommonBinsFlag = false;
          }

          // handle Nan's and inf's
          if ((std::isinf(first) != std::isinf(last)) ||
              (std::isnan(first) != std::isnan(last))) {
            m_isCommonBinsFlag = false;
          }
        }
      }
    }
    m_isCommonBinsFlagSet = true;
  }

  return m_isCommonBinsFlag;
}

/** Called by the algorithm MaskBins to mask a single bin for the first time,
 * algorithms that later propagate the
 *  the mask from an input to the output should call flagMasked() instead. Here
 * y-values and errors will be scaled
 *  by (1-weight) as well as the mask flags (m_masks) being updated. This
 * function doesn't protect the writes to the
 *  y and e-value arrays and so is not safe if called by multiple threads
 * working on the same spectrum. Writing to the mask set is marked parrallel