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

Last change on this file since 1655 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.6 KB
Line 
1// $Id: Local.cc 1487 2008-09-10 08:41:36Z jari $
2
3/*
4  Copyright (C) 2004 Peter Johansson
5  Copyright (C) 2005, 2006, 2007 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2008 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 "Local.h"
25#include "Kernel.h"
26#include "OneDimensionalWeighted.h"
27#include "yat/utility/Vector.h"
28#include "yat/utility/VectorView.h"
29
30#include <algorithm>
31#include <cassert>
32#include <iostream>
33#include <sstream>
34#include <stdexcept>
35
36namespace theplu {
37namespace yat {
38namespace regression {
39
40  Local::Local(OneDimensionalWeighted& r, Kernel& k)
41    : kernel_(&k), regressor_(&r)
42  {
43  }
44
45  Local::~Local(void)
46  {
47  }
48
49  void Local::add(const double x, const double y)
50  {
51    data_.push_back(std::make_pair(x,y));
52  }
53
54  void Local::fit(const size_t step_size, const size_t nof_points)
55  {
56    if (step_size==0){
57      std::stringstream ss;
58      ss << "yat::regression::Local: step_size must be larger than zero.";
59      throw std::runtime_error(ss.str());
60    }
61    if (nof_points<3){
62      std::stringstream ss;
63      ss << "yat::regression::Local: too few data points. "
64         << "At least 3 data points are needed to perform fitting.";
65      throw std::runtime_error(ss.str());
66    }
67    if (data_.size()<step_size){
68      std::stringstream ss;
69      ss << "yat::regression::Local: too large step_size "
70         << "step_size, " << step_size
71         << ", is larger than number of added data points " << data_.size();
72      throw std::runtime_error(ss.str());
73    }
74
75    size_t nof_fits=data_.size()/step_size;
76    x_.resize(nof_fits);
77    y_predicted_.resize(x_.size());
78    y_err_.resize(x_.size());
79    sort(data_.begin(), data_.end());
80
81    // copying data to 2 utility vectors ONCE to use views from
82    utility::Vector x(data_.size());
83    utility::Vector y(data_.size());
84    for (size_t j=0; j<x.size(); j++){
85      x(j)=data_[j].first;
86      y(j)=data_[j].second;
87    }
88
89    // looping over regression points and perform local regression
90    for (size_t i=0; i<nof_fits; i++) {
91      size_t max_index = static_cast<size_t>( (i+0.5)*step_size );
92      size_t min_index;
93      double width; // distance from middle of windo to border of window
94      double x_mid; // middle of window
95      // right border case
96      if (max_index > data_.size()-1){
97        min_index = max_index - nof_points + 1;
98        max_index = data_.size()-1;
99        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
100                  ( 2*(max_index-min_index)) );
101        x_mid = x(min_index)+width;
102      }
103      // normal middle case
104      else if (max_index > nof_points-1){
105        min_index = max_index - nof_points + 1;
106        width = (x(max_index)-x(min_index))/2;
107        x_mid = x(min_index)+width;
108      }
109      // left border case
110      else {
111        min_index = 0;
112        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
113                  ( 2*(max_index-min_index)) );
114        x_mid = x(max_index)-width;
115      }
116      assert(min_index<data_.size());
117      assert(max_index<data_.size());
118                               
119      utility::VectorView x_local(x, min_index, max_index-min_index+1);
120      utility::VectorView y_local(y, min_index, max_index-min_index+1);
121
122      // calculating weights
123      utility::Vector w(max_index-min_index+1);
124      for (size_t j=0; j<w.size(); j++)
125        w(j) = (*kernel_)( (x_local(j)- x_mid)/width );
126     
127      // fitting the regressor locally
128      regressor_->fit(x_local,y_local,w);
129      assert(i<y_predicted_.size());
130      assert(i<y_err_.size());
131      x_(i) = x(i*step_size);
132      y_predicted_(i) = regressor_->predict(x(i*step_size));
133      y_err_(i) = sqrt(regressor_->standard_error2(x(i*step_size)));
134    }
135  }
136
137
138  void Local::reset(void)
139  {
140    data_.clear();
141    x_.resize(0);
142    y_predicted_.resize(0);
143    y_err_.resize(0);
144  }
145
146
147  const utility::Vector& Local::x(void) const
148  {
149    return x_;
150  }
151
152  const utility::Vector& Local::y_predicted(void) const
153  {
154    return y_predicted_;
155  }
156
157  const utility::Vector& Local::y_err(void) const
158  {
159    return y_err_;
160  }
161
162  std::ostream& operator<<(std::ostream& os, const Local& r)
163  {
164    os << "# column 1: x\n"
165      << "# column 2: y\n"
166      << "# column 3: y_err\n";
167    for (size_t i=0; i<r.x().size(); i++) {
168      os << r.x()(i) << "\t" 
169         << r.y_predicted()(i) << "\t"
170         << r.y_err()(i) << "\n";
171    }   
172
173    return os;
174  }
175
176}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.