TrackCheckShip.cc

Go to the documentation of this file.

/*
 * (C) Copyright 2020 Met Office UK
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
 
#include "ufo/filters/TrackCheckShip.h"
 
#include <Eigen/Core>  // for easy array addition/subtraction
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <map>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
 
#include "eckit/config/Configuration.h"
#include "ioda/ObsDataVector.h"
#include "ioda/ObsSpace.h"
#include "oops/base/Variables.h"
#include "oops/util/DateTime.h"
#include "oops/util/Duration.h"
#include "oops/util/Logger.h"
#include "ufo/filters/TrackCheckShipDiagnostics.h"
#include "ufo/filters/TrackCheckShipParameters.h"
#include "ufo/utils/Constants.h"
#include "ufo/utils/RecursiveSplitter.h"
 
namespace ufo {
 
TrackCheckShip::TrackCheckShip(ioda::ObsSpace &obsdb, const eckit::Configuration &config,
                               std::shared_ptr<ioda::ObsDataVector<int> > flags,
                               std::shared_ptr<ioda::ObsDataVector<float> > obserr)
  : FilterBase(obsdb, config, flags, obserr)
{
  oops::Log::debug() << "TrackCheckShip: config = " << config << std::endl;
  options_.reset(new TrackCheckShipParameters());
  options_->deserialize(config);
  assert(options_->maxSpeed.value() > 0);
  if (options_->testingMode)
    diagnostics_.reset(new TrackCheckShipDiagnostics());
}
 
/// \brief Estimate of speed as calculated in Ops_CheckShipTrack.
///
/// Speed is calculated between \p obs1 and \p obs2, accounting for spatial/temporal
/// uncertainty using the resolution values stored in \p options.
double TrackCheckShip::speedEstimate(
    const TrackCheckShip::TrackObservation &obs1,
    const TrackCheckShip::TrackObservation &obs2,
    const TrackCheckShipParameters& options) {
  util::Duration temporalDistance = abs(obs1.getTime() -
                                        obs2.getTime());
  util::Duration tempRes = options.temporalResolution;
  auto dist = distance(obs1, obs2);
  auto spatialRes = options.spatialResolution;
  double speedEst = 0.0;
  if (dist > spatialRes) {
    speedEst = (dist - 0.5 * spatialRes) /
        std::max(temporalDistance.toSeconds(),
                 (tempRes.toSeconds()));
  }
  return speedEst;
}
 
/// \brief Returns the angle in degrees (rounded to the nearest .1 degree)
/// formed by the track segment going from observation \p a to \p b to \p c.
float TrackCheckShip::angle(const TrackCheckShip::TrackObservation &a,
                            const TrackCheckShip::TrackObservation &b,
                            const TrackCheckShip::TrackObservation &c) {
  auto locA = a.getLocation();
  auto locB = b.getLocation();
  auto locC = c.getLocation();
  Eigen::Vector3f disp1{locB[0]-locA[0], locB[1]-locA[1], locB[2]-locA[2]};
  Eigen::Vector3f disp2{locC[0]-locB[0], locC[1]-locB[1], locC[2]-locB[2]};
  auto retValue = std::acos(disp1.dot(disp2) / (disp1.norm() * disp2.norm()));
  retValue *= static_cast<float>(Constants::rad2deg);
  return std::round(retValue * 10.0f) * 0.1f;
}
 
const TrackCheckShipDiagnostics* TrackCheckShip::diagnostics() const
{
  return diagnostics_.get();
}
 
TrackCheckShip::TrackObservation::TrackObservation(
    double latitude, double longitude,
    const util::DateTime &time,
    std::shared_ptr<TrackStatistics> const& ts,
    std::shared_ptr<TrackCheckUtils::CheckCounter> const &checkCounter,
    size_t observationNumber)
  : obsLocationTime_(latitude, longitude, time), fullTrackStatistics_(ts),
    checkCounter_(checkCounter), observationNumber_(observationNumber) {}
 
void TrackCheckShip::TrackObservation::registerCheckResult(
    const TrackCheckUtils::CheckResult &result) {
  checkCounter_->registerCheckResult(result);
}
 
size_t TrackCheckShip::TrackObservation::getObservationNumber() const
{
  return observationNumber_;
}
 
// Required for the correct destruction of options_.
TrackCheckShip::~TrackCheckShip()
{}
 
/// Detects whether each observation in a track has been rejected, and marks
/// the corresponding space in \p isRejected as true if so.
/// \param trackObsIndicesBegin the begin iterator of the track
/// \param trackObsIndicesEnd the end iterator
/// \param validObsIds the vector marking all of the observations' global
/// positions
/// \param trackObservations the full vector of observations within
/// the single track
/// \param isRejected the boolean vector whose indices correspond to all of
/// the rejected observations in the full input dataset
void TrackCheckShip::flagRejectedTrackObservations(
    std::vector<size_t>::const_iterator trackObsIndicesBegin,
    std::vector<size_t>::const_iterator trackObsIndicesEnd,
    const std::vector<size_t> &validObsIds,
    const std::vector<TrackObservation> &trackObservations,
    std::vector<bool> &isRejected) const {
  auto trackObsIndexIt = trackObsIndicesBegin;
  auto trackObsIt = trackObservations.begin();
  for (; trackObsIndexIt != trackObsIndicesEnd; ++trackObsIndexIt, ++trackObsIt)
    if (trackObsIt->rejected())
      isRejected[validObsIds[*trackObsIndexIt]] = true;
}
 
void TrackCheckShip::print(std::ostream & os) const {
  os << "TrackCheckShip: config = " << config_ << std::endl;
}
 
 
/// \todo This section is not yet fully implemented. Current implementation includes separating
/// observations into tracks based on \p Station_Id, calculating distances, speeds, and
/// angles between observations, and incrementing track-specific counters should the
/// calculations produce unexpected values. Based on the counter values, the filter
/// may be stopped early. Simultaneous observations are then optionally ignored for later
/// re-inclusion or rejection. Using the above calculations and additional ones if necessary,
/// the observation directly before and/or after the fastest segment is chosen for rejection.
/// Re-checking of simultaneous observations is not yet implemented.
 
void TrackCheckShip::applyFilter(const std::vector<bool> & apply,
                                 const Variables & filtervars,
                                 std::vector<std::vector<bool>> & flagged) const {
  const std::vector<size_t> validObsIds = TrackCheckUtils::getValidObservationIds(apply, flags_);
 
  RecursiveSplitter splitter(validObsIds.size());
  TrackCheckUtils::groupObservationsByStation(validObsIds, splitter, config_, obsdb_);
  TrackCheckUtils::sortTracksChronologically(validObsIds, splitter, obsdb_);
 
  TrackCheckUtils::ObsGroupLocationTimes obsLocTime =
      TrackCheckUtils::collectObservationsLocations(obsdb_);
 
  std::vector<bool> isRejected(apply.size(), false);
  for (auto track : splitter.multiElementGroups()) {
    std::vector<TrackObservation> trackObservations = collectTrackObservations(
          track.begin(), track.end(), validObsIds, obsLocTime);
    std::vector<std::reference_wrapper<TrackObservation>> trackObservationsReferences;
    trackObservationsReferences.reserve(trackObservations.size());
    std::transform(trackObservations.begin(), trackObservations.end(),
                   std::back_inserter(trackObservationsReferences),
                   [](TrackObservation &obs) {
      return std::ref<TrackObservation>(obs); });
    calculateTrackSegmentProperties(trackObservationsReferences,
                                    CalculationMethod::FIRSTITERATION);
    if (!trackObservationsReferences.empty() &&
        this->options_->earlyBreakCheck &&
        TrackCheckShip::earlyBreak(trackObservationsReferences)) {
      continue;
    }
    if (options_->deferredCheckSimultaneous.value()) {
      trackObservationsReferences = removeSimultaneousObservations(trackObservationsReferences);
    } else {  // reinstating deferred observation is not yet implemented
      bool firstIterativeRemoval = true;
      while (trackObservationsReferences.size() >= 3) {
        // Initial loop: fastest (as determined by set of comparisons) observation removed
        // until all segments show slower speed than max threshold
        auto maxSpeedReferenceIterator = std::max_element(
              trackObservationsReferences.begin(), trackObservationsReferences.end(),
              [](TrackObservation a, TrackObservation b) {
            return a.getObservationStatistics().speed <
            b.getObservationStatistics().speed;});
        auto maxSpeedValue = maxSpeedReferenceIterator->get().getObservationStatistics().speed;
        if (maxSpeedValue <= (0.8 * options_->maxSpeed.value())) {
          break;
        } else if (maxSpeedValue < options_->maxSpeed.value()) {
          auto maxSpeedAngle = std::max(
                (maxSpeedReferenceIterator - 1)->get().getObservationStatistics().angle,
                maxSpeedReferenceIterator->get().getObservationStatistics().angle);
          if (maxSpeedAngle <= 90.0) {
            break;
          }
        }
        removeFaultyObservation(
              trackObservationsReferences, maxSpeedReferenceIterator, firstIterativeRemoval);
        firstIterativeRemoval = false;
        calculateTrackSegmentProperties(trackObservationsReferences, CalculationMethod::MAINLOOP);
      }
      flagRejectedTrackObservations(track.begin(), track.end(),
                                    validObsIds, trackObservations, isRejected);
    }
  }
  TrackCheckUtils::flagRejectedObservations(isRejected, flagged);
}
 
/// \returns a \p vector of \p TrackObservations that all hold a \p shared_ptr to an instance
/// of \p TrackStatistics, which holds all of the track-specific counters.
std::vector<TrackCheckShip::TrackObservation> TrackCheckShip::collectTrackObservations(
    std::vector<size_t>::const_iterator trackObsIndicesBegin,
    std::vector<size_t>::const_iterator trackObsIndicesEnd,
    const std::vector<size_t> &validObsIds,
    const TrackCheckUtils::ObsGroupLocationTimes &obsLocTime) const {
  std::vector<TrackObservation> trackObservations;
  trackObservations.reserve(trackObsIndicesEnd - trackObsIndicesBegin);
  std::shared_ptr<TrackStatistics> trackStatistics(new TrackStatistics());
  std::shared_ptr<TrackCheckUtils::CheckCounter> checkCounter(new TrackCheckUtils::CheckCounter);
  size_t observationNumber = 0;
  for (std::vector<size_t>::const_iterator it = trackObsIndicesBegin;
       it != trackObsIndicesEnd; ++it) {
    const size_t obsId = validObsIds[*it];
    trackObservations.push_back(TrackObservation(obsLocTime.latitudes[obsId],
                                                 obsLocTime.longitudes[obsId],
                                                 obsLocTime.datetimes[obsId], trackStatistics,
                                                 checkCounter,
                                                 observationNumber));
    observationNumber++;
  }
  return trackObservations;
}
 
/// \brief \returns true if at least half of the track segments have
/// incremented the relevant rejection counters
///
/// This filter is best for mostly-accurate observation tracks. If the track has
/// many "errors" (as indicated by the counters that are incremented before
/// any observations are removed), it will stop early by default.
bool TrackCheckShip::earlyBreak(const std::vector<std::reference_wrapper<TrackObservation>>
                                &trackObs) const {
  bool breakResult = false;
  const auto& trackStats = *(trackObs[0].get().getFullTrackStatistics());
  // if at least half of the track segments have a time difference of less than an hour
  // (if non-buoy), are faster than a configured maximum speed, or exhibit at least a 90
  // degree bend
  if ((2 * ((options_->inputCategory.value() != 1)  // 1 is input category of buoy
            * trackStats.numShort_ + trackStats.numFast_) + trackStats.numBends_)
      >= (trackObs.size() - 1)) {
    std::string stationId = "no station id provided";
    if (options_->stationIdVariable.value() != boost::none) {
      stationId = (options_->stationIdVariable.value().get()).variable();
    }
    oops::Log::warning() << "ShipTrackCheck: " << stationId << "\n" <<
                            "Time difference < 1 hour: " << trackStats.numShort_ << "\n" <<
                            "Fast: " << trackStats.numFast_ << "\n" <<
                            "Bends: " << trackStats.numBends_ << "\n" <<
                            "Track was not checked." << std::endl;
 
    breakResult = true;
  }
  if (options_->testingMode)
    diagnostics_->storeEarlyBreakResult(breakResult);
  return breakResult;
}
 
/// \brief Filters out all observations that have been marked simultaneous for
/// rechecking at the end.
std::vector<std::reference_wrapper<TrackCheckShip::TrackObservation>>
TrackCheckShip::removeSimultaneousObservations(
    const std::vector<std::reference_wrapper<TrackObservation>> &trackObs) const {
  std::vector<std::reference_wrapper<TrackObservation>> reducedTrackObs;
  std::copy_if(trackObs.begin(), trackObs.end(), std::back_inserter(reducedTrackObs),
               [](const TrackObservation &obs){
    return !(obs.getObservationStatistics().simultaneous);});
  calculateTrackSegmentProperties(reducedTrackObs, CalculationMethod::SIMULTANEOUSDEFERRAL);
  return reducedTrackObs;
}
 
/// \brief Chooses which of the observations surrounding the speediest segment to remove,
/// flagging it accordingly.
void TrackCheckShip::removeFaultyObservation(
    std::vector<std::reference_wrapper<TrackObservation>> &track,
    const std::vector<std::reference_wrapper<TrackObservation>>::iterator
    &observationAfterFastestSegment,
    bool firstIterativeRemoval) const {
  int errorCategory = 0;
  util::Duration four_days{"P4D"};
  auto rejectedObservation = observationAfterFastestSegment;
  // lambda function to "fail" an observation that should be rejected
  auto fail = [&rejectedObservation](
      const std::vector<std::reference_wrapper<TrackObservation>>::iterator &iter) {
    iter->get().registerCheckResult(TrackCheckUtils::CheckResult::FAILED);
    iter->get().registerSweepOutcome(true);
    rejectedObservation = iter;
  };
  auto neighborObservationStatistics = [&observationAfterFastestSegment](int index) ->
      const ObservationStatistics & {
    return (observationAfterFastestSegment + index)->get().getObservationStatistics();
  };
  auto meanSpeed = observationAfterFastestSegment->get().getFullTrackStatistics()->meanSpeed_;
  if (observationAfterFastestSegment == track.begin() + 1) {
    // Decide whether ob 0 or 1 agrees best with ob 2
    if (neighborObservationStatistics(0).speedAveraged <=
        options_->maxSpeed &&
        (neighborObservationStatistics(1).speed >
         options_->maxSpeed ||
         neighborObservationStatistics(1).angle > 45.0)) {
      fail(observationAfterFastestSegment);
      errorCategory = 2;
    } else {
      fail(observationAfterFastestSegment - 1);
      errorCategory = 1;
    }
  } else if (observationAfterFastestSegment == track.end() - 1) {
    if (neighborObservationStatistics(-1).speedAveraged <=
        options_->maxSpeed &&
        (neighborObservationStatistics(-1).speed >
         options_->maxSpeed ||
         neighborObservationStatistics(-2).angle > 45.0)) {
      fail(observationAfterFastestSegment - 1);
      errorCategory = 2;
    } else {
      fail(observationAfterFastestSegment);
      errorCategory = 1;
    }
  } else if (neighborObservationStatistics(-1).speed >
             options_->maxSpeed) {
    fail(observationAfterFastestSegment - 1);
    errorCategory = 4;
    //  Category 4: both segments surrounding observation have excessive speed
  } else if (neighborObservationStatistics(1).speed >
             options_->maxSpeed) {
    fail(observationAfterFastestSegment);
    errorCategory = 4;
  } else if (neighborObservationStatistics(0).speedAveraged >
             options_->maxSpeed) {
    fail(observationAfterFastestSegment - 1);
    errorCategory = 5;
    // Category 5: observation before fastest segment would still begin a fast segment
    // if observation after fastest segment were removed
  } else if (neighborObservationStatistics(-1).speedAveraged >
             options_->maxSpeed) {
    fail(observationAfterFastestSegment);
    errorCategory = 5;
    // Category 5: observation after fastest segment would still end a fast segment if
    // observation before fastest segment were removed
  } else if (neighborObservationStatistics(-1).angle >
             (45.0 + neighborObservationStatistics(0).angle)) {
    fail(observationAfterFastestSegment - 1);
    errorCategory = 6;
  } else if (neighborObservationStatistics(0).angle >
             45.0 + neighborObservationStatistics(-1).angle) {
    fail(observationAfterFastestSegment);
    errorCategory = 6;
  } else if (neighborObservationStatistics(-2).angle > 45.0 &&
             neighborObservationStatistics(-2).angle >
             neighborObservationStatistics(1).angle) {
    fail(observationAfterFastestSegment - 1);
    errorCategory = 7;
  } else if (neighborObservationStatistics(1).angle > 45.0) {
    fail(observationAfterFastestSegment);
    errorCategory = 7;
  } else if (neighborObservationStatistics(-1).speed <
             0.5 * std::min(
               neighborObservationStatistics(1).speed,
               meanSpeed)) {
    fail(observationAfterFastestSegment - 1);
    errorCategory = 8;
  } else if (neighborObservationStatistics(1).speed <
             0.5 * std::min(neighborObservationStatistics(-1).speed, meanSpeed)) {
    fail(observationAfterFastestSegment);
    errorCategory = 8;
  } else {
    double distanceSum = 0.0;
    distanceSum = std::accumulate(observationAfterFastestSegment - 1,
                    observationAfterFastestSegment + 2, distanceSum,
                    [](double summed, TrackObservation& obs) {
      return summed + obs.getObservationStatistics().distance;
    });
 
    double distancePrevObsOmitted =
          neighborObservationStatistics(-1).distanceAveraged +
        neighborObservationStatistics(1).distance;
    double distanceCurrentObsOmitted =
          neighborObservationStatistics(-1).distance +
        neighborObservationStatistics(0).distanceAveraged;
    util::Duration timeSum = (observationAfterFastestSegment + 1)->get().getTime() - (
          observationAfterFastestSegment - 2)->get().getTime();
    if (options_->testingMode.value()) {
      diagnostics_->storeDistanceSum(distanceSum);
      diagnostics_->storeDistancePrevObsOmitted(distancePrevObsOmitted);
      diagnostics_->storeDistanceCurrentObsOmitted(distanceCurrentObsOmitted);
      double timeDouble = timeSum.toSeconds();
      diagnostics_->storeTimeSum(timeDouble);
    }
    if (distancePrevObsOmitted < distanceCurrentObsOmitted - std::max(
          options_->spatialResolution.value(), 0.1 * distanceSum)) {
      fail(observationAfterFastestSegment - 1);
      errorCategory = 9;
    } else if (distanceCurrentObsOmitted < (
                 distancePrevObsOmitted - std::max(
                   options_->spatialResolution.value(), 0.1 * distanceSum))) {
      fail(observationAfterFastestSegment);
      errorCategory = 9;
    } else if (timeSum <= four_days && timeSum.toSeconds() > 0 &&
               std::min(distancePrevObsOmitted, distanceCurrentObsOmitted) > 0.0) {
      double previousSegmentDistanceProportion =
          // Prev segment dist/(prev segment distance + distAveragedCurrentObservation)
          neighborObservationStatistics(-1).distance /
          distanceCurrentObsOmitted;
      double previousObservationDistanceAveragedProportion =
          // Previous observation distAveraged / (prev obs distAveraged + next segment distance)
          neighborObservationStatistics(-1).
          distanceAveraged / distancePrevObsOmitted;
      double previousSegmentTimeProportion =
          static_cast<double>(((observationAfterFastestSegment - 1)->get().getTime() -
           (observationAfterFastestSegment - 2)->get().getTime()).toSeconds()) /
          timeSum.toSeconds();
      double previousAndFastestSegmentTimeProportion =
          static_cast<double>(((observationAfterFastestSegment->get().getTime()) -
            (observationAfterFastestSegment - 2)->get().getTime()).toSeconds()) /
          timeSum.toSeconds();
      if (options_->testingMode.value()) {
        diagnostics_->storePreviousSegmentDistanceProportion(previousSegmentDistanceProportion);
        diagnostics_->storePreviousObservationDistanceAveragedProportion(
              previousObservationDistanceAveragedProportion);
        diagnostics_->storePreviousSegmentTimeProportion(previousSegmentTimeProportion);
        diagnostics_->storePreviousAndFastestSegmentTimeProportion(
              previousAndFastestSegmentTimeProportion);
      }
      if (std::abs(previousSegmentDistanceProportion - previousSegmentTimeProportion) >
          0.1 + std::abs(previousObservationDistanceAveragedProportion -
                         previousAndFastestSegmentTimeProportion)) {
        // previous segment's spatial and temporal lengths are significantly more disproportionate
        // when compared to a larger portion of track
        // than the equivalents for the post-fastest segment
        fail(observationAfterFastestSegment - 1);
        errorCategory = 10;
      } else if (std::abs(previousObservationDistanceAveragedProportion -
                          previousAndFastestSegmentTimeProportion) > 0.1 +
                 std::abs(previousSegmentDistanceProportion - previousSegmentTimeProportion)) {
        // next segment after fastest has spatial and temporal lengths that are significantly more
        // disproportionate than the segment before the fastest
        fail(observationAfterFastestSegment);
        errorCategory = 10;
      } else {
        fail(observationAfterFastestSegment);
      }
      oops::Log::trace() << "CheckShipTrack: proportions " << previousSegmentDistanceProportion <<
                            " " << previousSegmentTimeProportion <<
                            " " << previousObservationDistanceAveragedProportion << " "
                         << previousAndFastestSegmentTimeProportion << " speeds: " << meanSpeed
                         << " " << neighborObservationStatistics(-1).speed << " " <<
                            neighborObservationStatistics(0).speed << " " <<
                            neighborObservationStatistics(1).speed << " [m/s]" << std::endl;
    }
    if (errorCategory == 9 || std::min(distancePrevObsOmitted, distanceCurrentObsOmitted) == 0.0) {
      oops::Log::trace() << "CheckShipTrack: Dist check, station id: " <<
                            options_->stationIdVariable.value().get().variable() << std::endl <<
                            " error category: " << errorCategory << std::endl <<
                            " distances: " << distanceSum * 0.001 << " " <<
                            distancePrevObsOmitted * 0.001 << " " <<
                            distanceCurrentObsOmitted * 0.001 << " " <<
                            (distancePrevObsOmitted - distanceCurrentObsOmitted) * 0.001 <<
                            " " << std::max(options_->spatialResolution.value(),
                                              0.1 * distanceSum) *
                            0.001 << "[km]" << std::endl;
    }
  }
  if (errorCategory == 0 || ((rejectedObservation->get().getObservationStatistics().
                              speedAveraged) >
                             options_->maxSpeed.value())) {
    oops::Log::trace() << "CheckShipTrack: cannot decide between observation and previous, " <<
                          "rejecting both." << std::endl;
    errorCategory += 100;
    if (options_->testingMode.value() && firstIterativeRemoval) {
      std::vector<size_t> observationNumbersAroundFastest{
        (observationAfterFastestSegment - 1)->get().getObservationNumber(),
            observationAfterFastestSegment->get().getObservationNumber()};
      diagnostics_->storeFirstIterativeRemovalInfo(
            std::make_pair(observationNumbersAroundFastest, errorCategory));
    }
    if (rejectedObservation == observationAfterFastestSegment) {
      fail(observationAfterFastestSegment - 1);
    } else {
      fail(observationAfterFastestSegment);
    }
    track.erase(observationAfterFastestSegment - 1, observationAfterFastestSegment);
  } else {
    if (options_->testingMode.value() && firstIterativeRemoval) {
      std::vector<size_t> rejectedObservationNumber{rejectedObservation->get().
            getObservationNumber()};
      diagnostics_->storeFirstIterativeRemovalInfo(
            std::make_pair(rejectedObservationNumber,
                           errorCategory));
    }
    track.erase(rejectedObservation);
  }
  oops::Log::trace() << "Error category: " << errorCategory << std::endl;
}
/// \todo Trace output will need to be changed to match that of OPS (indices, LWin)
 
/// \brief Calculates all of the statistics that require only two
/// adjacent \p TrackObservations, storing within the righthand observation.
///
/// This includes
/// distance between the two observations,
/// time difference between the observations, speed between the
/// observations, and if the
/// observations are recorded for the same time. Calls function to
/// increment track-wise counters based on
/// results.
void TrackCheckShip::TrackObservation::calculateTwoObservationValues(
    TrackObservation& prevObs, bool firstIteration,
    const TrackCheckShipParameters &options) {
  this->setDistance(TrackCheckShip::distance(prevObs, *this));
  (this->observationStatistics_.distance > options.spatialResolution) ?
        this->setSpeed(TrackCheckShip::speedEstimate(*this, prevObs, options)) :
        this->setSpeed(0.0);
  this->setTimeDifference(getTime() - prevObs.getTime());
  if (!(this->observationStatistics_.timeDifference.toSeconds())) {
    if (this->observationStatistics_.distance >
        options.spatialResolution) {
      prevObs.setSimultaneous(true);
    }
    setSimultaneous(true);
  }
  if (firstIteration) {
    adjustTwoObservationStatistics(options);
  }
}
 
/// Calculates all of the statistics that require three
/// adjacent \p TrackObservations,
/// storing within the middle observation. This includes
/// distance between two alternating observations,
/// speed between these alternating observations (if the middle
/// observation was not recorded), and the
/// angle formed by the three-observation track segment.
/// Calls function to increment track-wise counters based on results.
void TrackCheckShip::TrackObservation::calculateThreeObservationValues(
    const TrackObservation& prevObs, const TrackObservation& nextObs,
    bool firstIteration, const TrackCheckShipParameters &options) {
  this->setDistanceAveraged(TrackCheckShip::distance(prevObs, nextObs));
  this->setSpeedAveraged(speedEstimate(prevObs, nextObs, options));
  if (std::min(this->observationStatistics_.distance,
               nextObs.observationStatistics_.distance) >
      options.spatialResolution) {
    this->observationStatistics_.angle = angle(prevObs, *this, nextObs);
  }
  if (firstIteration) {
    adjustThreeObservationStatistics();
  }
}
 
/// This performs all of the necessary calculations based on the
/// observations' locations and timestamps by calling
/// \p calculateTwoObservationValues and \p calculateThreeObservationValues for the
/// non-edge-case observations.
void TrackCheckShip::calculateTrackSegmentProperties(
    const std::vector<std::reference_wrapper<TrackObservation>> &trackObservations,
    CalculationMethod calculationMethod) const {
  if (trackObservations.size()) {
    for (size_t obsIdx = 1; obsIdx < trackObservations.size(); obsIdx++) {
      TrackObservation &obs = trackObservations[obsIdx].get();
      TrackObservation &prevObs = trackObservations[obsIdx - 1].get();
      obs.calculateTwoObservationValues(prevObs, calculationMethod == FIRSTITERATION, *options_);
      if (obsIdx > 1) {
        const TrackObservation &prevPrevObs = trackObservations[obsIdx - 2].get();
        prevObs.calculateThreeObservationValues(prevPrevObs, obs,
                                                calculationMethod == FIRSTITERATION, *options_);
      }
      if (calculationMethod == FIRSTITERATION && (obsIdx == trackObservations.size() - 1)) {
        int potentialDenominator = trackObservations.size() - 1 -
            obs.getFullTrackStatistics()->numShort_ - obs.getFullTrackStatistics()->numFast_;
        (obs.getFullTrackStatistics()->meanSpeed_) = (obs.getFullTrackStatistics()->sumSpeed_) /
            std::max(1, potentialDenominator);
      }
    }
    if (options_->testingMode.value() && calculationMethod != MAINLOOP) {
      std::vector<TrackCheckShip::ObservationStatistics> obsStats;
      for (size_t obsIdx = 0; obsIdx < trackObservations.size(); ++obsIdx) {
        obsStats.push_back(trackObservations[obsIdx].get().getObservationStatistics());
      }
      auto trackStats = *(trackObservations[0].get().getFullTrackStatistics());
      if (calculationMethod == FIRSTITERATION)
        diagnostics_->storeInitialCalculationResults(std::make_pair(obsStats, trackStats));
      else if (calculationMethod == SIMULTANEOUSDEFERRAL)
        diagnostics_->storeCalculatedResultsSimultaneousDeferred(obsStats);
    }
  }
}
 
/// \brief Keeps track of 0-distanced, short, and fast track segments,
/// as well as incrementing \p sumSpeed_ for normal track segments.
void TrackCheckShip::TrackObservation::adjustTwoObservationStatistics
(const TrackCheckShipParameters &options) const {
  util::Duration hour{"PT1H"};
  if (getObservationStatistics().distance == 0.0)
    getFullTrackStatistics()->numDist0_++;
  if (getObservationStatistics().timeDifference < hour) {
    getFullTrackStatistics()->numShort_++;
  } else if (getObservationStatistics().speed >= options.maxSpeed) {
    getFullTrackStatistics()->numFast_++;
  } else {
    getFullTrackStatistics()->sumSpeed_ += getObservationStatistics().speed;
  }
}
 
/// \brief Increments number of bends if angle is greater than or equal to 90 degrees
void TrackCheckShip::TrackObservation::adjustThreeObservationStatistics() const {
  if (getObservationStatistics().angle >= 90.0)
    getFullTrackStatistics()->numBends_++;
}
 
const TrackCheckShip::ObservationStatistics &
TrackCheckShip::TrackObservation::getObservationStatistics() const
{
  return observationStatistics_;
}
 
const std::shared_ptr<TrackCheckShip::TrackStatistics>
TrackCheckShip::TrackObservation::getFullTrackStatistics() const
{
  return fullTrackStatistics_;
}
 
/// Sets the calling observation as simultaneous and increments the track-wise
/// \p numSimultaneous_ counter.
void TrackCheckShip::TrackObservation::setSimultaneous(bool simul) {
  assert(simul);  // reverting a simultaneous obs to normal is not yet implemented
  if (!this->observationStatistics_.simultaneous)
    this->fullTrackStatistics_->numSimultaneous_++;
  this->observationStatistics_.simultaneous = simul;
}
void TrackCheckShip::TrackObservation::setDistance(double dist) {
  this->observationStatistics_.distance = dist;
}
void TrackCheckShip::TrackObservation::setTimeDifference(util::Duration tDiff) {
  this->observationStatistics_.timeDifference = tDiff;
}
void TrackCheckShip::TrackObservation::setSpeed(double speed) {
  this->observationStatistics_.speed = speed;
}
void TrackCheckShip::TrackObservation::setAngle(double angle) {
  this->observationStatistics_.angle = angle;
}
void TrackCheckShip::TrackObservation::setDistanceAveraged(double distAvg) {
  this->observationStatistics_.distanceAveraged = distAvg;
}
void TrackCheckShip::TrackObservation::setSpeedAveraged(double speedAvg) {
  this->observationStatistics_.speedAveraged = speedAvg;
}
}  // namespace ufo