SimpleExtender.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2012.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Clemens Groepl $
// $Authors: $
// --------------------------------------------------------------------------

#ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_SIMPLEEXTENDER_H
#define OPENMS_TRANSFORMATIONS_FEATUREFINDER_SIMPLEEXTENDER_H

#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/FeaFiModule.h>
#include <OpenMS/MATH/STATISTICS/AveragePosition.h>

#include <queue>
#include <iostream>
#include <fstream>

namespace OpenMS
{

  template <class PeakType, class FeatureType>
  class SimpleExtender :
    public FeaFiModule<PeakType, FeatureType>,
    public FeatureFinderDefs
  {
public:
    typedef FeaFiModule<PeakType, FeatureType> Base;

    typedef typename Base::IntensityType IntensityType;
    typedef typename Base::CoordinateType CoordinateType;
    typedef DoubleReal ProbabilityType;

    SimpleExtender(const MSExperiment<PeakType> * map, FeatureMap<FeatureType> * features, FeatureFinder * ff) :
      Base(map, features, ff),
      last_pos_extracted_()
    {
      this->setName("SimpleExtender");

      this->defaults_.setValue("dist_mz_up", 6.0, "Maximum high m/z distance of peak in the region/boundary from the seed.");
      this->defaults_.setMinFloat("dist_mz_up", 0.0);
      this->defaults_.setValue("dist_mz_down", 2.0, "Maximum low m/z distance of peak in the region/boundary from the seed.");
      this->defaults_.setMinFloat("dist_mz_down", 0.0);
      this->defaults_.setValue("dist_rt_up", 5.0, "Maximum high RT distance of peak in the region/boundary from the seed.");
      this->defaults_.setMinFloat("dist_rt_up", 0.0);
      this->defaults_.setValue("dist_rt_down", 5.0, "Maximum low RT distance of peak in the region/boundary from the seed.");
      this->defaults_.setMinFloat("dist_rt_down", 0.0);

      // priority check is per default switched off
      // these values were used for the Myoglobin quantification project
      // DON'T REMOVE THIS
      this->defaults_.setValue("priority_thr", -0.1, "Minimum priority for data points to be included into the boundary of the feature (default 0.0). The priority of a data point is a function of its intensity and its distance to the last point included into the feature region. Setting this threshold to zero or a very small value is usually a good idea.", StringList::create("advanced"));

      this->defaults_.setValue("intensity_factor", 0.03, "Influences for intensity (ion count) threshold in the feature extension. We include only raw data points into this region if their intensity is larger than [intensity_factor * (intensity of the seed)].");
      this->defaults_.setMinFloat("intensity_factor", 0.0);
      this->defaults_.setMaxFloat("intensity_factor", 1.0);

      this->defaultsToParam_();
    }

    virtual ~SimpleExtender()
    {
    }

    void extend(const ChargedIndexSet & seed_region, ChargedIndexSet & result_region)
    {
      // empty region and boundary datastructures
      result_region.clear();
      priorities_.clear();
      running_avg_.clear();
      boundary_ = std::priority_queue<IndexWithPriority, std::vector<IndexWithPriority>, typename IndexWithPriority::PriorityLess>();

#ifdef DEBUG_FEATUREFINDER
      std::vector<IndexPair> debug_vector;
#endif
      // find maximum of region (seed)
      CoordinateType max_intensity = 0.0;
      IndexPair seed;

      for (IndexSet::const_iterator citer = seed_region.begin(); citer != seed_region.end(); ++citer)
      {
        if (this->getPeakIntensity(*citer) > max_intensity)
        {
          seed = *citer;
          max_intensity = this->getPeakIntensity(seed);
        }
      }

      // remember last extracted point (in this case the seed !)
      last_pos_extracted_[Peak2D::RT] = this->getPeakRt(seed);
      last_pos_extracted_[Peak2D::MZ] = this->getPeakMz(seed);

      // Add peaks received from seeder directly to boundary
      for (IndexSet::const_iterator citer = seed_region.begin(); citer != seed_region.end(); ++citer)
      {
        ProbabilityType priority = computePeakPriority_(*citer);
        priorities_[*citer] = priority;
        boundary_.push(IndexWithPriority(*citer, priority));
      }
      // pass on charge information
      result_region.charge = seed_region.charge;

      // re-compute intensity threshold
      intensity_threshold_ = (DoubleReal)(this->param_).getValue("intensity_factor") * this->getPeakIntensity(seed);

#ifdef DEBUG_FEATUREFINDER
      std::cout << "\n";
      std::cout << "Extending from " << this->getPeakRt(seed) << "/" << this->getPeakMz(seed) << std::endl;
      std::cout << "Intensity of seed " << this->getPeakIntensity(seed);
      std::cout << " (" << seed.first << "/" << seed.second << ")" << std::endl;
      std::cout << "Intensity_threshold: " << intensity_threshold_ << std::endl;
#endif

      while (!boundary_.empty())
      {
        // remove peak with highest priority
        const IndexPair  current_index = boundary_.top().index;
        boundary_.pop();

        //  check for corrupt index
        OPENMS_PRECONDITION(current_index.first < (*this->map_).size(), "Scan index outside of map!");
        OPENMS_PRECONDITION(current_index.second < (*this->map_)[current_index.first].size(), "Peak index outside of scan!");

        // remember last extracted peak
        last_pos_extracted_[Peak2D::RT] = this->getPeakRt(current_index);
        last_pos_extracted_[Peak2D::MZ] = this->getPeakMz(current_index);

        // Now we explore the neighbourhood of the current peak. Points in this area are included
        // into the boundary if their intensity is not too low and they are not too
        // far away from the seed.
        // Add position to the current average of positions weighted by intensity
        running_avg_.add(last_pos_extracted_, this->getPeakIntensity(current_index));

        // explore neighbourhood of current peak
        moveMzUp_(current_index);
        moveMzDown_(current_index);
        moveRtUp_(current_index);
        moveRtDown_(current_index);

        // set peak flags and add to boundary
        this->ff_->getPeakFlag(current_index) = USED;
#ifdef DEBUG_FEATUREFINDER
        debug_vector.push_back(current_index);
#endif
        result_region.insert(current_index);

      }       // end of while ( !boundary_.empty() )

#ifdef DEBUG_FEATUREFINDER
      std::cout << "Feature region size: " << result_region.size() << std::endl;
#endif

#ifdef DEBUG_FEATUREFINDER
      static UInt number = 1;
      writeDebugFile_(debug_vector, number++);
      debug_vector.clear();
#endif

      return;
    }     // end of extend

    struct IndexWithPriority
    {
      IndexWithPriority(const FeatureFinderDefs::IndexPair & i, DoubleReal p) :
        index(i),
        priority(p)
      {
      }

      IndexPair index;
      ProbabilityType priority;

      struct PriorityLess
      {
        inline bool operator()(const IndexWithPriority & x, const IndexWithPriority & y) const
        {
          return x.priority < y.priority;
        }

      };
    };

protected:

    virtual void updateMembers_()
    {
      dist_mz_up_ = this->param_.getValue("dist_mz_up");
      dist_mz_down_ = this->param_.getValue("dist_mz_down");
      dist_rt_up_ = this->param_.getValue("dist_rt_up");
      dist_rt_down_ = this->param_.getValue("dist_rt_down");
      priority_threshold_ = this->param_.getValue("priority_thr");
    }

    void writeDebugFile_(const std::vector<IndexPair> & peaks, UInt nr_feat)
    {
      String filename = String(nr_feat).fillLeft('0', 4) + "_Extension.dta2d";
      std::ofstream file(filename.c_str());
      for (Size i = 0; i < peaks.size(); ++i)
      {
        file << this->getPeakRt(peaks[i]) << " " << this->getPeakMz(peaks[i]) << " " << peaks.size() - i << std::endl;
      }
      file.close();
    }

    bool isTooFarFromCentroid_(const IndexPair & index)
    {
      //Corrupt index
      OPENMS_PRECONDITION(index.first < (*this->map_).size(), "Scan index outside of map!");
      OPENMS_PRECONDITION(index.second < (*this->map_)[index.first].size(), "Peak index outside of scan!");

      const DPosition<2> & curr_mean = running_avg_.getPosition();

      if (this->getPeakMz(index) > curr_mean[Peak2D::MZ] + dist_mz_up_   ||
          this->getPeakMz(index) < curr_mean[Peak2D::MZ] - dist_mz_down_ ||
          this->getPeakRt(index) > curr_mean[Peak2D::RT] + dist_rt_up_   ||
          this->getPeakRt(index) < curr_mean[Peak2D::RT] - dist_rt_down_)
      {
        //too far
        return true;
      }

      //close enough
      return false;
    }

    void moveMzUp_(const IndexPair & index)
    {
      try
      {
        IndexPair tmp = index;
        while (true)
        {
          this->getNextMz(tmp);
          if (isTooFarFromCentroid_(tmp)) break;
          checkNeighbour_(tmp);
        }
      }
      catch (NoSuccessor)
      {
      }
    }

    void moveMzDown_(const IndexPair & index)
    {
      try
      {
        IndexPair tmp = index;
        while (true)
        {
          this->getPrevMz(tmp);
          if (isTooFarFromCentroid_(tmp)) break;
          checkNeighbour_(tmp);
        }
      }
      catch (NoSuccessor)
      {
      }
    }

    void moveRtUp_(const IndexPair & index)
    {
      try
      {
        IndexPair tmp = index;

        while (true)
        {
          this->getNextRt(tmp);
          if (isTooFarFromCentroid_(tmp)) break;
          checkNeighbour_(tmp);
        }
      }
      catch (NoSuccessor)
      {
      }
    }

    void moveRtDown_(const IndexPair & index)
    {
      try
      {
        IndexPair tmp = index;
        while (true)
        {
          this->getPrevRt(tmp);
          if (isTooFarFromCentroid_(tmp)) break;
          checkNeighbour_(tmp);
        }
      }
      catch (NoSuccessor)
      {
      }
    }

    ProbabilityType computePeakPriority_(const IndexPair & index)
    {
      return (*this->map_)[index.first][index.second].getIntensity();
    }

    void checkNeighbour_(const IndexPair & index)
    {
      //Corrupt index
      OPENMS_PRECONDITION(index.first < (*this->map_).size(), "Scan index outside of map!");
      OPENMS_PRECONDITION(index.second < (*this->map_)[index.first].size(), "Peak index outside of scan!");

      // skip this point if its intensity is too low
      if (this->getPeakIntensity(index) <= intensity_threshold_)
      {
        return;
      }
      if (this->ff_->getPeakFlag(index) == UNUSED)
      {
        DoubleReal pr_new = computePeakPriority_(index);

        if (pr_new > priority_threshold_)
        {
          //std::map<IndexPair, DoubleReal>::iterator piter = priorities_.find(index);
          this->ff_->getPeakFlag(index) = USED;
          priorities_[index] = pr_new;
          boundary_.push(IndexWithPriority(index, pr_new));
        }
      }
    }

    Math::AveragePosition<2> running_avg_;

    std::map<IndexPair, ProbabilityType> priorities_;

    DPosition<2> last_pos_extracted_;

    std::priority_queue<IndexWithPriority, std::vector<IndexWithPriority>, typename IndexWithPriority::PriorityLess> boundary_;

    IntensityType intensity_threshold_;

    CoordinateType dist_mz_up_;
    CoordinateType dist_mz_down_;
    CoordinateType dist_rt_up_;
    CoordinateType dist_rt_down_;

    ProbabilityType priority_threshold_;

    ChargedIndexSet region_;

private:
    SimpleExtender();
    SimpleExtender & operator=(const SimpleExtender &);
    SimpleExtender(const SimpleExtender &);

  };
}
#endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_SIMPLEEXTENDER_H