is_picked_systematic.hpp

Go to the documentation of this file.

// (C) Copyright Renaud Detry   2007-2011.
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
 
#ifndef TRSL_IS_PICKED_SYSTEMATIC_HPP
#define TRSL_IS_PICKED_SYSTEMATIC_HPP
 
#include <trsl/common.hpp>
#include <trsl/weight_accessor.hpp>
 
#include <algorithm>
#include <functional>
#include <utility>
#include <cstdlib> // rand & random
#include <limits>
#include <cassert>
#include <boost/static_assert.hpp>
 
namespace trsl {
    
  template<
    typename ElementType,
    typename WeightType = double,
    typename WeightAccessor = mp_weight_accessor<WeightType, ElementType>
  > class is_picked_systematic
  {
  private:
    BOOST_STATIC_ASSERT((std::numeric_limits<WeightType>::is_integer == false));
  public:
    typedef ElementType element_type;
    typedef WeightType weight_type;
    typedef WeightAccessor weight_accessor_type;
    
    is_picked_systematic() :
      sampleSize_(0),
      populationWeight_(0)
      {
        initialize( 0 );
      }
    
    is_picked_systematic(size_t sampleSize,
                         WeightType populationWeight,
                         WeightAccessor const& wac = WeightAccessor()) :
      wac_(wac), sampleSize_(sampleSize),
      populationWeight_(populationWeight)
      {
        initialize( rand_gen::uniform_01<WeightType>() );
      }
 
    is_picked_systematic(size_t sampleSize,
                         WeightType populationWeight,
                         WeightType uniform01,
                         WeightAccessor const& wac = WeightAccessor()) :
      wac_(wac), sampleSize_(sampleSize),
      populationWeight_(populationWeight)
      {
        initialize(uniform01);
      }
 
    bool operator()(const ElementType & e)
      {
        if (sampleSize_ == 0) return false;
        
#ifdef TRSL_USE_SYSTEMATIC_INTUITIVE_ALGORITHM
        // This algorithm is the intuitive implementation of
        // systematic sampling, where one pictures a wheel with
        // N_SAMPLE spokes and distributes the population around the
        // wheel as segments of tire of length proportional to their
        // weight; the spokes point to picked elements.
        WeightType arrow = k_*step_;
        assert(cumulative_ <= arrow);
        if (arrow < cumulative_ + wac_(e))
        {
          k_++;
          return true;
        }
        cumulative_ += wac_(e);
        return false;
#else
        // This algorithm is a massaged version of the intuitive
        // algorithm.  Both algorithms are conceptually identical, but
        // this version is faster.
        assert(position_ >= 0);
        if (position_ < wac_(e))
        {
          position_ += step_;
          return true;
        }
        position_ -= wac_(e);
        return false;
#endif
      }
 
    bool is_first_pick(const ElementType & e) const
    {
      if (wac_(e) <= step_) return true;
#ifdef TRSL_USE_SYSTEMATIC_INTUITIVE_ALGORITHM
      return k_*step_ - cumulative_ < 2*step_;
#else
      return position_ < 2*step_;
#endif      
    }
 
    bool operator== (const is_picked_systematic<ElementType, WeightType, WeightAccessor> &p) const
      {
        if (sampleSize_ != p.sampleSize_ ||
            populationWeight_ != p.populationWeight_) return false;
 
#ifdef TRSL_USE_SYSTEMATIC_INTUITIVE_ALGORITHM
        if (cumulative_ != p.cumulative_ ||
            k_ != p.k_) return false;
#else
        if (position_ != p.position_) return false;
#endif
        
        return true;
      }
      
  private:
    void initialize(WeightType randomReal)
      {
        // If one could rely on
        //
        //   (x!=0.0)/0.0 == inf (with the sign of x)
        //
        // and on
        //
        //   0.0 * inf == NaN,
        //
        // neither of these throwing exceptions, then the following
        // block could be avoided, and the case sampleSize_ == 0 would
        // be implicitly handled by inf/NaN arithmetic.
        //
        // However, it would be necessary to maintain two different
        // implementations depending on
        // std::numeric_limits<double>::is_iec559, and even then I'm
        // not sure if inf/NaN arithmetic is reliable.
        if (sampleSize_ != 0)
          step_ = populationWeight_ / sampleSize_;
        else
          // A semantically correct value for step_ would be inf, but
          // we don't want to get into has_infinity trouble...
          // When sampleSize_==0, step_ should be considered undefined.
          step_ = 0;
        
        WeightType randomWeight = randomReal * step_;
#ifdef TRSL_USE_SYSTEMATIC_INTUITIVE_ALGORITHM
        cumulative_ = -randomWeight;
        k_ = 0;
#else
        position_ = randomWeight;
#endif
      }
 
  private:
    WeightAccessor wac_;
    size_t sampleSize_;
    WeightType populationWeight_;
    WeightType step_;
  
#ifdef TRSL_USE_SYSTEMATIC_INTUITIVE_ALGORITHM
    WeightType cumulative_;
    size_t k_;
#else
    WeightType position_;
#endif
  };
 
  template<class Iterator>
  bool is_first_pick(const Iterator& i)
  {
    return i.predicate().is_first_pick(*i);
  }
  
}
 
#endif // include guard