yat/regression/GSLInterpolation.cc

Code
Comments
Other
Rev Date Author Line
193 27 Oct 04 peter 1 // $Id$
193 27 Oct 04 peter 2
675 10 Oct 06 jari 3 /*
4359 23 Aug 23 peter 4   Copyright (C) 2004, 2005, 2006 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 5   Copyright (C) 2007 Peter Johansson
4359 23 Aug 23 peter 6   Copyright (C) 2008 Jari Häkkinen, Peter Johansson
2119 12 Dec 09 peter 7   Copyright (C) 2009 Jari Häkkinen
4207 26 Aug 22 peter 8   Copyright (C) 2012, 2022 Peter Johansson
193 27 Oct 04 peter 9
1437 25 Aug 08 peter 10   This file is part of the yat library, http://dev.thep.lu.se/yat
193 27 Oct 04 peter 11
675 10 Oct 06 jari 12   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 13   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 14   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 15   License, or (at your option) any later version.
675 10 Oct 06 jari 16
675 10 Oct 06 jari 17   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 18   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 19   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 20   General Public License for more details.
675 10 Oct 06 jari 21
675 10 Oct 06 jari 22   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 23   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 24 */
675 10 Oct 06 jari 25
2881 18 Nov 12 peter 26 #include <config.h>
2881 18 Nov 12 peter 27
1643 13 Dec 08 jari 28 #include "GSLInterpolation.h"
675 10 Oct 06 jari 29
1724 15 Jan 09 jari 30 #include "yat/utility/Exception.h"
1644 13 Dec 08 jari 31 #include "yat/utility/VectorBase.h"
1643 13 Dec 08 jari 32
1643 13 Dec 08 jari 33 #include <cassert>
1724 15 Jan 09 jari 34 #include <sstream>
1643 13 Dec 08 jari 35
193 27 Oct 04 peter 36 namespace theplu {
680 11 Oct 06 jari 37 namespace yat {
383 12 Aug 05 jari 38 namespace regression {
193 27 Oct 04 peter 39
703 18 Dec 06 jari 40
1643 13 Dec 08 jari 41   GSLInterpolation::GSLInterpolation(const gsl_interp_type* interpolator_type,
1644 13 Dec 08 jari 42                                      const utility::VectorBase& x,
1644 13 Dec 08 jari 43                                      const utility::VectorBase& y)
1724 15 Jan 09 jari 44     : evaluation_(0.0)
703 18 Dec 06 jari 45   {
1643 13 Dec 08 jari 46     assert(x.size()==y.size());
1643 13 Dec 08 jari 47     size_t size=x.size();
4187 18 Jul 22 peter 48     x_.reset(new double[size]);
4187 18 Jul 22 peter 49     std::copy(x.begin(), x.end(), x_.get());
4187 18 Jul 22 peter 50     y_.reset(new double[size]);
4187 18 Jul 22 peter 51     std::copy(y.begin(), y.end(), y_.get());
1644 13 Dec 08 jari 52     interpolator_ = gsl_interp_alloc(interpolator_type, size);
4187 18 Jul 22 peter 53     if (int status=gsl_interp_init(interpolator_, x_.get(), y_.get(), size)) {
1729 16 Jan 09 jari 54       std::stringstream ss;
1729 16 Jan 09 jari 55       ss << "GSLInterpolation::GSLInterpolation error code " << status;
1724 15 Jan 09 jari 56       // clean up
1724 15 Jan 09 jari 57       gsl_interp_free(interpolator_);
1724 15 Jan 09 jari 58       throw utility::GSL_error(ss.str());
1724 15 Jan 09 jari 59     }
1643 13 Dec 08 jari 60     accelerator_ = gsl_interp_accel_alloc();
703 18 Dec 06 jari 61   }
713 21 Dec 06 peter 62
718 26 Dec 06 jari 63
1643 13 Dec 08 jari 64   GSLInterpolation::~GSLInterpolation(void)
718 26 Dec 06 jari 65   {
1724 15 Jan 09 jari 66     gsl_interp_accel_free(accelerator_);
1643 13 Dec 08 jari 67     gsl_interp_free(interpolator_);
718 26 Dec 06 jari 68   }
718 26 Dec 06 jari 69
718 26 Dec 06 jari 70
1648 13 Dec 08 jari 71   double GSLInterpolation::evaluate(double x)
1650 14 Dec 08 jari 72   {
4187 18 Jul 22 peter 73     if (int status=gsl_interp_eval_e(interpolator_, x_.get(), y_.get(), x,
4187 18 Jul 22 peter 74                                      accelerator_, &evaluation_)) {
1729 16 Jan 09 jari 75       std::stringstream ss;
1730 16 Jan 09 jari 76       ss << "GSLInterpolation::evaluate for argument " << x << " fails.\n"
1730 16 Jan 09 jari 77          << "   Error code " << status;
1724 15 Jan 09 jari 78       throw utility::GSL_error(ss.str());
1724 15 Jan 09 jari 79     }
1650 14 Dec 08 jari 80     return evaluation_;
221 30 Dec 04 peter 81   }
221 30 Dec 04 peter 82
726 04 Jan 07 peter 83
1648 13 Dec 08 jari 84   double GSLInterpolation::evaluate_derivative(double x)
1648 13 Dec 08 jari 85   {
4187 18 Jul 22 peter 86     if (int status=gsl_interp_eval_deriv_e(interpolator_, x_.get(), y_.get(),
4187 18 Jul 22 peter 87                                            x, accelerator_, &evaluation_)) {
1729 16 Jan 09 jari 88       std::stringstream ss;
1730 16 Jan 09 jari 89       ss << "GSLInterpolation::evaluate_derivative for argument " << x
1730 16 Jan 09 jari 90          << " fails.\n   Error code " << status;
1724 15 Jan 09 jari 91       throw utility::GSL_error(ss.str());
1724 15 Jan 09 jari 92     }
1650 14 Dec 08 jari 93     return evaluation_;
1648 13 Dec 08 jari 94   }
1648 13 Dec 08 jari 95
1648 13 Dec 08 jari 96
1648 13 Dec 08 jari 97   double GSLInterpolation::evaluate_derivative2(double x)
1648 13 Dec 08 jari 98   {
4187 18 Jul 22 peter 99     if (int status=gsl_interp_eval_deriv2_e(interpolator_, x_.get(), y_.get(),
4187 18 Jul 22 peter 100                                             x, accelerator_, &evaluation_)) {
1729 16 Jan 09 jari 101       std::stringstream ss;
1730 16 Jan 09 jari 102       ss << "GSLInterpolation::evaluate_derivative2 for argument " << x
1730 16 Jan 09 jari 103          << " fails.\n   Error code " << status;
1724 15 Jan 09 jari 104       throw utility::GSL_error(ss.str());
1724 15 Jan 09 jari 105     }
1650 14 Dec 08 jari 106     return evaluation_;
1648 13 Dec 08 jari 107   }
1648 13 Dec 08 jari 108
1648 13 Dec 08 jari 109
1648 13 Dec 08 jari 110   double GSLInterpolation::evaluate_integral(double a, double b)
1648 13 Dec 08 jari 111   {
4187 18 Jul 22 peter 112     if (int status=gsl_interp_eval_integ_e(interpolator_, x_.get(), y_.get(),
4187 18 Jul 22 peter 113                                            a, b, accelerator_, &evaluation_)) {
1729 16 Jan 09 jari 114       std::stringstream ss;
1730 16 Jan 09 jari 115       ss << "GSLInterpolation::evaluate_integral for arguments a=" << a
1730 16 Jan 09 jari 116          << " and b=" << b << " fails.\n   Error code " << status;
1724 15 Jan 09 jari 117       throw utility::GSL_error(ss.str());
1724 15 Jan 09 jari 118     }
1650 14 Dec 08 jari 119     return evaluation_;
1648 13 Dec 08 jari 120   }
1648 13 Dec 08 jari 121
1648 13 Dec 08 jari 122
1650 14 Dec 08 jari 123   double GSLInterpolation::evaluation(void) const
1650 14 Dec 08 jari 124   {
1650 14 Dec 08 jari 125     return evaluation_;
1650 14 Dec 08 jari 126   }
1650 14 Dec 08 jari 127
1650 14 Dec 08 jari 128
1648 13 Dec 08 jari 129   unsigned int GSLInterpolation::min_size(void) const
1648 13 Dec 08 jari 130   {
1648 13 Dec 08 jari 131     return gsl_interp_min_size(interpolator_);
1648 13 Dec 08 jari 132   }
1648 13 Dec 08 jari 133
1648 13 Dec 08 jari 134
681 11 Oct 06 jari 135 }}} // of namespaces regression, yat, and theplu

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