source: trunk/yat/regression/Local.cc

Last change on this file was 2919, checked in by Peter, 9 years ago

update copyright years

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.7 KB
Line 
1// $Id: Local.cc 2919 2012-12-19 06:54:23Z peter $
2
3/*
4  Copyright (C) 2004 Peter Johansson
5  Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2010, 2011, 2012 Peter Johansson
7
8  This file is part of the yat library, http://dev.thep.lu.se/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 3 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with yat. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#include <config.h>
25
26#include "Local.h"
27#include "Kernel.h"
28#include "OneDimensionalWeighted.h"
29
30#include "yat/utility/Exception.h"
31#include "yat/utility/Vector.h"
32#include "yat/utility/VectorView.h"
33
34#include <algorithm>
35#include <cassert>
36#include <ostream>
37#include <sstream>
38#include <stdexcept>
39
40namespace theplu {
41namespace yat {
42namespace regression {
43
44  Local::Local(OneDimensionalWeighted& r, Kernel& k)
45    : kernel_(&k), regressor_(&r)
46  {
47  }
48
49  Local::~Local(void)
50  {
51  }
52
53  void Local::add(const double x, const double y)
54  {
55    data_.push_back(std::make_pair(x,y));
56  }
57
58  void Local::fit(const size_t step_size, const size_t nof_points)
59  {
60    if (step_size==0){
61      std::stringstream ss;
62      ss << "yat::regression::Local: step_size must be larger than zero.";
63      throw utility::runtime_error(ss.str());
64    }
65    if (nof_points<3){
66      std::stringstream ss;
67      ss << "yat::regression::Local: too few data points. "
68         << "At least 3 data points are needed to perform fitting.";
69      throw utility::runtime_error(ss.str());
70    }
71    if (data_.size()<step_size){
72      std::stringstream ss;
73      ss << "yat::regression::Local: too large step_size "
74         << "step_size, " << step_size
75         << ", is larger than number of added data points " << data_.size();
76      throw utility::runtime_error(ss.str());
77    }
78
79    size_t nof_fits=data_.size()/step_size;
80    x_.resize(nof_fits);
81    y_predicted_.resize(x_.size());
82    y_err_.resize(x_.size());
83    sort(data_.begin(), data_.end());
84
85    // copying data to 2 utility vectors ONCE to use views from
86    utility::Vector x(data_.size());
87    utility::Vector y(data_.size());
88    for (size_t j=0; j<x.size(); j++){
89      x(j)=data_[j].first;
90      y(j)=data_[j].second;
91    }
92
93    // looping over regression points and perform local regression
94    for (size_t i=0; i<nof_fits; i++) {
95      size_t max_index = static_cast<size_t>( (i+0.5)*step_size );
96      size_t min_index;
97      double width; // distance from middle of windo to border of window
98      double x_mid; // middle of window
99      // right border case
100      if (max_index > data_.size()-1){
101        min_index = max_index - nof_points + 1;
102        max_index = data_.size()-1;
103        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
104                  ( 2*(max_index-min_index)) );
105        x_mid = x(min_index)+width;
106      }
107      // normal middle case
108      else if (max_index > nof_points-1){
109        min_index = max_index - nof_points + 1;
110        width = (x(max_index)-x(min_index))/2;
111        x_mid = x(min_index)+width;
112      }
113      // left border case
114      else {
115        min_index = 0;
116        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
117                  ( 2*(max_index-min_index)) );
118        x_mid = x(max_index)-width;
119      }
120      assert(min_index<data_.size());
121      assert(max_index<data_.size());
122                               
123      utility::VectorView x_local(x, min_index, max_index-min_index+1);
124      utility::VectorView y_local(y, min_index, max_index-min_index+1);
125
126      // calculating weights
127      utility::Vector w(max_index-min_index+1);
128      for (size_t j=0; j<w.size(); j++)
129        w(j) = (*kernel_)( (x_local(j)- x_mid)/width );
130     
131      // fitting the regressor locally
132      regressor_->fit(x_local,y_local,w);
133      assert(i<y_predicted_.size());
134      assert(i<y_err_.size());
135      x_(i) = x(i*step_size);
136      y_predicted_(i) = regressor_->predict(x(i*step_size));
137      y_err_(i) = sqrt(regressor_->standard_error2(x(i*step_size)));
138    }
139  }
140
141
142  void Local::reset(void)
143  {
144    data_.clear();
145    x_.resize(0);
146    y_predicted_.resize(0);
147    y_err_.resize(0);
148  }
149
150
151  const utility::Vector& Local::x(void) const
152  {
153    return x_;
154  }
155
156  const utility::Vector& Local::y_predicted(void) const
157  {
158    return y_predicted_;
159  }
160
161  const utility::Vector& Local::y_err(void) const
162  {
163    return y_err_;
164  }
165
166  std::ostream& operator<<(std::ostream& os, const Local& r)
167  {
168    os << "# column 1: x\n"
169      << "# column 2: y\n"
170      << "# column 3: y_err\n";
171    for (size_t i=0; i<r.x().size(); i++) {
172      os << r.x()(i) << "\t" 
173         << r.y_predicted()(i) << "\t"
174         << r.y_err()(i) << "\n";
175    }   
176
177    return os;
178  }
179
180}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.