SaveNexusProcessed.cpp

// SaveNexusProcessed
// @author Ronald Fowler, based on SaveNexus
#include "MantidAPI/EnabledWhenWorkspaceIsType.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidAPI/IMDHistoWorkspace.h"
#include "MantidAPI/WorkspaceOpOverloads.h"
#include "MantidDataHandling/SaveNexusProcessed.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/OffsetsWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidNexus/NexusFileIO.h"
#include <nexus/NeXusFile.hpp>
#include <boost/regex.hpp>
#include <boost/shared_ptr.hpp>
#include <Poco/File.h>

using namespace Mantid::API;

namespace Mantid {
namespace DataHandling {

using namespace Kernel;
using namespace API;
using namespace DataObjects;
using Geometry::Instrument_const_sptr;

typedef NeXus::NexusFileIO::optional_size_t optional_size_t;

// Register the algorithm into the algorithm factory
DECLARE_ALGORITHM(SaveNexusProcessed)

/// Empty default constructor
SaveNexusProcessed::SaveNexusProcessed()
    : Algorithm(), m_timeProgInit(0.0), prog() {}

//-----------------------------------------------------------------------------------------------
/** Initialisation method.
 *
 */
void SaveNexusProcessed::init() {
  declareProperty(
      new WorkspaceProperty<Workspace>("InputWorkspace", "", Direction::Input),
      "Name of the workspace to be saved");
  // Declare required input parameters for algorithm
  std::vector<std::string> exts;
  exts.push_back(".nxs");
  exts.push_back(".nx5");
  exts.push_back(".xml");

  declareProperty(new FileProperty("Filename", "", FileProperty::Save, exts),
                  "The name of the Nexus file to write, as a full or relative\n"
                  "path");

  // Declare optional parameters (title now optional, was mandatory)
  declareProperty("Title", "", boost::make_shared<NullValidator>(),
                  "A title to describe the saved workspace");
  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
  mustBePositive->setLower(0);

  declareProperty("WorkspaceIndexMin", 0, mustBePositive,
                  "Index number of first spectrum to write, only for single\n"
                  "period data.");
  declareProperty("WorkspaceIndexMax", Mantid::EMPTY_INT(), mustBePositive,
                  "Index of last spectrum to write, only for single period\n"
                  "data.");
  declareProperty(new ArrayProperty<int>("WorkspaceIndexList"),
                  "List of spectrum numbers to read, only for single period\n"
                  "data.");

  declareProperty("Append", false, "Determines whether .nxs file needs to be\n"
                                   "over written or appended");

  declareProperty(
      "PreserveEvents", true,
      "For EventWorkspaces, preserve the events when saving (default).\n"
      "If false, will save the 2D histogram version of the workspace with the "
      "current binning parameters.");
  setPropertySettings(
      "PreserveEvents",
      new EnabledWhenWorkspaceIsType<EventWorkspace>("InputWorkspace", true));

  declareProperty(
      "CompressNexus", false,
      "For EventWorkspaces, compress the Nexus data field (default False).\n"
      "This will make smaller files but takes much longer.");
  setPropertySettings(
      "CompressNexus",
      new EnabledWhenWorkspaceIsType<EventWorkspace>("InputWorkspace", true));
}

/** Get the list of workspace indices to use
 *
 * @param spec :: returns the list of workspace indices
 * @param matrixWorkspace :: pointer to a MatrixWorkspace
 */
void SaveNexusProcessed::getSpectrumList(
    std::vector<int> &spec, MatrixWorkspace_const_sptr matrixWorkspace) {
  std::vector<int> spec_list = getProperty("WorkspaceIndexList");
  int spec_min = getProperty("WorkspaceIndexMin");
  int spec_max = getProperty("WorkspaceIndexMax");
  const bool list = !spec_list.empty();
  const bool interval = (spec_max != Mantid::EMPTY_INT());
  if (spec_max == Mantid::EMPTY_INT())
    spec_max = 0;
  const int numberOfHist =
      static_cast<int>(matrixWorkspace->getNumberHistograms());

  if (interval) {
    if (spec_max < spec_min || spec_max > numberOfHist - 1) {
      g_log.error("Invalid WorkspaceIndex min/max properties");
      throw std::invalid_argument("Inconsistent properties defined");
    }
    spec.reserve(1 + spec_max - spec_min);
    for (int i = spec_min; i <= spec_max; i++)
      spec.push_back(i);
    if (list) {
      for (size_t i = 0; i < spec_list.size(); i++) {
        int s = spec_list[i];
        if (s < 0)
          continue;
        if (s < spec_min || s > spec_max)
          spec.push_back(s);
      }
    }
  } else if (list) {
    spec_max = 0;
    spec_min = numberOfHist - 1;
    for (size_t i = 0; i < spec_list.size(); i++) {
      int s = spec_list[i];
      if (s < 0)
        continue;
      spec.push_back(s);
      if (s > spec_max)
        spec_max = s;
      if (s < spec_min)
        spec_min = s;
    }
  } else {
    spec_min = 0;
    spec_max = numberOfHist - 1;
    spec.reserve(1 + spec_max - spec_min);
    for (int i = spec_min; i <= spec_max; i++)
      spec.push_back(i);
  }
}

void SaveNexusProcessed::doExec(Workspace_sptr inputWorkspace,
                                Mantid::NeXus::NexusFileIO_sptr &nexusFile,
                                const bool keepFile,
                                optional_size_t entryNumber) {
  // TODO: Remove?
  NXMEnableErrorReporting();

  // Retrieve the filename from the properties
  m_filename = getPropertyValue("Filename");
  // m_entryname = getPropertyValue("EntryName");
  m_title = getPropertyValue("Title");
  // Do we prserve events?
  bool PreserveEvents = getProperty("PreserveEvents");

  MatrixWorkspace_const_sptr matrixWorkspace =
      boost::dynamic_pointer_cast<const MatrixWorkspace>(inputWorkspace);
  ITableWorkspace_const_sptr tableWorkspace =
      boost::dynamic_pointer_cast<const ITableWorkspace>(inputWorkspace);
  PeaksWorkspace_const_sptr peaksWorkspace =
      boost::dynamic_pointer_cast<const PeaksWorkspace>(inputWorkspace);
  OffsetsWorkspace_const_sptr offsetsWorkspace =
      boost::dynamic_pointer_cast<const OffsetsWorkspace>(inputWorkspace);
  if (peaksWorkspace)
    g_log.debug("We have a peaks workspace");
  // check if inputWorkspace is something we know how to save
  if (!matrixWorkspace && !tableWorkspace) {
    // get the workspace name for the error message
    std::string name = getProperty("InputWorkspace");

    // md workspaces should be saved using SaveMD
    if (bool(boost::dynamic_pointer_cast<const IMDEventWorkspace>(
            inputWorkspace)) ||
        bool(boost::dynamic_pointer_cast<const IMDHistoWorkspace>(
            inputWorkspace)))
      g_log.warning() << name << " can be saved using SaveMD\n";

    // standard error message
    std::stringstream msg;
    msg << "Workspace \"" << name
        << "\" not saved because it is not of a type we can presently save.";

    throw std::runtime_error(msg.str());
  }
  m_eventWorkspace =
      boost::dynamic_pointer_cast<const EventWorkspace>(matrixWorkspace);
  const std::string workspaceID = inputWorkspace->id();
  if ((workspaceID.find("Workspace2D") == std::string::npos) &&
      (workspaceID.find("RebinnedOutput") == std::string::npos) &&
      !m_eventWorkspace && !tableWorkspace && !offsetsWorkspace)
    throw Exception::NotImplementedError(
        "SaveNexusProcessed passed invalid workspaces. Must be Workspace2D, "
        "EventWorkspace, ITableWorkspace, or OffsetsWorkspace.");

  // Create progress object for initial part - depends on whether events are
  // processed
  if (PreserveEvents && m_eventWorkspace) {
    m_timeProgInit = 0.07; // Events processed 0.05 to 1.0
  } else {
    m_timeProgInit = 1.0; // All work is done in the initial part
  }
  Progress prog_init(this, 0.0, m_timeProgInit, 7);

  // If no title's been given, use the workspace title field
  if (m_title.empty())
    m_title = inputWorkspace->getTitle();

  // get the workspace name to write to file
  std::string wsName = inputWorkspace->getName();

  // If we don't want to append then remove the file if it already exists
  bool append_to_file = getProperty("Append");
  if (!append_to_file && !keepFile) {
    Poco::File file(m_filename);
    if (file.exists())
      file.remove();
  }

  nexusFile->resetProgress(&prog_init);
  nexusFile->openNexusWrite(m_filename, entryNumber);

  // Equivalent C++ API handle
  ::NeXus::File *cppFile = new ::NeXus::File(nexusFile->fileID);

  prog_init.reportIncrement(1, "Opening file");
  if (nexusFile->writeNexusProcessedHeader(m_title, wsName) != 0)
    throw Exception::FileError("Failed to write to file", m_filename);

  prog_init.reportIncrement(1, "Writing header");

  // write instrument data, if present and writer enabled
  if (matrixWorkspace) {
    // Save the instrument names, ParameterMap, sample, run
    matrixWorkspace->saveExperimentInfoNexus(cppFile);
    prog_init.reportIncrement(1, "Writing sample and instrument");

    // check if all X() are in fact the same array
    const bool uniformSpectra =
        API::WorkspaceHelpers::commonBoundaries(matrixWorkspace);

    // Retrieve the workspace indices (from params)
    std::vector<int> spec;
    this->getSpectrumList(spec, matrixWorkspace);

    prog_init.reportIncrement(1, "Writing data");
    // Write out the data (2D or event)
    if (m_eventWorkspace && PreserveEvents) {
      this->execEvent(nexusFile.get(), uniformSpectra, spec);
    } else if (offsetsWorkspace) {
      g_log.warning() << "Writing SpecialWorkspace2D ID=" << workspaceID
                      << "\n";
      nexusFile->writeNexusProcessedData2D(matrixWorkspace, uniformSpectra,
                                           spec, "offsets_workspace", true);
    } else {
      nexusFile->writeNexusProcessedData2D(matrixWorkspace, uniformSpectra,
                                           spec, "workspace", true);
    }

    cppFile->openGroup("instrument", "NXinstrument");
    matrixWorkspace->saveSpectraMapNexus(cppFile, spec, ::NeXus::LZW);
    cppFile->closeGroup();

  } // finish matrix workspace specifics

  if (peaksWorkspace) {
    // Save the instrument names, ParameterMap, sample, run
    peaksWorkspace->saveExperimentInfoNexus(cppFile);
    prog_init.reportIncrement(1, "Writing sample and instrument");
  }

  // peaks workspace specifics
  if (peaksWorkspace) {
    //  g_log.information("Peaks Workspace saving to Nexus would be done");
    //  int pNum = peaksWorkspace->getNumberPeaks();
    peaksWorkspace->saveNexus(cppFile);

  }                        // finish peaks workspace specifics
  else if (tableWorkspace) // Table workspace specifics
  {
    nexusFile->writeNexusTableWorkspace(tableWorkspace, "table_workspace");
  } // finish table workspace specifics

  // Switch to the Cpp API for the algorithm history
  if (trackingHistory()) {
    m_history->fillAlgorithmHistory(
        this, Mantid::Kernel::DateAndTime::getCurrentTime(), 0,
        Algorithm::g_execCount);
    if (!isChild()) {
      inputWorkspace->history().addHistory(m_history);
    }
    // this is a child algorithm, but we still want to keep the history.
    else if (isRecordingHistoryForChild() && m_parentHistory) {
      m_parentHistory->addChildHistory(m_history);
    }
  }

  inputWorkspace->history().saveNexus(cppFile);
  nexusFile->closeGroup();
  return;
}

//-----------------------------------------------------------------------------------------------
/** Executes the algorithm for a single workspace.
 *
 *  @throw runtime_error Thrown if algorithm cannot execute
 */
void SaveNexusProcessed::exec() {
  Workspace_sptr inputWorkspace = getProperty("InputWorkspace");

  // Then immediately open the file
  auto nexusFile = boost::make_shared<Mantid::NeXus::NexusFileIO>();

  // Perform the execution.
  doExec(inputWorkspace, nexusFile);
}

//-------------------------------------------------------------------------------------
/** Append out each field of a vector of events to separate array.
 *
 * @param events :: vector of TofEvent or WeightedEvent, etc.
 * @param offset :: where the first event goes in the array
 * @param tofs, weights, errorSquareds, pulsetimes :: arrays to write to.
 *        Must be initialized and big enough,
 *        or NULL if they are not meant to be written to.
 */
template <class T>
void SaveNexusProcessed::appendEventListData(std::vector<T> events,
                                             size_t offset, double *tofs,
                                             float *weights,
                                             float *errorSquareds,
                                             int64_t *pulsetimes) {
  // Do nothing if there are no events.
  if (events.empty())
    return;

  typename std::vector<T>::const_iterator it;
  typename std::vector<T>::const_iterator it_end = events.end();
  size_t i = offset;

  // Fill the C-arrays with the fields from all the events, as requested.
  for (it = events.begin(); it != it_end; it++) {
    if (tofs)
      tofs[i] = it->tof();
    if (weights)
      weights[i] = static_cast<float>(it->weight());
    if (errorSquareds)
      errorSquareds[i] = static_cast<float>(it->errorSquared());
    if (pulsetimes)
      pulsetimes[i] = it->pulseTime().totalNanoseconds();
    i++;
  }
}

//-----------------------------------------------------------------------------------------------
/** Execute the saving of event data.
 * This will make one long event list for all events contained.
 * */
void SaveNexusProcessed::execEvent(Mantid::NeXus::NexusFileIO *nexusFile,
                                   const bool uniformSpectra,
                                   const std::vector<int> spec) {
  prog = new Progress(this, m_timeProgInit, 1.0,
                      m_eventWorkspace->getNumberEvents() * 2);

  // Start by writing out the axes and crap
  nexusFile->writeNexusProcessedData2D(m_eventWorkspace, uniformSpectra, spec,
                                       "event_workspace", false);

  // Make a super long list of tofs, weights, etc.
  std::vector<int64_t> indices;
  indices.reserve(m_eventWorkspace->getNumberHistograms() + 1);
  // First we need to index the events in each spectrum
  size_t index = 0;
  for (int wi = 0;
       wi < static_cast<int>(m_eventWorkspace->getNumberHistograms()); wi++) {
    indices.push_back(index);
    // Track the total # of events
    index += m_eventWorkspace->getEventList(wi).getNumberEvents();
  }
  indices.push_back(index);

  // Initialize all the arrays
  int64_t num = index;
  double *tofs = NULL;
  float *weights = NULL;
  float *errorSquareds = NULL;
  int64_t *pulsetimes = NULL;

  // overall event type.
  EventType type = m_eventWorkspace->getEventType();
  bool writeTOF = true;
  bool writePulsetime = false;
  bool writeWeight = false;
  bool writeError = false;

  switch (type) {
  case TOF:
    writePulsetime = true;
    break;
  case WEIGHTED:
    writePulsetime = true;
    writeWeight = true;
    writeError = true;
    break;
  case WEIGHTED_NOTIME:
    writeWeight = true;
    writeError = true;
    break;
  }

  // --- Initialize the combined event arrays ----
  if (writeTOF)
    tofs = new double[num];
  if (writeWeight)
    weights = new float[num];
  if (writeError)
    errorSquareds = new float[num];
  if (writePulsetime)
    pulsetimes = new int64_t[num];

  // --- Fill in the combined event arrays ----
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int wi = 0;
       wi < static_cast<int>(m_eventWorkspace->getNumberHistograms()); wi++) {
    PARALLEL_START_INTERUPT_REGION
    const DataObjects::EventList &el = m_eventWorkspace->getEventList(wi);

    // This is where it will land in the output array.
    // It is okay to write in parallel since none should step on each other.
    size_t offset = indices[wi];

    switch (el.getEventType()) {
    case TOF:
      appendEventListData(el.getEvents(), offset, tofs, weights, errorSquareds,
                          pulsetimes);
      break;
    case WEIGHTED:
      appendEventListData(el.getWeightedEvents(), offset, tofs, weights,
                          errorSquareds, pulsetimes);
      break;
    case WEIGHTED_NOTIME:
      appendEventListData(el.getWeightedEventsNoTime(), offset, tofs, weights,
                          errorSquareds, pulsetimes);
      break;
    }
    prog->reportIncrement(el.getNumberEvents(), "Copying EventList");

    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  /*Default = DONT compress - much faster*/
  bool CompressNexus = getProperty("CompressNexus");

  // Write out to the NXS file.
  nexusFile->writeNexusProcessedDataEventCombined(m_eventWorkspace, indices,
                                                  tofs, weights, errorSquareds,
                                                  pulsetimes, CompressNexus);

  // Free mem.
  delete[] tofs;
  delete[] weights;
  delete[] errorSquareds;
  delete[] pulsetimes;
}

//-----------------------------------------------------------------------------------------------
/** virtual method to set the non workspace properties for this algorithm
 *  @param alg :: pointer to the algorithm
 *  @param propertyName :: name of the property
 *  @param propertyValue :: value  of the property
 *  @param perioidNum :: period number
 */
void SaveNexusProcessed::setOtherProperties(IAlgorithm *alg,
                                            const std::string &propertyName,
                                            const std::string &propertyValue,
                                            int perioidNum) {
  if (!propertyName.compare("Append")) {
    if (perioidNum != 1) {
      alg->setPropertyValue(propertyName, "1");
    } else
      alg->setPropertyValue(propertyName, propertyValue);
  } else
    Algorithm::setOtherProperties(alg, propertyName, propertyValue, perioidNum);
}

/**
 Overriden process groups.
 */
bool SaveNexusProcessed::processGroups() {
  // Then immediately open the file
  auto nexusFile = boost::make_shared<Mantid::NeXus::NexusFileIO>();

  /* Unless we have explicity been asked to append to the file. We should assume
  that we can remove any existing
  files of the same name prior to processing. */
  bool append_to_file = this->getProperty("Append");
  if (!append_to_file) {
    const std::string filename = getPropertyValue("Filename");
    Poco::File file(filename);
    if (file.exists()) {
      file.remove();
    }
  }

  // Only the input workspace property can take group workspaces. Therefore
  // index = 0.
  std::vector<Workspace_sptr> &thisGroup = m_groups[0];
  if (!thisGroup.empty()) {
    for (size_t entry = 0; entry < m_groupSize; entry++) {
      Workspace_sptr ws = thisGroup[entry];
      this->doExec(ws, nexusFile, true /*keepFile*/, entry);
      std::stringstream buffer;
      buffer << "Saving group index " << entry;
      m_log.information(buffer.str());
    }
  }

  nexusFile->closeNexusFile();
  // We finished successfully.
  setExecuted(true);
  notificationCenter().postNotification(
      new FinishedNotification(this, isExecuted()));

  return true;
}

} // namespace DataHandling
} // namespace Mantid