Skip to content
Snippets Groups Projects
TimeSeriesProperty.cpp 87.2 KiB
Newer Older
// 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
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/EmptyValues.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/TimeSplitter.h"
#include <json/value.h>
#include <nexus/NeXusFile.hpp>
#include <boost/regex.hpp>
#include <numeric>
namespace Mantid {
using namespace Types::Core;
namespace Kernel {
namespace {
/// static Logger definition
Logger g_log("TimeSeriesProperty");
} // namespace
/**
 * Constructor
 *  @param name :: The name to assign to the property
 */
template <typename TYPE>
TimeSeriesProperty<TYPE>::TimeSeriesProperty(const std::string &name)
    : Property(name, typeid(std::vector<TimeValueUnit<TYPE>>)), m_values(),
      m_size(), m_propSortedFlag(), m_filterApplied() {}

Tom Titcombe's avatar
Tom Titcombe committed
/**
 * Constructor
 * @param name :: The name to assign to the property
 * @param times :: A vector of DateAndTime objects
 * @param values :: A vector of TYPE
 */
template <typename TYPE>
TimeSeriesProperty<TYPE>::TimeSeriesProperty(
    const std::string &name, const std::vector<Types::Core::DateAndTime> &times,
    const std::vector<TYPE> &values)
    : TimeSeriesProperty(name) {
  addValues(times, values);
}

/// Virtual destructor
template <typename TYPE> TimeSeriesProperty<TYPE>::~TimeSeriesProperty() {}

/**
 * "Virtual" copy constructor
 */
template <typename TYPE>
TimeSeriesProperty<TYPE> *TimeSeriesProperty<TYPE>::clone() const {
  return new TimeSeriesProperty<TYPE>(*this);
}
/**
 * "Virutal copy constructor with a time shift
 * @param timeShift :: a time shift in seconds
 */
template <typename TYPE>
Property *
TimeSeriesProperty<TYPE>::cloneWithTimeShift(const double timeShift) const {
  auto timeSeriesProperty = this->clone();
  auto values = timeSeriesProperty->valuesAsVector();
  auto times = timeSeriesProperty->timesAsVector();
  // Shift the time
  for (auto it = times.begin(); it != times.end(); ++it) {
    // There is a known issue which can cause cloneWithTimeShift to be called
    // with a large (~9e+9 s) shift. Actual shifting is capped to be ~4.6e+19
    // seconds in DateAndTime::operator+=
    (*it) += timeShift;
  }
  timeSeriesProperty->clear();
  timeSeriesProperty->addValues(times, values);
  return timeSeriesProperty;
}

/** Return time series property, containing time derivative of current property.
 * The property itself and the returned time derivative become sorted by time
 * and the derivative is calculated in seconds^-1. (e.g. dValue/dT where
 * dT=t2-t1 is time difference in seconds for subsequent time readings and
 * dValue=Val1-Val2 is difference in subsequent values)
 *
 */
template <typename TYPE>
std::unique_ptr<TimeSeriesProperty<double>>
TimeSeriesProperty<TYPE>::getDerivative() const {
  if (this->m_values.size() < 2) {
    throw std::runtime_error("Derivative is not defined for a time-series "
                             "property with less then two values");
  auto it = this->m_values.begin();
  int64_t t0 = it->time().totalNanoseconds();
  TYPE v0 = it->value();

  it++;
  auto timeSeriesDeriv = std::make_unique<TimeSeriesProperty<double>>(
      this->name() + "_derivative");
  timeSeriesDeriv->reserve(this->m_values.size() - 1);
  for (; it != m_values.end(); it++) {
    TYPE v1 = it->value();
    int64_t t1 = it->time().totalNanoseconds();
    if (t1 != t0) {
      double deriv = 1.e+9 * (double(v1 - v0) / double(t1 - t0));
      auto tm = static_cast<int64_t>((t1 + t0) / 2);
      timeSeriesDeriv->addValue(Types::Core::DateAndTime(tm), deriv);
    t0 = t1;
    v0 = v1;
  }
  return timeSeriesDeriv;
}
/** time series derivative specialization for string type */
template <>
std::unique_ptr<TimeSeriesProperty<double>>
TimeSeriesProperty<std::string>::getDerivative() const {
  throw std::runtime_error(
      "Time series property derivative is not defined for strings");
/**
 * Return the memory used by the property, in bytes
 * */
template <typename TYPE>
size_t TimeSeriesProperty<TYPE>::getMemorySize() const {
  // Rough estimate
  return m_values.size() * (sizeof(TYPE) + sizeof(DateAndTime));
}
/**
 * Just returns the property (*this) unless overridden
 *  @param rhs a property that is merged in some descendent classes
 *  @return a property with the value
 */
template <typename TYPE>
TimeSeriesProperty<TYPE> &TimeSeriesProperty<TYPE>::merge(Property *rhs) {
  return operator+=(rhs);
}
/**
 * Add the value of another property
 * @param right the property to add
 * @return the sum
 */
template <typename TYPE>
TimeSeriesProperty<TYPE> &TimeSeriesProperty<TYPE>::
operator+=(Property const *right) {
  auto const *rhs = dynamic_cast<TimeSeriesProperty<TYPE> const *>(right);

  if (rhs) {
    if (this->operator!=(*rhs)) {
      m_values.insert(m_values.end(), rhs->m_values.begin(),
                      rhs->m_values.end());
      m_propSortedFlag = TimeSeriesSortStatus::TSUNKNOWN;
    } else {
      // Do nothing if appending yourself to yourself. The net result would be
      // the same anyway
      ;
    }

    // Count the REAL size.
    m_size = static_cast<int>(m_values.size());

  } else
    g_log.warning() << "TimeSeriesProperty " << this->name()
                    << " could not be added to another property of the same "
                       "name but incompatible type.\n";

  return *this;
}
/**
 * Deep comparison.
 * @param right The other property to compare to.
 * @return true if the are equal.
 */
template <typename TYPE>
bool TimeSeriesProperty<TYPE>::
operator==(const TimeSeriesProperty<TYPE> &right) const {

  if (this->name() != right.name()) // should this be done?
  {
    return false;
  }
  if (this->m_size != right.m_size) {
    return false;
  }
  if (this->realSize() != right.realSize()) {
    return false;
  } else {
    const std::vector<DateAndTime> lhsTimes = this->timesAsVector();
    const std::vector<DateAndTime> rhsTimes = right.timesAsVector();
    if (!std::equal(lhsTimes.begin(), lhsTimes.end(), rhsTimes.begin())) {
      return false;
    const std::vector<TYPE> lhsValues = this->valuesAsVector();
    const std::vector<TYPE> rhsValues = right.valuesAsVector();
    if (!std::equal(lhsValues.begin(), lhsValues.end(), rhsValues.begin())) {
      return false;
    }
  }
/**
 * Deep comparison.
 * @param right The other property to compare to.
 * @return true if the are equal.
 */
template <typename TYPE>
bool TimeSeriesProperty<TYPE>::operator==(const Property &right) const {
  auto rhs_tsp = dynamic_cast<const TimeSeriesProperty<TYPE> *>(&right);
  if (!rhs_tsp)
    return false;
  return this->operator==(*rhs_tsp);
}
/**
 * Deep comparison (not equal).
 * @param right The other property to compare to.
 * @return true if the are not equal.
 */
template <typename TYPE>
bool TimeSeriesProperty<TYPE>::
operator!=(const TimeSeriesProperty<TYPE> &right) const {
  return !(*this == right);
}
/**
 * Deep comparison (not equal).
 * @param right The other property to compare to.
 * @return true if the are not equal.
 */
template <typename TYPE>
bool TimeSeriesProperty<TYPE>::operator!=(const Property &right) const {
  return !(*this == right);
}
/**
 * Set name of the property
 */
template <typename TYPE>
void TimeSeriesProperty<TYPE>::setName(const std::string &name) {
  m_name = name;
}
/** Filter out a run by time. Takes out any TimeSeriesProperty log entries
 *outside of the given
 *  absolute time range.
 *  Be noticed that this operation is not reversible.
 *
 *  Use case 1: if start time of the filter fstart is in between t1 and t2 of
 *the TimeSeriesProperty,
 *              then, the new start time is fstart and the value of the log is
 *the log value @ t1
 *
 *  Use case 2: if the start time of the filter in on t1 or before log start
 *time t0, then
 *              the new start time is t1/t0/filter start time.
 *
 * EXCEPTION: If there is only one entry in the list, it is considered to mean
 * "constant" so the value is kept even if the time is outside the range.
 *
 * @param start :: Absolute start time. Any log entries at times >= to this time
 *are kept.
 * @param stop  :: Absolute stop time. Any log entries at times < than this time
 *are kept.
 */
template <typename TYPE>
LamarMoore's avatar
LamarMoore committed
void TimeSeriesProperty<TYPE>::filterByTime(
    const Types::Core::DateAndTime &start,
    const Types::Core::DateAndTime &stop) {

  // 1. Do nothing for single (constant) value
  if (m_values.size() <= 1)
    return;

  typename std::vector<TimeValueUnit<TYPE>>::iterator iterhead, iterend;

  // 2. Determine index for start and remove  Note erase is [...)
  int istart = this->findIndex(start);
  if (istart >= 0 && static_cast<size_t>(istart) < m_values.size()) {
    // "start time" is behind time-series's starting time
    iterhead = m_values.begin() + istart;

Lamar Moore's avatar
Lamar Moore committed
    // False - The filter time is on the mark.  Erase [begin(),  istart)
    // True - The filter time is larger than T[istart]. Erase[begin(), istart)
    // ...
    //       filter start(time) and move istart to filter startime
    bool useprefiltertime = !(m_values[istart].time() == start);

    // Remove the series
    m_values.erase(m_values.begin(), iterhead);

    if (useprefiltertime) {
      m_values[0].setTime(start);
    }
  } else {
    // "start time" is before/after time-series's starting time: do nothing
  // 3. Determine index for end and remove  Note erase is [...)
  int iend = this->findIndex(stop);
  if (static_cast<size_t>(iend) < m_values.size()) {
    if (m_values[iend].time() == stop) {
      // Filter stop is on a log.  Delete that log
      iterend = m_values.begin() + iend;
    } else {
      // Filter stop is behind iend. Keep iend
      iterend = m_values.begin() + iend + 1;
    }
    // Delete from [iend to mp.end)
    m_values.erase(iterend, m_values.end());
  }
  // 4. Make size consistent
  m_size = static_cast<int>(m_values.size());
}
/**
 * Filter by a range of times. If current property has a single value it remains
 * unaffected
 * @param splittervec :: A list of intervals to split filter on
 */
template <typename TYPE>
void TimeSeriesProperty<TYPE>::filterByTimes(
    const std::vector<SplittingInterval> &splittervec) {
  // 1. Sort

  // 2. Return for single value
  if (m_values.size() <= 1) {
    return;
  }
  // 3. Prepare a copy
  std::vector<TimeValueUnit<TYPE>> mp_copy;
  g_log.debug() << "DB541  mp_copy Size = " << mp_copy.size()
                << "  Original MP Size = " << m_values.size() << "\n";
  // 4. Create new
  for (const auto &splitter : splittervec) {
    Types::Core::DateAndTime t_start = splitter.start();
    Types::Core::DateAndTime t_stop = splitter.stop();
    int tstartindex = findIndex(t_start);
    if (tstartindex < 0) {
      // The splitter is not well defined, and use the first
      tstartindex = 0;
    } else if (tstartindex >= int(m_values.size())) {
      // The splitter is not well defined, adn use the last
      tstartindex = int(m_values.size()) - 1;
    int tstopindex = findIndex(t_stop);
    if (tstopindex < 0) {
      tstopindex = 0;
    } else if (tstopindex >= int(m_values.size())) {
      tstopindex = int(m_values.size()) - 1;
    } else {
      if (t_stop == m_values[size_t(tstopindex)].time() &&
          size_t(tstopindex) > 0) {
        tstopindex--;
    /* Check */
    if (tstartindex < 0 || tstopindex >= int(m_values.size())) {
      g_log.warning() << "Memory Leak In SplitbyTime!\n";
    }
    if (tstartindex == tstopindex) {
      TimeValueUnit<TYPE> temp(t_start, m_values[tstartindex].value());
      mp_copy.emplace_back(temp);
      mp_copy.emplace_back(t_start, m_values[tstartindex].value());
      for (auto im = size_t(tstartindex + 1); im <= size_t(tstopindex); ++im) {
        mp_copy.emplace_back(m_values[im].time(), m_values[im].value());
  g_log.debug() << "DB530  Filtered Log Size = " << mp_copy.size()
                << "  Original Log Size = " << m_values.size() << "\n";
  // 5. Clear
  m_values.clear();
  m_values = mp_copy;
  mp_copy.clear();
  m_size = static_cast<int>(m_values.size());
}
 * Split this time series property by time intervals to multiple time series
 * property according to number of distinct splitters' indexes, such as 0 and 1
 *
 * NOTE: If the input TSP has a single value, it is assumed to be a constant
 *  and so is not split, but simply copied to all outputs.
 *
 * @param splitter :: a TimeSplitterType object containing the list of intervals
 *                     and destinations.
 * @param outputs  :: A vector of output TimeSeriesProperty pointers of the same
 * @param isPeriodic :: whether the log (this TSP) is periodic. For example
 *                    proton-charge is periodic log.
 */
template <typename TYPE>
void TimeSeriesProperty<TYPE>::splitByTime(
    std::vector<SplittingInterval> &splitter, std::vector<Property *> outputs,
    bool isPeriodic) const {
  // 0. Sort if necessary

  if (outputs.empty())
    return;

  std::vector<TimeSeriesProperty<TYPE> *> outputs_tsp;

  size_t numOutputs = outputs.size();
  // 1. Clear the outputs before you start
  for (size_t i = 0; i < numOutputs; i++) {
    auto *myOutput = dynamic_cast<TimeSeriesProperty<TYPE> *>(outputs[i]);
    if (myOutput) {
      outputs_tsp.emplace_back(myOutput);
      if (this->m_values.size() == 1) {
        // Special case for TSP with a single entry = just copy.
        myOutput->m_values = this->m_values;
        myOutput->m_size = 1;
      } else {
        myOutput->m_values.clear();
        myOutput->m_size = 0;
      }
    } else {
      outputs_tsp.emplace_back(nullptr);
  // 2. Special case for TSP with a single entry = just copy.
  if (this->m_values.size() == 1)
    return;

  // 3. We will be iterating through all the entries in the the map/vector
  size_t i_property = 0;

  //    And at the same time, iterate through the splitter
  auto itspl = splitter.begin();
  size_t counter = 0;
  g_log.debug() << "[DB] Number of time series entries = " << m_values.size()
                << ", Number of splitters = " << splitter.size() << "\n";
  while (itspl != splitter.end() && i_property < m_values.size()) {
    // Get the splitting interval times and destination
    DateAndTime start = itspl->start();
    DateAndTime stop = itspl->stop();

    int output_index = itspl->index();
    // output workspace index is out of range. go to the next splitter
    if (output_index < 0 || output_index >= static_cast<int>(numOutputs))
      continue;

    TimeSeriesProperty<TYPE> *myOutput = outputs_tsp[output_index];
    // skip if the input property is of wrong type
    if (!myOutput) {
      ++itspl;
      ++counter;

    // Skip the events before the start of the time
    while (i_property < m_values.size() && m_values[i_property].time() < start)
      ++i_property;
    if (i_property == m_values.size()) {
      // i_property is out of the range. Then use the last entry
      myOutput->addValue(m_values[i_property - 1].time(),
                         m_values[i_property - 1].value());

      ++itspl;
      ++counter;
      break;
    }

    // The current entry is within an interval. Record them until out
    if (m_values[i_property].time() > start && i_property > 0 && !isPeriodic) {
      // Record the previous oneif this property is not exactly on start time
      //   and this entry is not recorded
      size_t i_prev = i_property - 1;
      if (myOutput->size() == 0 ||
          m_values[i_prev].time() != myOutput->lastTime())
        myOutput->addValue(m_values[i_prev].time(), m_values[i_prev].value());
    }

    // Loop through all the entries until out.
    while (i_property < m_values.size() && m_values[i_property].time() < stop) {

      // Copy the log out to the output
      myOutput->addValue(m_values[i_property].time(),
                         m_values[i_property].value());
      ++i_property;
    // Go to the next interval
    ++itspl;
    ++counter;
    // But if we reached the end, then we are done.
    if (itspl == splitter.end())
      break;
    // No need to keep looping through the filter if we are out of events
    if (i_property == this->m_values.size())
  } // Looping through entries in the splitter vector
  // Make sure all entries have the correct size recorded in m_size.
  for (std::size_t i = 0; i < numOutputs; i++) {
    auto *myOutput = dynamic_cast<TimeSeriesProperty<TYPE> *>(outputs[i]);
    if (myOutput) {
      myOutput->m_size = myOutput->realSize();
/// Split this TimeSeriresProperty by a vector of time with N entries,
/// and by the target workspace index defined by target_vec
/// Requirements:
/// vector outputs must be defined before this method is called
template <typename TYPE>
void TimeSeriesProperty<TYPE>::splitByTimeVector(
    std::vector<DateAndTime> &splitter_time_vec, std::vector<int> &target_vec,
    std::vector<TimeSeriesProperty *> outputs) {
  // check target vector to make it a set
  std::set<int> target_set;
  for (auto target : target_vec)
    target_set.insert(target);

  if (splitter_time_vec.size() != target_vec.size() + 1) {
    std::stringstream errss;
    errss << "Try to split TSP " << this->m_name
          << ": Input time vector's size " << splitter_time_vec.size()
          << " does not match (one more larger than) taget "
             "workspace index vector's size "
          << target_vec.size() << "\n";
    throw std::runtime_error(errss.str());
  }
  // return if the output vector TimeSeriesProperties is not defined
  // sort if necessary
  sortIfNecessary();

  // work on m_values, m_size, and m_time
  std::vector<Types::Core::DateAndTime> tsp_time_vec = this->timesAsVector();

  // go over both filter time vector and time series property time vector
  size_t index_splitter = 0;
  // tsp_time is start time of time series property
  DateAndTime tsp_time = tsp_time_vec[index_tsp_time];
  DateAndTime split_start_time = splitter_time_vec[index_splitter];
  DateAndTime split_stop_time = splitter_time_vec[index_splitter + 1];
  // move splitter index such that the first entry of TSP is before the stop
  // time of a splitter
  bool continue_search = true;
  bool no_entry_in_range = false;
  std::vector<DateAndTime>::iterator splitter_iter;
  splitter_iter = std::lower_bound(splitter_time_vec.begin(),
                                   splitter_time_vec.end(), tsp_time);
  if (splitter_iter == splitter_time_vec.begin()) {
    // do nothing as the first TimeSeriesProperty entry's time is before any
    // splitters
  } else if (splitter_iter == splitter_time_vec.end()) {
    // already search to the last splitter which is still earlier than first TSP
    // entry
    no_entry_in_range = true;
  } else {
    // calculate the splitter's index (now we check the stop time)
    index_splitter = splitter_iter - splitter_time_vec.begin() - 1;
    split_start_time = splitter_time_vec[index_splitter];
    split_stop_time = splitter_time_vec[index_splitter + 1];
  // move along the entries to find the entry inside the current splitter
  bool first_splitter_after_last_entry(false);
  if (!no_entry_in_range) {
    std::vector<DateAndTime>::iterator tsp_time_iter;
    tsp_time_iter = std::lower_bound(tsp_time_vec.begin(), tsp_time_vec.end(),
                                     split_start_time);
    if (tsp_time_iter == tsp_time_vec.end()) {
      // the first splitter's start time is LATER than the last TSP entry, then
      // there won't be any
      // TSP entry to be split into any target splitter.
      no_entry_in_range = true;
      first_splitter_after_last_entry = true;
    } else {
      // first splitter start time is between tsp_time_iter and the one before
      // it.
      // so the index for tsp_time_iter is the first TSP entry in the splitter
      index_tsp_time = tsp_time_iter - tsp_time_vec.begin();
  } else {
    // no entry in range is true, which corresponding to the previous case
    // "already search to the last splitter which is still earlier than first
    // TSP
    // entry"
  if (no_entry_in_range && first_splitter_after_last_entry) {
    // initialize all the splitters with the last value
    DateAndTime last_entry_time = this->lastTime();
    TYPE last_entry_value = this->lastValue();
    for (size_t i = 0; i < outputs.size(); ++i)
      outputs[i]->addValue(last_entry_time, last_entry_value);
  }
  // now it is the time to put TSP's entries to corresponding
  continue_search = !no_entry_in_range;
  size_t outer_while_counter = 0;
  bool partial_target_filled(false);
  while (continue_search) {
    int target = target_vec[index_splitter];

    // get the first entry index (overlap)
    if (index_tsp_time > 0)
      --index_tsp_time;
    // add the continous entries to same target time series property
    const size_t tspTimeVecSize = tsp_time_vec.size();
    bool continue_add = true;
    while (continue_add) {
      size_t inner_while_counter = 0;
      if (index_tsp_time == tspTimeVecSize) {
        // last entry. quit all loops
        continue_add = false;
        continue_search = false;
        partial_target_filled = true;
      // add current entry
      if (outputs[target]->size() == 0 ||
          outputs[target]->lastTime() < tsp_time_vec[index_tsp_time]) {
        // avoid to add duplicate entry
        outputs[target]->addValue(m_values[index_tsp_time].time(),
                                  m_values[index_tsp_time].value());
      }

      const size_t nextTspIndex = index_tsp_time + 1;
      if (nextTspIndex < tspTimeVecSize) {
        if (tsp_time_vec[nextTspIndex] > split_stop_time) {
          // next entry is out of this splitter: add the next one and quit
          if (outputs[target]->lastTime() < m_values[nextTspIndex].time()) {
            // avoid the duplicate cases occurred in fast frequency issue
            outputs[target]->addValue(m_values[nextTspIndex].time(),
                                      m_values[nextTspIndex].value());
          }
          // FIXME - in future, need to find out WHETHER there is way to
          // skip the
          // rest without going through the whole sequence
          continue_add = false;
    } // END-WHILE continue add

    // make splitters to advance to next
    ++index_splitter;
    if (index_splitter == splitter_time_vec.size() - 1) {
      // already last splitters
    } else {
      split_start_time = split_stop_time;
      split_stop_time = splitter_time_vec[index_splitter + 1];
  // Still in 'continue search'-while-loop.  But the TSP runs over before
  // splitters.
  // Therefore, the rest of the chopper must have one more entry added!
  if (partial_target_filled) {
    // fill the target
    std::set<int> fill_target_set;
    for (size_t isplitter = index_splitter;
         isplitter < splitter_time_vec.size() - 1; ++isplitter) {
      int target_i = target_vec[isplitter];
      if (fill_target_set.find(target_i) == fill_target_set.end()) {
        if (outputs[target_i]->size() == 0 ||
            outputs[target_i]->lastTime() != m_values.back().time())
          outputs[target_i]->addValue(m_values.back().time(),
                                      m_values.back().value());
        fill_target_set.insert(target_i);
        // quit loop if it goes over all the targets
        if (fill_target_set.size() == target_set.size())
          break;
      }
    }
  }

  // Add a debugging check such that there won't be any time entry with zero log
  for (size_t i = 0; i < outputs.size(); ++i) {
    if (outputs[i]->size() == 0) {
      std::stringstream errss;
      errss << i << "-th split-out term (out of " << outputs.size()
            << " total output TSP) of '" << m_name << "'' has "
            << outputs[i]->size() << " size, whose first entry is at "
            << this->firstTime().toSimpleString();
      g_log.debug(errss.str());
// The makeFilterByValue & expandFilterToRange methods generate a bunch of
// warnings when the template type is the wider integer types
// (when it's being assigned back to a double such as in a call to minValue or
// firstValue)
// However, in reality these methods are only used for TYPE=int or double (they
// are only called from FilterByLogValue) so suppress the warnings
#ifdef _WIN32
#pragma warning(push)
#pragma warning(disable : 4244)
#pragma warning(disable : 4804) // This one comes about for TYPE=bool - again
                                // the method is never called for this type
#endif
#if defined(__GNUC__) && !(defined(__INTEL_COMPILER))
#pragma GCC diagnostic ignored "-Wconversion"
/**
 * Fill a TimeSplitterType that will filter the events by matching
 * log values >= min and <= max. Creates SplittingInterval's where
 * times match the log values, and going to index==0.
 * This method is used by the FilterByLogValue algorithm.
 *
 * @param split :: Splitter that will be filled.
 * @param min :: min value
 * @param max :: max value
 * @param TimeTolerance :: offset added to times in seconds (default: 0)
 * @param centre :: Whether the log value time is considered centred or at the
 *beginning (the default).
 */
template <typename TYPE>
void TimeSeriesProperty<TYPE>::makeFilterByValue(
    std::vector<SplittingInterval> &split, double min, double max,
    double TimeTolerance, bool centre) const {
  const bool emptyMin = (min == EMPTY_DBL());
  const bool emptyMax = (max == EMPTY_DBL());

  if (!emptyMin && !emptyMax && max < min) {
    std::stringstream ss;
    ss << "TimeSeriesProperty::makeFilterByValue: 'max' argument must be "
          "greater than 'min' "
       << "(got min=" << min << " max=" << max << ")";
    throw std::invalid_argument(ss.str());
  }
  // If min or max were unset ("empty") in the algorithm, set to the min or max
  // value of the log
  if (emptyMin)
    min = minValue();
  if (emptyMax)
    max = maxValue();

  // Make sure the splitter starts out empty
  split.clear();

  // Do nothing if the log is empty.
  if (m_values.empty())
    return;

  // 1. Sort

  // 2. Do the rest
  bool lastGood(false);
  time_duration tol = DateAndTime::durationFromSeconds(TimeTolerance);
  int numgood = 0;
  DateAndTime t;
  DateAndTime start, stop;

  for (size_t i = 0; i < m_values.size(); ++i) {
    const DateAndTime lastTime = t;
    // The new entry
    t = m_values[i].time();
    TYPE val = m_values[i].value();

    // A good value?
    const bool isGood = ((val >= min) && (val <= max));
    if (isGood)
      numgood++;

    if (isGood != lastGood) {
      // We switched from bad to good or good to bad

      if (isGood) {
        // Start of a good section. Subtract tolerance from the time if
        // boundaries are centred.
        start = centre ? t - tol : t;
      } else {
        // End of the good section. Add tolerance to the LAST GOOD time if
        // boundaries are centred.
        // Otherwise, use the first 'bad' time.
        stop = centre ? lastTime + tol : t;
        split.emplace_back(start, stop, 0);
        // Reset the number of good ones, for next time
        numgood = 0;
      }
      lastGood = isGood;
  if (numgood > 0) {
    // The log ended on "good" so we need to close it using the last time we
    // found
    stop = t + tol;
    split.emplace_back(start, stop, 0);
/** Function specialization for TimeSeriesProperty<std::string>
 *  @throws Kernel::Exception::NotImplementedError always
 */
template <>
void TimeSeriesProperty<std::string>::makeFilterByValue(
    std::vector<SplittingInterval> & /*split*/, double /*min*/, double /*max*/,
    double /*TimeTolerance*/, bool /*centre*/) const {
  throw Exception::NotImplementedError("TimeSeriesProperty::makeFilterByValue "
                                       "is not implemented for string "
                                       "properties");
}
/** If the first and/or last values in a log are between min & max, expand and
 * existing TimeSplitter
 *  (created by makeFilterByValue) if necessary to cover the full TimeInterval
 * given.
 *  This method is used by the FilterByLogValue algorithm.
 *  @param split The splitter to modify if necessary
 *  @param min   The minimum 'good' value
 *  @param max   The maximum 'good' value
 *  @param range The full time range that we want this splitter to cover
 */
template <typename TYPE>
void TimeSeriesProperty<TYPE>::expandFilterToRange(
    std::vector<SplittingInterval> &split, double min, double max,
    const TimeInterval &range) const {
  const bool emptyMin = (min == EMPTY_DBL());
  const bool emptyMax = (max == EMPTY_DBL());

  if (!emptyMin && !emptyMax && max < min) {
    std::stringstream ss;
    ss << "TimeSeriesProperty::expandFilterToRange: 'max' argument must be "
          "greater than 'min' "
       << "(got min=" << min << " max=" << max << ")";
    throw std::invalid_argument(ss.str());
  }
  // If min or max were unset ("empty") in the algorithm, set to the min or max
  // value of the log
  if (emptyMin)
    min = minValue();
  if (emptyMax)
    max = maxValue();

  // Assume everything before the 1st value is constant
  double val = firstValue();
  if ((val >= min) && (val <= max)) {
    TimeSplitterType extraFilter;
    extraFilter.emplace_back(range.begin(), firstTime(), 0);
    // Include everything from the start of the run to the first time measured
    // (which may be a null time interval; this'll be ignored)
    split = split | extraFilter;
  }
  // Assume everything after the LAST value is constant
  val = lastValue();
  if ((val >= min) && (val <= max)) {
    TimeSplitterType extraFilter;
    extraFilter.emplace_back(lastTime(), range.end(), 0);
    // Include everything from the start of the run to the first time measured
    // (which may be a null time interval; this'll be ignored)
    split = split | extraFilter;
  }
}
/** Function specialization for TimeSeriesProperty<std::string>
 *  @throws Kernel::Exception::NotImplementedError always
 */
template <>
void TimeSeriesProperty<std::string>::expandFilterToRange(
    std::vector<SplittingInterval> & /*split*/, double /*min*/, double /*max*/,
    const TimeInterval & /*range*/) const {
  throw Exception::NotImplementedError("TimeSeriesProperty::makeFilterByValue "
                                       "is not implemented for string "
                                       "properties");
}
/** Calculates the time-weighted average of a property.
 *  @return The time-weighted average value of the log.
 */
template <typename TYPE>
double TimeSeriesProperty<TYPE>::timeAverageValue() const {
  double retVal = 0.0;
  try {
    const auto &filter = getSplittingIntervals();
    retVal = this->averageValueInFilter(filter);
  } catch (std::exception &) {
    // just return nan
    retVal = std::numeric_limits<double>::quiet_NaN();
  }
  return retVal;
}

/** Calculates the time-weighted average of a property in a filtered range.
 *  This is written for that case of logs whose values start at the times given.
 *  @param filter The splitter/filter restricting the range of values included
 *  @return The time-weighted average value of the log in the range within the
 * filter.
 */
template <typename TYPE>
double TimeSeriesProperty<TYPE>::averageValueInFilter(
    const std::vector<SplittingInterval> &filter) const {
  // TODO: Consider logs that aren't giving starting values.

  // First of all, if the log or the filter is empty, return NaN
  if (realSize() == 0 || filter.empty()) {
    return std::numeric_limits<double>::quiet_NaN();
  }
  // If there's just a single value in the log, return that.
  if (realSize() == 1) {
    return static_cast<double>(m_values.front().value());
  }

  double numerator(0.0), totalTime(0.0);
  // Loop through the filter ranges
  for (const auto &time : filter) {
    // Calculate the total time duration (in seconds) within by the filter
    totalTime += time.duration();

    // Get the log value and index at the start time of the filter
    int index;
    double value = getSingleValue(time.start(), index);
    DateAndTime startTime = time.start();
    while (index < realSize() - 1 && m_values[index + 1].time() < time.stop()) {
      ++index;
      numerator +=
          DateAndTime::secondsFromDuration(m_values[index].time() - startTime) *
          value;
      startTime = m_values[index].time();
      value = static_cast<double>(m_values[index].value());
    }

    // Now close off with the end of the current filter range
    numerator +=
        DateAndTime::secondsFromDuration(time.stop() - startTime) * value;
  // 'Normalise' by the total time
  return numerator / totalTime;
}