yat/regression/Local.cc

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

Generated on Sat Oct 14 02:23:16 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)