FDiffFilter.hpp

Go to the documentation of this file.

//==============================================================================
//
//  This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
//
//  The GNSSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 3.0 of the License, or
//  any later version.
//
//  The GNSSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GNSSTk; if not, write to the Free Software Foundation,
//  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
//
//  This software was developed by Applied Research Laboratories at the
//  University of Texas at Austin.
//  Copyright 2004-2022, The Board of Regents of The University of Texas System
//
//==============================================================================
 
//==============================================================================
//
//  This software was developed by Applied Research Laboratories at the
//  University of Texas at Austin, under contract to an agency or agencies
//  within the U.S. Department of Defense. The U.S. Government retains all
//  rights to use, duplicate, distribute, disclose, or release this software.
//
//  Pursuant to DoD Directive 523024
//
//  DISTRIBUTION STATEMENT A: This software has been approved for public
//                            release, distribution is unlimited.
//
//==============================================================================
 
#ifndef FDIFF_FILTER_INCLUDE
#define FDIFF_FILTER_INCLUDE
 
#include "RobustStats.hpp"
#include "Stats.hpp"
#include "StatsFilterHit.hpp"
#include "StringUtils.hpp"
 
#include <vector>
 
namespace gnsstk
{
   //---------------------------------------------------------------------------------
   // TODO
 
   //---------------------------------------------------------------------------------
   template <class T> class IterativeFDiffFilter;
 
   //---------------------------------------------------------------------------------
   template <class T> class FDiffFilter
   {
   private:
      friend class IterativeFDiffFilter<T>;
 
      // embedded class -----------------------------------------
      class Analysis
      {
       public:
         Analysis()
         {
            index  = 0;
            sloInd = 0;
            diff = sigN = T(0);
         }
         // member data
         unsigned int
            index;   
         T diff;     
         T sigN;     
         T sloN;     
         int sloInd; 
      };             // end class Analysis
 
      std::vector<Analysis> Avec;
 
      std::vector<FilterHit<T>> results;
 
      // member data ---------------------------------------
      int osw, osp;  
      bool noxdata;  
      bool noflags;  
 
      const std::vector<T>
         &xdata;     
      const std::vector<T> &data; 
      const std::vector<int>
         &flags;     
      int ilimit;    
 
      T fdlim;            
      T siglim;           
      T Esiglim;          
      unsigned int Nwind; 
      unsigned int Nsig;  
      bool doSmall;       
      bool keepSigIndex;  
      std::vector<unsigned int>
         sigIndexes;      
 
      T medSlope, madSlope; 
 
   public:
      // member functions ---------------------------------------
      FDiffFilter(const std::vector<T>& x, const std::vector<T>& d,
                  const std::vector<int>& f)
         : data(d), xdata(x), flags(f)
      {
         fdlim        = T(0.8);
         siglim       = T(0.3);
         Esiglim      = T(0);
         noxdata      = (xdata.size() == 0);
         noflags      = (flags.size() == 0);
         osw          = 8;
         osp          = 3;
         Nsig         = 0;
         doSmall      = true;
         keepSigIndex = false;
      }
 
      inline void setWidth(unsigned int w) { Nwind = w; }
 
      inline unsigned int getWidth() { return Nwind; }
 
      inline void setLimit(T val) { fdlim = val; }
 
      inline T getLimit() { return fdlim; }
 
      inline void setSigma(T val) { siglim = val; }
 
      inline T getSigma() { return siglim; }
 
      inline bool doSmallSlips(bool doit)
      {
         doSmall = doit;
         return doit;
      }
 
      inline bool doSmallSlips() { return doSmall; }
 
      inline bool indexHighSigmas(bool doit)
      {
         keepSigIndex = doit;
         return doit;
      }
 
      inline bool indexingHighSigmas() { return keepSigIndex; }
 
      inline void setw(int w) { osw = w; }
 
      inline void setprecision(int p) { osp = p; }
 
      inline std::vector<FilterHit<T>> getResults() { return results; }
 
      inline unsigned int getNhighSigma() { return Nsig; }
 
      inline std::vector<unsigned int> getHighSigmaIndex()
      {
         return sigIndexes;
      }
 
      int filter(const size_t index = 0, int npts = -1);
 
      void computeRobustSigmaLimit(int& N, T& new_siglim);
 
      int analysis();
 
      void dump(std::ostream& s, std::string tag = std::string());
 
   }; // end class FDiffFilter
 
   //---------------------------------------------------------------------------------
   template <class T> int FDiffFilter<T>::filter(const size_t i0, int dsize)
   {
      if (dsize == -1)
      {
         dsize = data.size() - i0;
      }
 
      // largest allowed index is ilimit-1
      ilimit = dsize + i0;
 
      // is there enough data? if not return -1
      int i(i0), n(0);
      if (!noflags)
      {
         while (i < ilimit && n < 2)
         {
            if (flags[i] == 0)
            {
               n++;
            }
            i++;
         }
         if (i == ilimit)
         {
            return -1;
         }
      }
      else if (dsize < 2)
      {
         return -1;
      }
 
      // if xdata or flags is there, make sure its big enough
      if (!noflags && flags.size() - i0 < dsize)
      {
         return -3;
      }
 
      // generate the analysis vector
      Avec.clear();
 
      // compute stats on sigmas and data in a sliding window of width Nwind
      gnsstk::TwoSampleStats<T> fstats; // stats on the first diffs in window
      gnsstk::TwoSampleStats<T> dstats; // stats on the data in window
      std::vector<T> slopes;           // store slopes, for robust stats
 
      // loop over all data, computing first difference and stats in sliding
      // window
      Nsig = 0;
      int iprev(-1), j, islope(0);
      // start at the first good point
      n = 0;
      i = i0;
      if (!noflags)
      {
         while (i < ilimit && flags[i])
            i++;
      }
      while (i < ilimit)
      {
         Analysis A;
         A.index = i;
         n++; // count data points, just for approx slope when n==2
 
         // add to stats on data (for slope)
         dstats.Add(xdata[i], data[i]);
 
         // compute first difference, accounting for slope of data
         if (iprev == -1)
         {
            A.diff = T(0);
         }
         else
         {
            // get approx slope at first pt (often slope==0 here => OUT;
            // index=0)
            if (n == 2)
            {
               Avec[islope].sloN = // first slope = 2nd slope = d(data)/dx
                  (data[i] - data[iprev]) / (xdata[i] - xdata[iprev]);
            }
            // index of latest good slope
            A.sloInd = islope;
            // compute the difference = change in data - correction for slope
            A.diff = data[i] - data[iprev] -
                     Avec[islope].sloN * (xdata[i] - xdata[iprev]);
            // add diff to stats
            fstats.Add(xdata[i], A.diff);
         }
 
         // remove old data from stats buffers if full
         j = Avec.size() - Nwind; // index of earliest of the Nwind points
         if (fstats.N() > Nwind)
         {
            fstats.Subtract(xdata[Avec[j].index], Avec[j].diff);
         }
         if (dstats.N() > Nwind)
         {
            dstats.Subtract(xdata[Avec[j].index], data[Avec[j].index]);
         }
 
         // NO A.sigN = fstats.SigmaYX();     // sigma first diff, given slope
         // in fdiffs
         A.sigN = fstats.StdDevY(); // sigma first diff
         A.sloN = dstats.Slope();   // slope of data
         if (A.sigN > siglim)
         {
            Nsig++; // count it if sigma is high
            if (keepSigIndex)
            {
               sigIndexes.push_back(i);
            }
         }
         else
         {
            islope = Avec.size(); // keep this, the most recent good slope
            // keep slopes for roubst stats
            slopes.push_back(A.sloN);
         }
 
         Avec.push_back(A);
         iprev = i;
         i++;
         if (!noflags)
         {
            while (i < ilimit && flags[i])
               i++;
         }
      }
 
      // compute robust stats on slopes
      madSlope = gnsstk::Robust::MedianAbsoluteDeviation(
         &slopes[0], slopes.size(), medSlope, false);
 
      return Avec.size();
 
   } // end FDiffFilter::filter()
 
   //---------------------------------------------------------------------------------
   template <class T>
   void FDiffFilter<T>::computeRobustSigmaLimit(int& N, T& new_siglim)
   {
      new_siglim = siglim;
      if (Avec.size() < 2)
      {
         N = -1;
         return;
      }
 
      // compute high-outlier limit of sigmas using robust stats
      // compute quartiles - must sort first
      unsigned int i;
      std::vector<T> sd; // put sigmas in temp vector
      for (i = 0; i < Avec.size(); i++)
         sd.push_back(Avec[i].sigN);
 
      T Q1, Q3;
      gnsstk::QSort(&sd[0], sd.size());                     // sort
      gnsstk::Robust::Quartiles(&sd[0], sd.size(), Q1, Q3); // get Quartiles
 
      // compute new sigma limit ; outlier limit (high) 2.5Q3-1.5Q1
      new_siglim = 2.5 * Q3 - 1.5 * Q1;
      if (new_siglim <= siglim)
      {
         N = Nsig;
         return;
      }
 
      // count up hi-sigma points
      for (N = 0, i = 0; i < Avec.size(); i++)
         if (Avec[i].sigN > new_siglim)
         {
            N++;
         }
   }
 
   //---------------------------------------------------------------------------------
   template <class T> int FDiffFilter<T>::analysis()
   {
      // bool debug(true);
 
      // don't clear results - may be a case where caller wants to keep old ones
 
      // loop over analysis vector
      // outliers have >= Nwind big sigmas
      //      unless its the first point, then there are Nwind-2 bad sigmas
      // slips have exactly Nwind big sigmas
      bool isBad(false);
      int nbad(0), nout, n1st, Nw(Nwind), k; // DO NOT use unsigned
      unsigned int i, j, bad0;
      for (i = 0; i < Avec.size(); i++)
      {
         // if(debug)
         //   std::cout << "Loop " << i << " " << xdata[Avec[i].index] <<
         //   std::endl;
 
         if (Avec[i].sigN > siglim)
         {
            if (!isBad)
            {
               bad0  = i;
               isBad = true;
            }
            nbad++;
         }
 
         else if (isBad)
         {
            nout = int(bad0) + nbad - Nw;
            n1st = int(i) - Nw;
            // if(debug) std::cout << "FDFa: isBad at i " << i
            //      << " index " << Avec[i].index << " x " <<
            //      xdata[Avec[i].index]
            //      << " bad0 " << bad0 << " nbad " << nbad << " nout " << nout
            //      << " Nwind " << Nw << " 2Nw-2 " << 2*Nw-2 << std::endl;
 
            if (nout > 0 && nbad > Nw)
            { // outliers - more than a slip
               // if(debug) for(j=bad0; j<nout; j++)
               //   std::cout << "FDFa: OUT " << Avec[j].index
               //            << " " << xdata[Avec[j].index] << std::endl;
               j = bad0;
               FilterHit<T> fdfr;
               fdfr.type  = FilterHit<T>::outlier;
               fdfr.index = Avec[j].index;
               fdfr.npts  = nbad - Nw;
               fdfr.dx =
                  xdata[Avec[int(j) + nbad - Nw].index] - xdata[Avec[j].index];
               results.push_back(fdfr);
            }
 
            else if (nbad == Nw)
            { // slip
               j = bad0;
               // if(debug) std::cout << "FDFa: SLIP " << Avec[j].index <<
               // std::fixed
               //   << " " << xdata[Avec[j].index] << std::setprecision(osp)
               //   << " " << Avec[j].diff << std::endl;
               FilterHit<double> fdfr;
               fdfr.type  = FilterHit<T>::slip;
               fdfr.index = Avec[j].index;
               fdfr.sigma = Avec[j].sigN;
 
               // find number of missing pts before the slip
               k = Avec[j].index - 1; // j is an Avec[] index
               if (noflags && k > 0)
               {
                  k--; // NB k is int
               }
               else
               {
                  while (k > 0 && flags[k] != 0)
                     k--; // k is the data[] index
               }
               fdfr.dx   = xdata[Avec[j].index] - xdata[k];
               fdfr.npts = Avec[j].index - k;
 
               // get the step = first difference of data across slip
               fdfr.step = Avec[j].diff;
               // score here is just step/fdlim as a percentage (max 100)
               fdfr.score =
                  (unsigned int)(0.5 + 100. * ::fabs(fdfr.step) / fdlim);
               if (fdfr.score > 100)
               {
                  fdfr.score = 100;
               }
 
               // BUT if step is of order few*MADslope*dx, then probably a false
               // positive
               if (fdfr.score == 100 &&
                   ::fabs(fdfr.step) < 3 * madSlope * fdfr.dx)
               {
                  // if(debug)
                  // std::cout << "FDFa: F+ slip (|diff| < 3*madSlope*dx) at
                  // index "
                  //<< Avec[j].index << std::fixed << std::setprecision(osp)
                  //<< ": " << Avec[j].diff << " < 3 * " << madSlope
                  //<< " * " << fdfr.dx << " = " << 3*madSlope*fdfr.dx <<
                  //std::endl;
                  fdfr.score = -100;
               }
 
               // save the hit
               if (doSmall || fdfr.score == 100)
               {
                  results.push_back(fdfr);
               }
            }
 
            else if (int(i) <= 2 * Nw - 2 && n1st >= 0)
            { // first <Nwind pts outliers
               // if(debug) std::cout << "FDFa: isBad at small i " << i <<
               // std::endl; if(debug) for(int n=0; n<n1st; n++)
               //   std::cout << "FDFa OUT " << Avec[n].index << " "
               //            << xdata[Avec[n].index] << " small i " << i <<
               //            std::endl;
               FilterHit<double> fdfr;
               fdfr.type  = FilterHit<T>::outlier;
               fdfr.index = Avec[0].index;
               fdfr.npts  = n1st;
               fdfr.dx    = xdata[Avec[n1st].index] - xdata[Avec[0].index];
               results.push_back(fdfr);
            }
 
            isBad = false;
            nbad  = 0;
         }
      }
 
      // catch outliers at the very end
      if (isBad)
      { // bad at the very end
         // if(debug) for(j=bad0; j<Avec.size(); j++)
         //   std::cout << "FDFa: OUT " << Avec[j].index
         //            << " " << xdata[Avec[j].index] << std::endl;
         j = bad0;
         FilterHit<double> fdfr;
         fdfr.type  = FilterHit<T>::outlier;
         fdfr.index = Avec[j].index;
         fdfr.npts  = Avec.size() - j;
         fdfr.dx    = xdata[Avec[Avec.size() - 1].index] - xdata[Avec[j].index];
         results.push_back(fdfr);
      }
 
      return results.size();
 
   } // end FDiffFilter::analysis()
 
   //---------------------------------------------------------------------------------
   template <class T>
   void FDiffFilter<T>::dump(std::ostream& os, std::string tag)
   {
      size_t i, j, k, iprev;
      int w(osw > 5 ? osw + 1 : 5);
      os << "#" << tag << " FDiffFilter::dump() with limit " << std::fixed
         << std::setprecision(osp) << fdlim << " sigma limit " << siglim
         << std::setprecision(osp + 2) << " med,mad slope " << medSlope << " "
         << madSlope << std::setprecision(osp)
         << (noxdata ? " (xdata is index)" : "") << "\n#" << tag << " "
         << std::setw(2) << "i"
         << " " << std::setw(w) << "xd"
         << " " << std::setw(3) << "flg"
         << " " << std::setw(w) << "data"
         << " " << std::setw(w) << "fdif"
         << " " << std::setw(w) << "sig"
         << " " << std::setw(w) << "slope"
         << " " << std::setw(w) << "slp_u"
         << " " << std::setw(w) << "sl*dx"
         << " " << std::setw(w) << "slu*dx"
         << " " << std::setw(5)
         << "dx"
         //<< " " << std::setw(w) << "sigslope*dx"
         << std::endl;
 
      T dt(0);
      std::string sdif, ssig, sldx, slop, slou, sludx; //,sigslo;
      const size_t N(Avec.size());
      for (i = 0, j = 0, k = 0, iprev = -1; i < ilimit; i++)
      {
         // find j where Avec[j].index == i
         while (j < N && Avec[j].index < i)
            j++;
         bool haveAvec(Avec[j].index == i);
         if (haveAvec)
         {
            sdif = gnsstk::StringUtils::asString(Avec[j].diff, osp);
            ssig = gnsstk::StringUtils::asString(Avec[j].sigN, osp);
            if (iprev > -1)
            {
               dt = (noxdata ? T(i - iprev) : xdata[i] - xdata[iprev]);
            }
            slop = gnsstk::StringUtils::asString(Avec[j].sloN, osp);
            slou = gnsstk::StringUtils::asString(Avec[Avec[j].sloInd].sloN, osp);
            sldx = gnsstk::StringUtils::asString(Avec[j].sloN * dt, osp);
            sludx = gnsstk::StringUtils::asString(Avec[Avec[j].sloInd].sloN * dt,
                                                 osp);
            // sigslo = gnsstk::StringUtils::asString(madSlope*dt,osp);
         }
         else
         {
            sdif = ssig = sldx = sludx = slou = "?"; // sigslo = "?";
         }
         os << tag << std::fixed << std::setprecision(osp) << " "
            << std::setw(3) << i << " " << std::setw(w)
            << (noxdata ? T(i) : xdata[i]) << " " << std::setw(3)
            << (noflags ? 0 : flags[i]) << " " << std::setw(w) << data[i] << " "
            << std::setw(w) << sdif << " " << std::setw(w) << ssig << " "
            << std::setw(w) << slop << " " << std::setw(w) << slou << " "
            << std::setw(w) << sldx << " " << std::setw(w) << sludx << " "
            << std::setprecision(2) << std::setw(5)
            << dt
            //<< " " << std::setw(w) << sigslo
            << (haveAvec && Avec[j].sigN > siglim ? " SIG" : "")
            << (haveAvec ? "" : " NA");
         if (k < results.size() && haveAvec && i == results[k].index)
         {
            os << " " << results[k].asString()
               << (results[k].type == FilterHit<T>::slip &&
                         results[k].score < 100
                      ? " SMALL"
                      : "");
            k++;
         }
         os << std::endl;
         if (haveAvec)
         {
            iprev = Avec[j].index;
         }
      }
   } // end FDiffFilter<T>::dump()
 
   //---------------------------------------------------------------------------------
   //---------------------------------------------------------------------------------
   template <class T> class IterativeFDiffFilter
   {
   private:
      std::vector<FilterHit<T>> results;
 
      // member data ---------------------------------------
      // first those passes to FDiffFilter
      int osw, osp;          
 
      const std::vector<T>
         &xdata;             
      const std::vector<T>
         &data;              
      const std::vector<int>
         &flags;             
 
      T fdlim;               
      T siglim;              
      unsigned int Nwind;    
      unsigned int Nsig;     
      bool doSmall;          
 
      // and those unique to this class
      unsigned int itermax;  
      T Esiglim;             
      std::ostream &logstrm; 
      bool resetSigma;       
      T siguse;              
      bool keepSigIndex;     
      std::vector<unsigned int>
         sigIndexes;         
      bool verbose;          
      std::string label;     
 
   public:
      // member functions ---------------------------------------
      IterativeFDiffFilter(const std::vector<T>& x, const std::vector<T>& d,
                           const std::vector<int>& f,
                           std::ostream& os = std::cout)
         : data(d), xdata(x), flags(f), logstrm(os)
      {
         keepSigIndex = resetSigma = verbose = false;
         doSmall                             = true;
         itermax                             = 3;
         label                               = std::string();
         // rest are defaults of FDiffFilter
         FDiffFilter<T> fdf(x, d, f);
         osw = fdf.osw;
         osp = fdf.osp;
      }
 
      inline void setWidth(unsigned int w) { Nwind = w; }
 
      inline unsigned int getWidth() { return Nwind; }
 
      inline void setLimit(T val) { fdlim = val; }
 
      inline T getLimit() { return fdlim; }
 
      inline void setSigma(T val) { siglim = val; }
 
      inline T getSigma() { return siglim; }
 
      inline bool doSmallSlips(bool doit)
      {
         doSmall = doit;
         return doit;
      }
 
      inline bool doSmallSlips() { return doSmall; }
 
      inline bool doResetSigma(bool doit)
      {
         resetSigma = doit;
         return doit;
      }
 
      inline bool indexHighSigmas(bool doit)
      {
         keepSigIndex = doit;
         return doit;
      }
 
      inline bool indexingHighSigmas() { return keepSigIndex; }
 
      inline bool doVerbose(bool doit)
      {
         verbose = doit;
         return doit;
      }
 
      inline void setLabel(const std::string doit) { label = doit; }
 
      inline std::string getLabel() { return label; }
 
      inline void setw(int w) { osw = w; }
 
      inline void setprecision(int p) { osp = p; }
 
      inline T getEstimatedSigmaLimit() { return Esiglim; }
 
      inline T getUsedSigmaLimit() { return siguse; }
 
      std::vector<FilterHit<T>> getResults() { return results; }
 
      inline unsigned int getNhighSigma() { return Nsig; }
 
      inline std::vector<unsigned int> getHighSigmaIndex()
      {
         return sigIndexes;
      }
 
      int analysis();
 
      int editArrays(std::vector<T>& data, std::vector<int>& flags,
                     bool doInt = true, const int badFlag = 1);
 
   }; // end class IterativeFDiffFilter
 
   //---------------------------------------------------------------------------------
   template <class T> int IterativeFDiffFilter<T>::analysis()
   {
      // size of the data arrays
      const unsigned int size(xdata.size() < data.size() ? xdata.size()
                                                         : data.size());
      if (size <= 2)
      {
         return -1;
      }
 
      // save all results in each iteration; use results during each iteration
      std::vector<std::vector<FilterHit<T>>> Results;
 
      // handle input data ---------------------------------
      // use a copy of the data within the iteration loop
      std::vector<T> Tdata(data);
      // copy flags, create if needed
      std::vector<int> Tflags;
      if (flags.size() < size)
      {
         Tflags = std::vector<int>(size, 0);
      }
      else
      {
         Tflags = std::vector<int>(flags);
      }
 
      // analysis ------------------------------------------
      unsigned int iter, i, j, k, nr(0);
      siguse = siglim;
 
      // iterate over (filter / compute new sigma limit / analysis)
      iter = 1;
      while (iter <= itermax)
      {
         // must redefine filter each time since arrays (const in fdf) change
         FDiffFilter<T> fdf(xdata, Tdata, Tflags);
         fdf.setWidth(Nwind);       // fdf.Nwind
         fdf.setLimit(fdlim);       // fdf.fdlim
         fdf.setSigma(siguse);      // fdf.siglim
         fdf.setprecision(osp);     // fdf.osp
         fdf.setw(osw);             // fdf.osw
         fdf.doSmallSlips(doSmall); // fdf.doSmall
         fdf.indexHighSigmas(iter == itermax && keepSigIndex);
 
         // filter the data -----------
         i = fdf.filter();
         // if(verbose) logstrm << "# " << label << " filter returned " << i <<
         // std::endl;
         if (i <= 2)
         {
            logstrm << "Not enough data, abort.\n";
            return -1;
         }
 
         if (iter == itermax && keepSigIndex)
         {
            sigIndexes = fdf.getHighSigmaIndex();
         }
 
         // compute outlier limit from robust stats, and count outliers
         int N;
         fdf.computeRobustSigmaLimit(N, Esiglim);
         if (verbose)
         {
            logstrm << "# " << label << " Estimated sigma limit " << std::fixed
                    << std::setprecision(osp) << Esiglim
                    << " and used sigma limit " << siguse << " (input was "
                    << siglim << ") with " << N << " hi-sigma points "
                    << std::endl;
         }
 
         // reset sigma limit in filter, but not if its too large
         // Esiglim/siglim < 3.0            // not too big
         // Esiglim/siglim > 0.1            // not too small
         if (resetSigma)
         {
            if (Esiglim > siglim)
            {
               if (Esiglim / siglim < 3.0)
               {
                  siguse = Esiglim;
               }
               else
               {
                  siguse = 3 * siglim;
               }
            }
            else if (Esiglim < siglim)
            {
               if (Esiglim / siglim > 0.1)
               {
                  siguse = Esiglim;
               }
               else
               {
                  siguse = 0.1 * siglim;
               }
            }
            fdf.setSigma(siguse); // use the new sigma limit
         }
 
         // analysis ------------------
         fdf.analysis(); // get the outliers and slips
 
         // consider results ----------
         results = fdf.getResults();
         std::vector<unsigned int>
            eraseIndex; // index in results[] to be erased
 
         // loop over results: mark duplicates to be erased, mark outliers, fix
         // slips
         for (i = 0; i < results.size(); i++)
         {
            if (verbose)
            {
               logstrm << "# " << label << " Result " << ++nr << std::fixed
                       << std::setprecision(osp) << " "
                       << xdata[results[i].index] << " "
                       << results[i].asString(osp) << std::endl;
            }
 
            // mark outliers
            if (results[i].type == FilterHit<T>::outlier)
            {
               k = results[i].index;
               for (j = 0; j < results[i].npts; j++)
                  Tflags[k + j] = 1;
            }
 
            // slips: handle duplicates, and edit data
            if (results[i].type == FilterHit<T>::slip)
            {
 
               // search previous results, if they exist, for duplicate slips
               bool skip(false);
               for (k = 0; k < Results.size(); k++)
               {
 
                  for (j = 0; j < Results[k].size(); j++)
                  { // loop over prev results
                     FilterHit<T> &oldres(Results[k][j]);
 
                     if (oldres.type != FilterHit<T>::slip) // not a slip
                     {
                        continue;
                     }
                     if (results[i].index != oldres.index) // different index
                     {
                        continue;
                     }
 
                     // if both slips are non-small: add new to previous and
                     // delete new
                     if (results[i].score == 100 && oldres.score == 100)
                     {
                        // logstrm << " Combine slips " <<
                        // results[i].asString(osp)
                        // << " and " << oldres.asString(osp) << std::endl;
                        T step      = oldres.step + results[i].step;
                        oldres.step = step;
                        // keep the later sigma - outliers have been removed
                        oldres.sigma = results[i].sigma;
                        oldres.score =
                           (unsigned int)(0.5 + 100. * ::fabs(step) / fdlim);
                        if (oldres.score > 100)
                        {
                           oldres.score = 100;
                        }
                        // logstrm << " New combined slip is "
                        // << oldres.asString(osp) << std::endl;
                     }
 
                     // if both slips are small, they should be identical =>
                     // delete new
                     // ...so, if only one is small, go on
                     else if (results[i].score == 100 || oldres.score == 100)
                     {
                        continue;
                     }
 
                     // else both are small: delete new result
                     skip = true;
                     break;
                  }
               }
 
               // save the index, to be deleted later
               if (skip)
               {
                  eraseIndex.push_back(i);
                  continue;
               }
 
               // if slip is too small to fix, go on
               if (::fabs(results[i].step) < fdlim)
               {
                  // if(verbose) logstrm << "# " << label << " Slip too small: "
                  // << std::fixed << std::setprecision(osp) << results[i].step
                  //<< " < " << fdlim << std::endl;
                  continue;
               }
 
               // fix the slip in the Tdata
               int islip(results[i].step +
                         (results[i].step >= 0.0 ? 0.5 : -0.5));
               if (islip)
               {
                  if (verbose)
                  {
                     logstrm << "# " << label << " Fix slip " << islip << " "
                             << results[i].asString(osp) << std::endl;
                  }
 
                  // NB if slips in Tdata not fixed, its going to affect later
                  // runs of the analysis
                  for (j = results[i].index; j < Tdata.size(); j++)
                     Tdata[j] -= islip;
               }
            } // end if slip
 
         } // end loop over results
 
         // erase marked slips
         if (eraseIndex.size() > 0)
         {
            // sort(eraseIndex.begin(),eraseIndex.end());
            for (i = eraseIndex.size() - 1; i >= 0; i--)
            {
               // logstrm << " Erase " << results[eraseIndex[i]].asString(osp)
               //<< std::endl;
               results.erase(results.begin() + eraseIndex[i]); // erase vector
               if (i == 0)
               {
                  break; // i is unsigned
               }
            }
         }
 
         // save to Results
         Nsig = fdf.getNhighSigma();
         Results.push_back(results);
         results.clear();
 
         // dump this analysis
         if (verbose)
         {
            fdf.dump(logstrm,
                     "FIX" + gnsstk::StringUtils::asString(iter - 1) + label);
         }
 
         ++iter; // next iteration
 
      } // end iteration loop
 
      // copy all Results into results
      results.clear();
      for (k = 0; k < Results.size(); k++)
         for (j = 0; j < Results[k].size(); j++)
            results.push_back(Results[k][j]);
 
      // scan for remaining small slips (result of combination of two large)
      if (!doSmall)
      {
         std::vector<unsigned int>
            eraseIndex; // index in results[] to be erased
         for (i = 0; i < results.size(); i++)
         {
            if (results[i].type != FilterHit<T>::slip)
            {
               continue;
            }
            // logstrm << " Final check " << results[i].asString(osp) <<
            // std::endl;
            if (results[i].score == 100)
            {
               continue;
            }
            eraseIndex.push_back(i);
         } // end loop over results
         if (eraseIndex.size() > 0)
         {
            for (i = eraseIndex.size() - 1; i >= 0; i--)
            {
               // logstrm << " Final Erase " <<
               // results[eraseIndex[i]].asString(osp)
               //<< std::endl;
               results.erase(results.begin() + eraseIndex[i]); // erase vector
               if (i == 0)
               {
                  break; // i is unsigned
               }
            }
         }
         // logstrm << " End Final check" << std::endl;
      }
 
      return results.size();
   } // end IterativeFDiffFilter::analysis()
 
   //---------------------------------------------------------------------------------
   template <class T>
   int IterativeFDiffFilter<T>::editArrays(std::vector<T>& data,
                                           std::vector<int>& flags,
                                           bool doInt, const int badFlag)
   {
      int nedit(0);
      unsigned int i, j, n;
 
      for (i = 0; i < results.size(); i++)
      {
         if (results[i].type == FilterHit<T>::BOD ||
             results[i].type == FilterHit<T>::other)
         {
            continue;
         }
 
         else if (results[i].type == FilterHit<T>::outlier)
         {
            j = results[i].index;
            // TD should this loop account for flags[j] already set?
            for (n = 0; n < results[i].npts; nedit++, n++, j++)
               flags[j] = badFlag;
         }
 
         else if (results[i].type == FilterHit<T>::slip)
         {
            double slip(results[i].step);
            if (doInt)
            {
               int islip(slip + (slip >= 0.0 ? 0.5 : -0.5));
               if (islip == 0)
               {
                  continue;
               }
               slip = double(islip);
            }
            for (j = results[i].index; j < data.size(); j++)
               data[j] -= slip;
            nedit++;
         }
      }
 
      return nedit;
   } // end int IterativeFDiffFilter::editArrays()
 
   //---------------------------------------------------------------------------------
   //---------------------------------------------------------------------------------
} // namespace gnsstk
#endif // #define FDIFF_FILTER_INCLUDE