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

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

fitting a Local with no data is impossible so added a throw in that case

  • 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 1472 2008-09-03 16:46:34Z 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    if (data_.size()<step_size){
70      std::stringstream ss;
71      ss << "yat::regression::Local: too large step_size "
72         << "step_size, " << step_size
73         << ", is larger than number of added data points " << data_.size();
74      throw std::runtime_error(ss.str());
75    }
76
77    size_t nof_fits=data_.size()/step_size;
78    x_.resize(nof_fits);
79    y_predicted_.resize(x_.size());
80    y_err_.resize(x_.size());
81    sort(data_.begin(), data_.end());
82
83    // copying data to 2 utility vectors ONCE to use views from
84    utility::Vector x(data_.size());
85    utility::Vector y(data_.size());
86    for (size_t j=0; j<x.size(); j++){
87      x(j)=data_[j].first;
88      y(j)=data_[j].second;
89    }
90
91    // looping over regression points and perform local regression
92    for (size_t i=0; i<nof_fits; i++) {
93      size_t max_index = static_cast<size_t>( (i+0.5)*step_size );
94      size_t min_index;
95      double width; // distance from middle of windo to border of window
96      double x_mid; // middle of window
97      // right border case
98      if (max_index > data_.size()-1){
99        min_index = max_index - nof_points + 1;
100        max_index = data_.size()-1;
101        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
102                  ( 2*(max_index-min_index)) );
103        x_mid = x(min_index)+width;
104      }
105      // normal middle case
106      else if (max_index > nof_points-1){
107        min_index = max_index - nof_points + 1;
108        width = (x(max_index)-x(min_index))/2;
109        x_mid = x(min_index)+width;
110      }
111      // left border case
112      else {
113        min_index = 0;
114        width = ( (( x(max_index)-x(0) )*(nof_points-1)) / 
115                  ( 2*(max_index-min_index)) );
116        x_mid = x(max_index)-width;
117      }
118      assert(min_index<data_.size());
119      assert(max_index<data_.size());
120                               
121      utility::VectorView x_local(x, min_index, max_index-min_index+1);
122      utility::VectorView y_local(y, min_index, max_index-min_index+1);
123
124      // calculating weights
125      utility::Vector w(max_index-min_index+1);
126      for (size_t j=0; j<w.size(); j++)
127        w(j) = (*kernel_)( (x_local(j)- x_mid)/width );
128     
129      // fitting the regressor locally
130      regressor_->fit(x_local,y_local,w);
131      assert(i<y_predicted_.size());
132      assert(i<y_err_.size());
133      x_(i) = x(i*step_size);
134      y_predicted_(i) = regressor_->predict(x(i*step_size));
135      y_err_(i) = sqrt(regressor_->standard_error2(x(i*step_size)));
136    }
137  }
138
139
140  void Local::reset(void)
141  {
142    data_.clear();
143    x_.resize(0);
144    y_predicted_.resize(0);
145    y_err_.resize(0);
146  }
147
148
149  const utility::Vector& Local::x(void) const
150  {
151    return x_;
152  }
153
154  const utility::Vector& Local::y_predicted(void) const
155  {
156    return y_predicted_;
157  }
158
159  const utility::Vector& Local::y_err(void) const
160  {
161    return y_err_;
162  }
163
164  std::ostream& operator<<(std::ostream& os, const Local& r)
165  {
166    os << "# column 1: x\n"
167      << "# column 2: y\n"
168      << "# column 3: y_err\n";
169    for (size_t i=0; i<r.x().size(); i++) {
170      os << r.x()(i) << "\t" 
171         << r.y_predicted()(i) << "\t"
172         << r.y_err()(i) << "\n";
173    }   
174
175    return os;
176  }
177
178}}} // of namespaces regression, yat, and theplu
Note: See TracBrowser for help on using the repository browser.