Changeset 216 for trunk/src/RegressionLocal.cc

Timestamp:

Dec 29, 2004, 10:29:54 AM (18 years ago)

Author:

Peter

Message:

RegressionLocal? added

File:

• trunk/src/RegressionLocal.cc (modified) (2 diffs)

trunk/src/RegressionLocal.cc

@@ -8 +8 @@
 #include "vector.h"
 
+#include "algorithm"
 
 namespace theplu {
@@ -14 +15 @@
   RegressionLocal::RegressionLocal(const gslapi::vector& x, 
                                    const gslapi::vector& y, 
-                                   const Regression& r, 
-                                   const RegressionKernel& k)
-    : kernel_(&k), regression_(&r)
+                                   Regression& r, 
+                                   RegressionKernel& k,
+                                   const gslapi::vector& estimated_x)
+    : kernel_(&k), regression_(&r), estimated_x_(estimated_x), 
+      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());
+  }
+
+  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 when the vector class is updated.
+      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;
+      // 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
+      regression_->fit_weighted(x,y,w);
+      regression_->estimate(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