source: trunk/src/RegressionLocal.cc @ 221

// $Id: RegressionLocal.cc 221 2004-12-30 22:36:25Z peter $

#include "RegressionLocal.h"

#include "Regression.h"
#include "RegressionKernel.h"
#include "RegressionLinear.h"
#include "vector.h"

#include "algorithm"

namespace theplu {
namespace statistics {

  RegressionLocal::RegressionLocal(const gslapi::vector& x, 
                                   const gslapi::vector& y, 
                                   Regression& r, 
                                   RegressionKernel& k,
                                   const size_t nof_predictions)
    : kernel_(&k), regression_(&r)
  {
    estimated_x_ = gslapi::vector(nof_predictions);
    estimated_y_ = estimated_x_;
    estimated_y_err_ = estimated_x_;
    //sorting data with respect to x
    // Peter, could this be done without all the copying
    std::vector<std::pair<double, double> > data;
    for (size_t i=0; i<x.size(); i++){
      std::pair<double, double> tmp(x(i),y(i));
      data_.push_back(tmp);
    }
    sort(data_.begin(), data_.end());
    for (size_t i=0; i<nof_predictions; i++)
      estimated_x_(i) =
        data_[static_cast<size_t>((i+0.5)/nof_predictions*data_.size())].first;
  }

  void RegressionLocal::fit(const double f)
  {
    // number of points on each side
    size_t nof_data_points = static_cast<size_t>(f/2*data_.size());
    size_t index_min = 0;
    size_t index_max = 0;
    for (size_t i=0; i<estimated_x_.size(); i++) {
      // extreme left case
      if (data_[nof_data_points-1].first>estimated_x_(i)){
        index_min = 0;
        index_max = 2*nof_data_points-1;
      }
      
      // extreme right case
      else if (data_[data_.size()-estimated_x_.size()-1].first<
               estimated_x_(i)){
        index_min = data_.size()-2*nof_data_points;
        index_max = data_.size()-1;
      }
      
      // normal case - window in the middle
      else {
        for (size_t j=nof_data_points; j<data_.size(); j++)
          if (data_[j].first>estimated_x_(i)) {
            index_min = j-nof_data_points;
            index_max = j+nof_data_points-1;
            break;
          }
      }
      
      // copying data 
      // Peter, this is a stupid solution which could be
      // improved by using a moving view (in the vector class).
      gslapi::vector x(index_max-index_min+1);
      gslapi::vector y(index_max-index_min+1);
      gslapi::vector w(index_max-index_min+1);
      for (size_t j=0; j<x.size(); j++)
        x(j)=data_[index_min+j].first;
      for (size_t j=0; j<y.size(); j++){
        y(j)=data_[index_min+j].second;
        if (y(j)>10)
          std::cout << x(j) << "\t" << y(j) << std::endl;
      }
      // calculating weights
      for (size_t j=0; j<y.size(); j++)
        w(j) = kernel_->weight( 2*(x(j)-data_[index_min].first)/
                                (data_[index_max].first-
                                 data_[index_min].first)-1);

      
      // fitting the regressor locally
      std::cout << "x: " << estimated_x_(i) << std::endl;
      std::cout << "x.begin: " << x(0) << std::endl;
      std::cout << "x.end: " << x(x.size()-1) << std::endl;
      regression_->fit(x,y,w);
      std::cout << "estimating value\n";
      regression_->predict(estimated_x_(i),
                            estimated_y_(i),
                            estimated_y_err_(i));
    }
  }
  
  std::ostream& operator<< (std::ostream& s, const RegressionLocal& r)
  {
    gslapi::vector x=r.estimated_x();
    gslapi::vector y(r.estimated_y());
    gslapi::vector y_err(r.estimated_y_err());
    for (size_t i=0; i<x.size(); i++){
      s << x(i) << "\t";
      s << y(i) << "\t";
      s << y_err(i) << "\n";
    }
    return s;
  }

}} // of namespace cpptools and namespace theplu