ExportTimeSeriesLog.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 "MantidAlgorithms/ExportTimeSeriesLog.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/IEventList.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/Events.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/Instrument.h"
#include "MantidHistogramData/Histogram.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/System.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/UnitFactory.h"
#include <algorithm>
#include <fstream>
#include <sstream>
 
using namespace Mantid;
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::HistogramData;
using Mantid::Types::Core::DateAndTime;
 
using namespace std;
 
namespace Mantid::Algorithms {
 
DECLARE_ALGORITHM(ExportTimeSeriesLog)
 
 
void ExportTimeSeriesLog::init() {
  declareProperty(
      std::make_unique<API::WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "Anonymous", Direction::InOut),
      "Name of input Matrix workspace containing the log to export. ");
 
  declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "Dummy", Direction::Output),
                  "Name of the workspace containing the log events in Export. ");
 
  declareProperty("CalculateFirstDerivative", false,
                  "If specified then the first derivative of exported data "
                  "will be calcualted and put to spectrum 1.");
 
  declareProperty("LogName", "", "Log's name to filter events.");
 
  std::vector<std::string> units{"Seconds", "Nano Seconds"};
  declareProperty("UnitOfTime", "Seconds", std::make_shared<Kernel::StringListValidator>(units),
                  "StartTime, StopTime and DeltaTime can be given in various unit."
                  "The unit can be 'Seconds' or 'Nanoseconds' from run start time."
                  "They can also be defined as 'Percentage' of total run time.");
 
  declareProperty("StartTime", EMPTY_DBL(),
                  "Relative starting time of the output series. "
                  "Its unit is determined by property UnitOfTime.");
 
  declareProperty("StopTime", EMPTY_DBL(),
                  "Relative stopping time of the output series."
                  "Its unit is determined by property UnitOfTime.");
 
  declareProperty("OutputAbsoluteTime", false, "If true, the output times will be absolute time to 1990.01.01.");
 
  declareProperty("NumberEntriesExport", EMPTY_INT(),
                  "Number of entries of the log to be exported.  Default is all entries.");
 
  declareProperty("IsEventWorkspace", true,
                  "If set to true, output workspace "
                  "is EventWorkspace.  Otherwise, it "
                  "is Workspace2D.");
}
 
void ExportTimeSeriesLog::exec() {
  // Get properties
  m_inputWS = this->getProperty("InputWorkspace");
  string logname = getProperty("LogName");
 
  double start_time = getProperty("StartTime");
  double stop_time = getProperty("StopTime");
  std::string time_unit = getProperty("UnitOfTime");
  bool exportEpochTime = getProperty("OutputAbsoluteTime");
 
  int numberoutputentries = getProperty("NumberEntriesExport");
  bool outputeventworkspace = getProperty("IsEventWorkspace");
 
  bool cal1stderiv = getProperty("CalculateFirstDerivative");
 
  // Call the main
  exportLog(logname, time_unit, start_time, stop_time, exportEpochTime, outputeventworkspace, numberoutputentries,
            cal1stderiv);
 
  // calcualte first derivative
  if (cal1stderiv)
    calculateFirstDerivative(outputeventworkspace);
 
  // set up the sample log values for meta information
  setupMetaData(logname, time_unit, exportEpochTime);
 
  // 3. Output
  setProperty("OutputWorkspace", m_outWS);
}
 
void ExportTimeSeriesLog::exportLog(const std::string &logname, const std::string &timeunit, const double &starttime,
                                    const double &stoptime, const bool exportepoch, bool outputeventws, int numentries,
                                    bool cal_first_deriv) {
 
  // Get log, time, and etc.
  std::vector<Types::Core::DateAndTime> times;
  std::vector<double> values;
 
  if (!logname.empty()) {
    // Log
    auto *tlog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(m_inputWS->run().getProperty(logname));
    if (!tlog) {
      std::stringstream errmsg;
      errmsg << "TimeSeriesProperty Log " << logname << " does not exist in workspace " << m_inputWS->getName();
      g_log.error(errmsg.str());
      throw std::invalid_argument(errmsg.str());
    }
    times = tlog->timesAsVector();
    values = tlog->valuesAsVector();
  } else {
    throw std::runtime_error("Log name cannot be left empty.");
  }
 
  // Get start time, stop time and unit factor
  double timeunitfactor = 1.;
  if (timeunit == "Seconds")
    timeunitfactor = 1.E-9;
 
  // Get index for start time
  size_t i_start = 0;
  size_t i_stop = times.size() - 1;
  // Rule out the case that start time is behind last log entry
  bool i_start_cal = false;
  if (starttime != EMPTY_DBL()) {
    int64_t timerangens = times.back().totalNanoseconds() - times.front().totalNanoseconds();
    double timerange = static_cast<double>(timerangens) * timeunitfactor;
    g_log.debug() << "Time range is " << timerange << ", Start time is " << starttime << "\n";
    if (timerange < starttime) {
      i_start = times.size() - 1;
      i_start_cal = true;
    }
  }
 
  if ((!i_start_cal) && (starttime != EMPTY_DBL() || stoptime != EMPTY_DBL())) {
    bool export_partial = calculateTimeSeriesRangeByTime(times, starttime, i_start, stoptime, i_stop, timeunitfactor);
    if (!export_partial)
      throw std::runtime_error("Unable to find proton_charge for run start time. "
                               "Failed to get partial time series.");
  }
 
  // Determine number of export log
  if (numentries == EMPTY_INT()) {
    numentries = static_cast<int>(times.size());
  } else if (numentries <= 0) {
    stringstream errmsg;
    errmsg << "For Export Log, NumberEntriesExport must be greater than 0.  "
              "Input = "
           << numentries;
    g_log.error(errmsg.str());
    throw std::runtime_error(errmsg.str());
  } else if (static_cast<size_t>(numentries) > times.size()) {
    numentries = static_cast<int>(times.size());
  }
 
  // Create otuput workspace
  if (outputeventws) {
    setupEventWorkspace(i_start, i_stop, numentries, times, values, exportepoch);
  } else {
    size_t nspec(1);
    if (cal_first_deriv)
      nspec = 2;
    setupWorkspace2D(i_start, i_stop, numentries, times, values, exportepoch, timeunitfactor, nspec);
  }
}
 
void ExportTimeSeriesLog::setupWorkspace2D(const size_t &start_index, const size_t &stop_index, int numentries,
                                           vector<DateAndTime> &times, vector<double> values, const bool &epochtime,
                                           const double &timeunitfactor, size_t nspec) {
  // Determine time shift
  int64_t timeshift(0);
  if (!epochtime) {
    // relative time
    Types::Core::DateAndTime runstart(m_inputWS->run().getProperty("run_start")->value());
    timeshift = runstart.totalNanoseconds();
  }
 
  // Determine the size
  size_t outsize = stop_index - start_index + 1;
  if (outsize > static_cast<size_t>(numentries))
    outsize = static_cast<size_t>(numentries);
 
  // Create 2D workspace
  m_outWS = create<Workspace2D>(nspec, Points(outsize));
 
  auto &vecX = m_outWS->mutableX(0);
  auto &vecY = m_outWS->mutableY(0);
  auto &vecE = m_outWS->mutableE(0);
 
  for (size_t index = 0; index < outsize; ++index) {
    size_t i_time = index + start_index;
    // safety check
    if (i_time >= times.size()) {
      std::stringstream errss;
      errss << "It shouldn't happen that the index is out of boundary."
            << "start index = " << start_index << ", output size = " << outsize << ", index = " << index << "\n";
      throw std::runtime_error(errss.str());
    }
 
    int64_t dtns = times[i_time].totalNanoseconds() - timeshift;
    vecX[index] = static_cast<double>(dtns) * timeunitfactor;
    vecY[index] = values[i_time];
    vecE[index] = 0.0;
  }
 
  Axis *xaxis = m_outWS->getAxis(0);
  xaxis->setUnit("Time");
}
 
void ExportTimeSeriesLog::setupEventWorkspace(const size_t &start_index, const size_t &stop_index, int numentries,
                                              vector<DateAndTime> &times, vector<double> values,
                                              const bool &epochtime) {
  Types::Core::DateAndTime runstart(m_inputWS->run().getProperty("run_start")->value());
 
  // Get some stuff from the input workspace
  const size_t numberOfSpectra = 1;
  // determine output size
  size_t outsize = stop_index - start_index + 1;
  if (outsize > static_cast<size_t>(numentries))
    outsize = static_cast<size_t>(numentries);
 
  std::shared_ptr<EventWorkspace> outEventWS =
      create<EventWorkspace>(*m_inputWS, numberOfSpectra, HistogramData::BinEdges(2));
  m_outWS = outEventWS;
 
  // Create the output event list (empty)
  EventList &outEL = outEventWS->getSpectrum(0);
  outEL.switchTo(WEIGHTED_NOTIME);
 
  // Allocate all the required memory
  outEL.reserve(outsize);
  outEL.clearDetectorIDs();
 
  int64_t time_shift_ns(0);
  if (!epochtime) {
    // relative time
    time_shift_ns = runstart.totalNanoseconds();
  }
 
  for (size_t i = 0; i < outsize; i++) {
    Types::Core::DateAndTime tnow = times[i + start_index];
    int64_t dt = tnow.totalNanoseconds() - time_shift_ns;
 
    // convert to microseconds
    double dtmsec = static_cast<double>(dt) / 1000.0;
    outEL.addEventQuickly(WeightedEventNoTime(dtmsec, values[i + start_index], values[i + start_index]));
  }
  // Ensure thread-safety
  outEventWS->sortAll(TOF_SORT, nullptr);
 
  // Now, create a default X-vector for histogramming, with just 2 bins.
  std::vector<WeightedEventNoTime> &events = outEL.getWeightedEventsNoTime();
  outEventWS->setBinEdges(0, HistogramData::BinEdges{events.begin()->tof(), events.rbegin()->tof()});
}
 
bool ExportTimeSeriesLog::calculateTimeSeriesRangeByTime(std::vector<Types::Core::DateAndTime> &vec_times,
                                                         const double &rel_start_time, size_t &i_start,
                                                         const double &rel_stop_time, size_t &i_stop,
                                                         const double &time_factor) {
  // Initialize if there is something wrong.
  i_start = 0;
  i_stop = vec_times.size() - 1;
 
  // Check existence of proton_charge as run start
  Types::Core::DateAndTime run_start(0);
  if (m_inputWS->run().hasProperty("proton_charge")) {
    auto ts = dynamic_cast<TimeSeriesProperty<double> *>(m_inputWS->run().getProperty("proton_charge"));
    if (nullptr == ts) {
      throw std::runtime_error("Found the run property proton_charge but failed to interpret it "
                               "as a time series property of double values (failed dynamic cast).");
    }
    run_start = ts->nthTime(0);
  } else {
    g_log.warning("Property proton_charge does not exist so it is unable to "
                  "determine run start time. "
                  "StartTime and StopTime are ignored.  TimeSeriesProperty is "
                  "exported in full length.");
    return false;
  }
 
  // Get time 0
  if (rel_start_time != EMPTY_DBL()) {
    int64_t start_time_ns = run_start.totalNanoseconds() + static_cast<int64_t>(rel_start_time / time_factor);
    Types::Core::DateAndTime start_time(start_time_ns);
    i_start = static_cast<size_t>(std::lower_bound(vec_times.begin(), vec_times.end(), start_time) - vec_times.begin());
 
    // try to record the log value before starting time
    if (i_start >= 1)
      --i_start;
  }
 
  if (rel_stop_time != EMPTY_DBL()) {
    int64_t stop_time_ns = run_start.totalNanoseconds() + static_cast<int64_t>(rel_stop_time / time_factor);
    Types::Core::DateAndTime stop_time(stop_time_ns);
    i_stop = static_cast<size_t>(std::lower_bound(vec_times.begin(), vec_times.end(), stop_time) - vec_times.begin());
  }
 
  return true;
}
 
void ExportTimeSeriesLog::calculateFirstDerivative(bool is_event_ws) {
  if (is_event_ws) {
    g_log.error("It is not supported to calculate first derivative if the "
                "output is an EventWorkspace.");
    return;
  }
 
  // calcualte output
  size_t datasize = m_outWS->mutableX(1).size();
  auto vecX = m_outWS->mutableX(0);
  auto vecY = m_outWS->mutableY(0);
  auto &derivX = m_outWS->mutableX(1);
  auto &derivY = m_outWS->mutableY(1);
  if (vecY.size() != datasize)
    throw std::runtime_error("Output workspace 2D is not supposed to have "
                             "different size of X and Y.");
 
  std::stringstream errmsg_ss;
  for (size_t i = 1; i < datasize - 1; ++i) {
    // set up X
    derivX[i] = vecX[i];
    // set up Y
    double dx = vecX[i + 1] - vecX[i];
    if (dx <= 0) {
      errmsg_ss << "Entry " << i << ": " << vecX[i] << " >= " << vecX[i + 1] << "\n";
      derivY[i] = 1.;
    } else {
      derivY[i] = (vecY[i + 1] - vecY[i]) / dx;
    }
  }
 
  // last value
  derivX[datasize - 1] = vecX[datasize - 1];
  derivY.back() = 0.;
 
  // error message
  std::string errmsg = errmsg_ss.str();
  if (!errmsg.empty())
    g_log.error(errmsg);
 
  return;
}
 
void ExportTimeSeriesLog::setupMetaData(const std::string &log_name, const std::string &time_unit,
                                        const bool &export_epoch) {
 
  m_outWS->mutableRun().addProperty("SampleLogName", log_name, true);
  m_outWS->mutableRun().addProperty("TimeUnit", time_unit, true);
 
  std::string is_epoch("0");
  if (export_epoch)
    is_epoch = "1";
  m_outWS->mutableRun().addProperty("IsEpochTime", is_epoch, true);
}
 
} // namespace Mantid::Algorithms