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

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

adding reset function in regression::Local

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