source: trunk/yat/regression/Local.cc @ 1009

Last change on this file since 1009 was 1009, checked in by Peter, 13 years ago

merging branch peter-dev into trunk delta 1008:994

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.1 KB
Line 
1// $Id: Local.cc 1009 2008-02-01 04:15:44Z peter $
2
3/*
4  Copyright (C) 2004 Peter Johansson
5  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
6
7  This file is part of the yat library, http://trac.thep.lu.se/yat
8
9  The yat library is free software; you can redistribute it and/or
10  modify it under the terms of the GNU General Public License as
11  published by the Free Software Foundation; either version 2 of the
12  License, or (at your option) any later version.
13
14  The yat library is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  General Public License for more details.
18
19  You should have received a copy of the GNU General Public License
20  along with this program; if not, write to the Free Software
21  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
22  02111-1307, USA.
23*/
24
25#include "Local.h"
26#include "Kernel.h"
27#include "OneDimensionalWeighted.h"
28#include "yat/utility/vector.h"
29#include "yat/utility/vectorView.h"
30
31#include <algorithm>
32#include <cassert>
33#include <iostream>
34
35namespace theplu {
36namespace yat {
37namespace regression {
38
39  Local::Local(OneDimensionalWeighted& r, Kernel& k)
40    : kernel_(&k), regressor_(&r)
41  {
42  }
43
44  Local::~Local(void)
45  {
46  }
47
48  void Local::add(const double x, const double y)
49  {
50    data_.push_back(std::make_pair(x,y));
51  }
52
53  void Local::fit(const size_t step_size, const size_t nof_points)
54  {
55    if (step_size==0 || nof_points<3){
56      std::cerr << "yat::regression::Local "
57                << "Parameters invalid. Fitting ignored." << std::endl;
58      return;
59    }
60
61    size_t nof_fits=data_.size()/step_size;
62    x_ = utility::vector(nof_fits);
63    y_predicted_ = utility::vector(x_.size());
64    y_err_ = utility::vector(x_.size());
65    sort(data_.begin(), data_.end());
66
67    // coying data to 2 utility vectors ONCE to use views from
68    utility::vector x(data_.size());
69    utility::vector y(data_.size());
70    for (size_t j=0; j<x.size(); j++){
71      x(j)=data_[j].first;
72      y(j)=data_[j].second;
73    }
74
75    // looping over regression points and perform local regression
76    for (size_t i=0; i<nof_fits; i++) {
77      size_t max_index = static_cast<size_t>( (i+0.5)*step_size );
78      size_t min_index;
79      double width; // distance from middle of windo to border of window
80      double x_mid; // middle of window
81      // right border case
82      if (max_index > data_.size()-1){
83        min_index = max_index - nof_points + 1;
84        max_index = data_.size()-1;
85        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
86                  ( 2*(max_index-min_index)) );
87        x_mid = x(min_index)+width;
88      }
89      // normal middle case
90      else if (max_index > nof_points-1){
91        min_index = max_index - nof_points + 1;
92        width = (x(max_index)-x(min_index))/2;
93        x_mid = x(min_index)+width;
94      }
95      // left border case
96      else {
97        min_index = 0;
98        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
99                  ( 2*(max_index-min_index)) );
100        x_mid = x(max_index)-width;
101      }
102      assert(min_index<data_.size());
103      assert(max_index<data_.size());
104                               
105      utility::vectorView x_local(x, min_index, max_index-min_index+1);
106      utility::vectorView y_local(y, min_index, max_index-min_index+1);
107
108      // calculating weights
109      utility::vector w(max_index-min_index+1);
110      for (size_t j=0; j<w.size(); j++)
111        w(j) = (*kernel_)( (x_local(j)- x_mid)/width );
112     
113      // fitting the regressor locally
114      regressor_->fit(x_local,y_local,w);
115      assert(i<y_predicted_.size());
116      assert(i<y_err_.size());
117      y_predicted_(i) = regressor_->predict(x(i*step_size));
118      y_err_(i) = sqrt(regressor_->standard_error2(x(i*step_size)));
119    }
120  }
121
122  const utility::vector& Local::x(void) const
123  {
124    return x_;
125  }
126
127  const utility::vector& Local::y_predicted(void) const
128  {
129    return y_predicted_;
130  }
131
132  const utility::vector& Local::y_err(void) const
133  {
134    return y_err_;
135  }
136
137  std::ostream& operator<<(std::ostream& os, const Local& r)
138  {
139    os << "# column 1: x\n"
140      << "# column 2: y\n"
141      << "# column 3: y_err\n";
142    for (size_t i=0; i<r.x().size(); i++) {
143      os << r.x()(i) << "\t" 
144         << r.y_predicted()(i) << "\t"
145         << r.y_err()(i) << "\n";
146    }   
147
148    return os;
149  }
150
151}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.