RdChannelLinearPredictorRFISuppressor.cc

Go to the documentation of this file.

#include "RdChannelLinearPredictorRFISuppressor.h"
#include "LinearPredictor.h"

#include <cmath>
#include <vector>
#include <algorithm>

#include <fwk/CentralConfig.h>

#include <utl/config.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/HannWindow.h>
#include <utl/StringCompare.h>
#include <utl/AugerException.h>

#include <evt/Event.h>
#include <evt/ShowerRRecDataQuantities.h>
#include <revt/StationRecData.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Channel.h>
#include <revt/Station.h>


#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "

using namespace revt;
using namespace fwk;
using namespace utl;
using namespace std;

namespace RdChannelLinearPredictorRFISuppressor {

  string
  SampleData::HelpOutput() {
    ostringstream helpoutput;
    helpoutput << "Useful values are:\n"
               << "  fBinning / ns = " << fBinning / nanosecond << "\n"
               << "All following units are samples: \n"
               << "  fNumCoefficients          = " << fNumCoefficients << "\n"
               << "  fDelayLine                = " << fDelayLine << "\n"
               << "  fTrainRegionStart         = " << fTrainRegionStart << "\n"
               << "  fTrainRegionStop          = " << fTrainRegionStop << "\n"
               << "  fMinimumTrainRegionLength = " << fMinimumTrainRegionLength << "\n"
               << "  fTraceLength              = " << fTraceLength << "\n"
               << "Occupied samples are the samples that can't be used because they are "
                  "in/close to the SignalSearchWindow or the NoiseWindow \n"
               << "  fOccupiedRegion1Start     = " << fOccupiedRegion1Start << "\n"
               << "  fOccupiedRegion1Stop      = " << fOccupiedRegion1Stop << "\n"
               << "  fOccupiedRegion2Start     = " << fOccupiedRegion2Start << "\n"
               << "  fOccupiedRegion2Stop      = " << fOccupiedRegion2Stop << "\n"
               << "  fTrainRegionLength        = " << fTrainRegionLength << "\n";
    return helpoutput.str();
  }


  unsigned int
  RdChannelLinearPredictorRFISuppressor::TranslateTimeToSample(double time, double binning, char const *name) {
    const double floatSamples = time / binning;
    const int signedintSamples = round(time/binning);
    const unsigned int intSamples = signedintSamples > 0 ? abs(signedintSamples) : 0;
    if ( fabs(floatSamples - intSamples) > 1E-10) {
      ostringstream warning;
      warning << "The configuration variable " << name << " with value "
              << time / nanosecond << "ns doesn't seem to be a multiple of the time-binning ("
              << binning / nanosecond << "ns) of this trace. Rounding to the nearest value of "
              << intSamples << "samples. This is not necessarily a problem.";
      INFOIntermediate(warning);
    }
    return intSamples;
  }


  bool
  RdChannelLinearPredictorRFISuppressor::TestOverlap(int a1, int a2, int b1, int b2) {
    // This method returns true if the region [a1...a2] and [b1...b2] overlap
    if (a1 > b1 && a1 < b2)
      return true;
    if (a2 > b1 && a2 < b2)
      return true;
    if (b1 > a1 && b1 < a2)
      return true;
    if (b2 > a1 && b2 < a2)
      return true;

    return false;
  }

  VModule::ResultFlag
  RdChannelLinearPredictorRFISuppressor::Init()
  {
    INFO("RdChannelLinearPredictorRFISuppressor::Init()");

    // Read in the configurations of the xml file
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelLinearPredictorRFISuppressor");
    topBranch.GetChild("SegmentLength").GetData(fSegmentLength);
    topBranch.GetChild("DelayLine").GetData(fDelayLine);
    topBranch.GetChild("TrainRegionStartOverride").GetData(fTrainRegionStartOverride);
    topBranch.GetChild("TrainRegionStopOverride").GetData(fTrainRegionStopOverride);
    topBranch.GetChild("MinimumTrainRegionLength").GetData(fMinimumTrainRegionLength);
    topBranch.GetChild("WindowMargins").GetData(fWindowMargins);
    topBranch.GetChild("FudgeFactor").GetData(fFudgeFactor);
    topBranch.GetChild("NumChannels").GetData(fNumChannels);
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("Optimization").GetData(fOptimization);

    if ( !( fDelayLine + fSegmentLength < (fTrainRegionStopOverride - fTrainRegionStartOverride)*2 ) ) {
      InvalidConfigurationException("The train-region must contain considerably more "
      "samples than those that are spanned by the coefficients and the delay line.");
    };

    return eSuccess;
  }

  void
  RdChannelLinearPredictorRFISuppressor::FillSampleDataWithRelevantValuesForFirstChannel(SampleData &samples, ChannelTimeSeries const &timeseries, Station const &station) {
    // The struct SampleData is filled with the relevant values and the TrainRegion is determined
    samples.fBinning                  = timeseries.GetBinning();
    samples.fTraceLength              = timeseries.GetSize();
    samples.fNumCoefficients          = TranslateTimeToSample(fSegmentLength, samples.fBinning, "SegmentLength");
    samples.fDelayLine                = TranslateTimeToSample(fDelayLine, samples.fBinning, "DelayLine");
    samples.fTrainRegionStart         = TranslateTimeToSample(fTrainRegionStartOverride, samples.fBinning, "RegionStart");
    samples.fTrainRegionStop          = TranslateTimeToSample(fTrainRegionStopOverride, samples.fBinning, "RegionStop");
    samples.fMinimumTrainRegionLength = TranslateTimeToSample(fMinimumTrainRegionLength, samples.fBinning, "MinimumTrainRegionLength");
    // If there is no override i.e. both sampleTrainRegionStart and sampleTrainRegionStop are zero
    // this part of the code tries to find the best train region outside the SignalSearchWindow and
    // the NoiseWindow
    // This module works on the channel level but some information from the station is necessary to determine and check the right train-region
    StationRecData const &recStation = station.GetRecData();
    double NoiseWindowStart = recStation.GetParameter(eNoiseWindowStart) - fWindowMargins;
    double NoiseWindowStop = recStation.GetParameter(eNoiseWindowStop) + fWindowMargins;
    double SignalSearchWindowStart = recStation.GetParameter(eSignalSearchWindowStart) - fWindowMargins;
    double SignalSearchWindowStop = recStation.GetParameter(eSignalSearchWindowStop) + fWindowMargins;
    if ( NoiseWindowStart < SignalSearchWindowStart ) {
      samples.fOccupiedRegion1Start = TranslateTimeToSample(NoiseWindowStart, samples.fBinning, "NoiseWindowStart");
      samples.fOccupiedRegion1Stop  = TranslateTimeToSample(NoiseWindowStop, samples.fBinning, "NoiseWindowStop");
      samples.fOccupiedRegion2Start = TranslateTimeToSample(SignalSearchWindowStart, samples.fBinning, "SignalSearchWindowStart");
      samples.fOccupiedRegion2Stop  = TranslateTimeToSample(SignalSearchWindowStop, samples.fBinning, "SignalSearchWindowStop");
    } else {
      samples.fOccupiedRegion2Start = TranslateTimeToSample(NoiseWindowStart, samples.fBinning, "NoiseWindowStart");
      samples.fOccupiedRegion2Stop  = TranslateTimeToSample(NoiseWindowStop, samples.fBinning, "NoiseWindowStop");
      samples.fOccupiedRegion1Start = TranslateTimeToSample(SignalSearchWindowStart, samples.fBinning, "SignalSearchWindowStart");
      samples.fOccupiedRegion1Stop  = TranslateTimeToSample(SignalSearchWindowStop, samples.fBinning, "SignalSearchWindowStop");
    }

    // This code is written to automatically pick the best piece of available noise
    if (samples.fTrainRegionStart == 0 and samples.fTrainRegionStop == 0) {
      unsigned int Range1Start = timeseries.GetStart();
      unsigned int Range1Stop  = samples.fOccupiedRegion1Start;
      unsigned int Range2Start = samples.fOccupiedRegion1Stop;
      unsigned int Range2Stop  = samples.fOccupiedRegion2Start;
      unsigned int Range3Start = samples.fOccupiedRegion2Stop;
      unsigned int Range3Stop  = timeseries.GetStop() + 1;
      int Range1Length = int(Range1Stop) - int(Range1Start); // negative values are possible but if all but the region that is chosen should definitiely not be negative
      int Range2Length = int(Range2Stop) - int(Range2Start);
      int Range3Length = int(Range3Stop) - int(Range3Start);
      if ((Range1Length > Range2Length) && (Range1Length > Range3Length)) {
        samples.fTrainRegionStart = Range1Start;
        samples.fTrainRegionStop = Range1Stop;
      } else if ((Range2Length > Range1Length) && (Range2Length > Range3Length)) {
        samples.fTrainRegionStart = Range2Start;
        samples.fTrainRegionStop = Range2Stop;
      } else {
        samples.fTrainRegionStart = Range3Start;
        samples.fTrainRegionStop = Range3Stop;
      }
    }
    samples.fTotalTraceLength = timeseries.GetStop() - timeseries.GetStart() + 1;
    samples.fTrainRegionLength = samples.fTrainRegionStop - samples.fTrainRegionStart;
  }


  VModule::ResultFlag
  RdChannelLinearPredictorRFISuppressor::TestSampleDataForConsistency(SampleData &samples, ChannelTimeSeries const &timeseries) {
    // Test for overlaps of the regions. No regions should overlap
    if ( TestOverlap(samples.fTrainRegionStart, samples.fTrainRegionStop, samples.fOccupiedRegion1Start, samples.fOccupiedRegion1Stop)
      || TestOverlap(samples.fTrainRegionStart, samples.fTrainRegionStop, samples.fOccupiedRegion2Start, samples.fOccupiedRegion2Stop)
    ) {
      ostringstream error;
      error << "There is overlap of the TrainRegion with the SignalSearchWindow or the SignalNoiseWindow "
        "Maybe you have set the train-region to a wrong manual value?). Here are some values to help you" << endl
        << samples.HelpOutput() << endl;
      ERROR(error.str());
      return eFailure;
    }
    // Test that that the length of the Train region is larger or equal to the minimum specified length
    if ( samples.fTrainRegionLength < int(samples.fMinimumTrainRegionLength) ) {
      ostringstream error;
      error << "The number of samples in the train region (" << samples.fTrainRegionLength <<
        ") is lower than the minimum number of train-region samples (" <<
        samples.fMinimumTrainRegionLength << "). Here are some values to help you:" << endl
        << samples.HelpOutput() << endl;
      ERROR(error.str());
      return eFailure;
    }
    // Test that the train-region-start sample is in the valid part of the time-series
    if ( samples.fTrainRegionStart < timeseries.GetStart() ) {
      ostringstream error;
      error << "Lowest sample of the noise region (" << samples.fTrainRegionStart <<
        ") is lower than the lowest valid sample of the trace (" <<
        timeseries.GetStart() << "). Here are some values to help you:" << endl
        << samples.HelpOutput() << endl;
      ERROR(error.str());
      return eFailure;
    }
    // Test that the train-region-stop sample is in the valid part of the time-series
    if ( samples.fTrainRegionStop > timeseries.GetStop() + 1 ) {
      ostringstream error;
      error << "Highest sample of the noise region (" << samples.fTrainRegionStop - 1 <<
        ") is higher than the highest valid sample of the trace (" <<
        timeseries.GetStop() << "). Here are some values to help you:" << endl
        << samples.HelpOutput() << endl;
      ERROR(error.str());
      return eFailure;
    }
    // Test that the binning of the channels is the same
    if ( samples.fBinning != timeseries.GetBinning() ) {
      ERROR("Something bad is going on. The traces that I am asked to work with have "
        "different binning. I can't work with that.");
      return eFailure;
    }
    // Test that the length of the channels is the same (if not then a warning is emitted an attempt is made to continue)
    if ( samples.fTraceLength != timeseries.GetSize() ) {
      WARNING("Something bad might be going on. The traces of one station not have the "
        "same length. I am adujsting the total length to the shortest trace.");
      if ( timeseries.GetSize() < samples.fTraceLength ) {
        samples.fTraceLength = timeseries.GetSize();
      }
    }
    // Test that the the train region length is at least larger than the delay-line plus number of coefficients
    // in general, for the algorithm to work well, it should be at least several times larger but if it is
    // than this value then the algorithm would fail completely.
    if ( !( int(samples.fDelayLine + samples.fNumCoefficients) < samples.fTrainRegionLength ) ) {
      ostringstream error;
      error << "This is not good: based on this sampling "
        "the train-region does not contain more samples than those that are spanned "
        "by the coefficients and the delay line. It would be impossible to calculate "
        "the covariances and the coefficients. Here are some values to help you:" << endl
        << samples.HelpOutput() << endl;
      ERROR(error.str());
      return eFailure;
    }
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelLinearPredictorRFISuppressor::Run(evt::Event& event)
  {
    INFO("RdChannelLinearPredictorRFISuppressor::Run()");

    // Check whether there are any events at all //
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();

    // Initialization of the values such as the binning and the regions //
    // loop through stations and for every station through every channel
    for (auto& station : rEvent.StationsRange()) {

      unsigned int countActiveChannels = 0;
      // Create a new struct instance of SampleData which will contain the relevant information
      // about the samples
      SampleData samples;

      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        ChannelTimeSeries& timeseries(channel.GetChannelTimeSeries());
        if (countActiveChannels == 0) {
          // Most of the information about the time-series is extracted from a single channel
          // for instance the binning must be the same for all channels, otherwise the predictor
          // won't work
          FillSampleDataWithRelevantValuesForFirstChannel(samples, timeseries, station);
        }
        // The values are determined from the first channel tested and for consistency
        // for all channels here. If, for instance, the binning is not the same for all
        // channels then an error is emitted here.
        // In addition some other checks are done such as to test for acceptable size
        // of the regions and to test that regions do not overlap
        VModule::ResultFlag test_samples_result = TestSampleDataForConsistency(samples, timeseries);
        if (test_samples_result != eSuccess) {
          return test_samples_result;
        }
        ++countActiveChannels;
      }
      if ( countActiveChannels == 0 ) {
        WARNING("No active channels were found. Skipping.");
        return eContinueLoop;
      }
      if ( countActiveChannels != fNumChannels && fNumChannels != 0 ) {
        ostringstream error;
        error << "The value of fNumChannels (" << fNumChannels << ") is not the same as the number of active channels ("
              << countActiveChannels << ")" << endl;
        ERROR(error.str());
        return eFailure;
      }

      // Some useful Output //
      OUT(2) << "--- Working on event with Id=" << rEvent.GetHeader().GetId()
             << " and station with Id=" << station.GetId() << "---" << endl;
      OUT(2) << "Setting up the linear predictor with:" << endl
             << "      - train region         : [" << samples.fTrainRegionStart << ", " << samples.fTrainRegionStop << ") samples" << endl
             << "      - number of coefficients: " << samples.fNumCoefficients << " samples" << endl
             << "      - delay line            : " << samples.fDelayLine << " samples" << endl
             << "      - total trace length    : " << samples.fTotalTraceLength << " samples" << endl;
      OUT(4) << samples.HelpOutput() << endl;

      // Here we apply the algorithm and use the LinearPredictor class //
      // We prepare the noise-region and the predict-region of the traces and put them into c++ vectors.
      unsigned int channelNum = 0;
      vector< vector<double> > noiseRegion(countActiveChannels, vector<double>(samples.fTrainRegionLength));
      vector< vector<double> > predictRegion(countActiveChannels, vector<double>(samples.fTraceLength));

      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        ChannelTimeSeries &timeseries(channel.GetChannelTimeSeries());
        for (unsigned int index = 0; int(index) < samples.fTrainRegionLength; ++index) {
          noiseRegion[channelNum][index] = timeseries[index + samples.fTrainRegionStart];
        }
        for (unsigned int index = 0; index < samples.fTraceLength; ++index) {
          predictRegion[channelNum][index] = timeseries[index];
        }
        ++channelNum;
      }

      // Here we train the linear predictor with the noise-region and
      // let it make a prediction based on the predict-region
      LinearPredictor lp(samples.fNumCoefficients, samples.fDelayLine, fNumChannels, fOptimization);

      OUT(2) << "Training the linear predictor..." << endl;
      lp.train(noiseRegion, fFudgeFactor);
      vector< vector<double> > prediction(lp.predict(predictRegion));

      OUT(2) << "Subtracting the predicted signal from the traces..." << endl;
      // We now have the predicted RFI values. These are subtracted from the
      // trace to yield a trace that contains less RFI. (i.e. all 'predictable'
      // RFI is removed)
      channelNum = 0;
      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        ChannelTimeSeries &timeseries(channel.GetChannelTimeSeries());
        for (unsigned int index = 0; index < samples.fTraceLength; ++index) {
          if ( index < samples.fDelayLine + samples.fNumCoefficients + 1 ) {
            timeseries[index] = 0;
          } else {
            timeseries[index] -= prediction[channelNum][index];
          }
        }
        ++channelNum;
      }
    }
    return eSuccess;
  }

  VModule::ResultFlag
  RdChannelLinearPredictorRFISuppressor::Finish()
  {
    // Put any termination or cleanup code here.
    // This method is called once at the end of the run.

    INFO("RdChannelLinearPredictorRFISuppressor::Finish()");

    return eSuccess;
  }

}