xerus
a general purpose tensor library
simpleNumerics.cpp
Go to the documentation of this file.
1 // Xerus - A General Purpose Tensor Library
2 // Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf.
3 //
4 // Xerus is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU Affero General Public License as published
6 // by the Free Software Foundation, either version 3 of the License,
7 // or (at your option) any later version.
8 //
9 // Xerus is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU Affero General Public License for more details.
13 //
14 // You should have received a copy of the GNU Affero General Public License
15 // along with Xerus. If not, see <http://www.gnu.org/licenses/>.
16 //
17 // For further information on Xerus visit https://libXerus.org
18 // or contact us at contact@libXerus.org.
19 
26 
27 #include <xerus/misc/math.h>
29 #include <xerus/misc/exceptions.h>
30 #include <xerus/misc/namedLogger.h>
31 #include <xerus/misc/check.h>
32 #include <xerus/misc/internal.h>
33 
34 namespace xerus {
35  namespace misc {
36 
37  double integrate(const std::function<double(double)>& _f, const double _a, const double _b, double _eps, const uint _minIter, const uint _maxIter, const uint _branchFactor, const uint _maxRecursion, const bool _relativeError) {
38  REQUIRE(_minIter > 1, "");
39  REQUIRE(_branchFactor > 1, "");
40  double min = std::min(_a,_b);
41  double max = std::max(_a,_b);
42  if (_relativeError) { _eps = std::max(_eps, std::numeric_limits<double>::epsilon()); }
43  std::vector<double> iterants;
44  iterants.push_back((max - min)*(_f(min) + _f(max))/2);
45  double h = (max-min);
46  double sum;
47  double error = 1;
48  double maxVal = std::abs(iterants[0]);
49  for (uint iter=0; iter<_maxIter; ++iter) {
50  sum = 0;
51  for (double x = min+h/2; x<max; x += h) {
52  double f = _f(x);
53  sum += f;
54  maxVal = std::max(maxVal, std::abs(f));
55  }
56 
57  h /= 2;
58  sum *= h;
59  sum += iterants.back()/2;
60  iterants.push_back(sum);
61  double oldIt0 = iterants[0];
62  for (size_t k=0; k<iterants.size()-1; ++k) {
63  size_t i=iterants.size()-1-k;
64  iterants[i-1] = iterants[i] + (iterants[i]-iterants[i-1])/(misc::pow(2.0, 2*(k+1))-1);
65  }
66  if (_relativeError) {
67  error = std::abs((iterants[0]-oldIt0)/oldIt0);
68  if (std::isnan(error)) { error = std::abs(iterants[0]-oldIt0); }
69  } else {
70  error = std::abs(iterants[0]-oldIt0);
71  }
72  if (iter >= _minIter && error <= _eps) {
73  return (_a>_b?-1:1)*iterants[0];
74  }
75  }
76  if (_maxRecursion == 0) {
77  return (_a>_b?-1:1)*iterants[0];
78  }
79  // divide and conquer
80  // this dynamicaly divides those parts that did not converge easily
81  LOG(debug, (_relativeError?"relative":"absolute") << " error " << error << " is larger than eps = " << _eps);
82  sum = 0;
83  h = (max-min) / _branchFactor;
84  double newEps = (_relativeError
85  ? std::max(std::abs(iterants[0]), maxVal)*_eps
86  : std::max(
87  _eps, // /std::sqrt(_branchFactor),
88  std::sqrt(_branchFactor)*std::numeric_limits<double>::epsilon()*std::max(std::abs(iterants[0]), maxVal)
89  )
90  );
91  for (uint i=0; i<_branchFactor; ++i) {
92  sum += integrate(_f, min+i*h, min+(i+1)*h, newEps, _minIter, _maxIter, _branchFactor, _maxRecursion-1, false);
93  }
94  return (_a>_b?-1:1)*sum;
95  }
96 
97 
98  double integrate_segmented(const std::function<double(double)> &_f, const double _a, const double _b, const double _segmentation, const double _eps, const uint _minIter, const uint _maxIter, const uint _branchFactor, const uint _maxRecursion) {
99  const double min = std::min(_a, _b);
100  const double max = std::max(_a, _b);
101  double res = 0;
102  for (double x = min; x < max; x += _segmentation) {
103  res += integrate(_f, x, std::min(x+_segmentation, max), _eps, _minIter, _maxIter, _branchFactor, _maxRecursion);
104  }
105  return (_a>_b?-1:1)*res;
106  }
107 
108 
109  double find_root_bisection(const std::function<double(double)> &_f, double _min, double _max, const double _epsilon) {
110  const double nm = std::max(_min, _max);
111  _min = std::min(_min, _max);
112  _max = nm;
113 
114  double fmin = _f(_min);
115  double fmax = _f(_max);
116 
117  REQUIRE(fmin * fmax <= 0, "bisection requires inputs to both sides of the root (and no NaNs)");
118 
119  if (std::fpclassify(fmin) == FP_ZERO) {
120  return _min;
121  }
122 
123  if (std::fpclassify(fmax) == FP_ZERO) {
124  return _max;
125  }
126 
127  while (_max-_min > _epsilon) {
128  double mid = (_max+_min)/2.0;
129  double fmid = _f(mid);
130  if (std::fpclassify(fmid) == FP_ZERO) {
131  return mid;
132  }
133  REQUIRE(std::isfinite(fmid), "invalid function value f("<<mid<<") = " << fmid << " reached in bisection");
134  if (fmin * fmid < 0) {
135 // fmax = fmid; // never needed again
136  _max = mid;
137  } else {
138  fmin = fmid;
139  _min = mid;
140  }
141  }
142 
143  return (_max + _min) / 2.0;
144  }
145 
146  // - - - - - - - - - - - - - - - Polynomial - - - - - - - - - - - - - - -
147 
148  Polynomial::Polynomial() {}
149 
150  Polynomial::Polynomial(std::vector<double> _coeff) : coefficients(std::move(_coeff)) {}
151 
152  size_t Polynomial::terms() const {
153  return coefficients.size();
154  }
155 
156  Polynomial& Polynomial::operator+=(const Polynomial &_rhs) {
157  coefficients.resize(std::max(terms(), _rhs.terms()));
158 
159  for(size_t i = 0; i < _rhs.terms(); ++i) {
160  coefficients[i] += _rhs.coefficients[i];
161  }
162 
163  return *this;
164  }
165 
166  Polynomial& Polynomial::operator-=(const Polynomial &_rhs) {
167  coefficients.resize(std::max(terms(), _rhs.terms()));
168 
169  for(size_t i = 0; i < _rhs.terms(); ++i) {
170  coefficients[i] -= _rhs.coefficients[i];
171  }
172 
173  return *this;
174  }
175 
176 
177  Polynomial& Polynomial::operator/=(double _rhs) {
178  for (size_t i=0; i<terms(); ++i) {
179  coefficients[i] /= _rhs;
180  }
181  return *this;
182  }
183 
184  Polynomial& Polynomial::operator*=(double _rhs) {
185  for (size_t i=0; i<terms(); ++i) {
186  coefficients[i] *= _rhs;
187  }
188  return *this;
189  }
190 
191  Polynomial Polynomial::operator*(const Polynomial &_rhs) const {
192  Polynomial result;
193  result.coefficients.resize(_rhs.terms()+terms()-1);
194  for (size_t i=0; i<terms(); ++i) {
195  for (size_t j=0; j<_rhs.terms(); ++j) {
196  result.coefficients[i+j] += coefficients[i]*_rhs.coefficients[j];
197  }
198  }
199  return result;
200  }
201 
202  Polynomial Polynomial::operator*(double _rhs) const {
203  Polynomial result;
204  result.coefficients.resize(terms());
205  for (size_t i=0; i<terms(); ++i) {
206  result.coefficients[i] = coefficients[i]*_rhs;
207  }
208  return result;
209  }
210 
211  double Polynomial::operator()(double x) const {
212  double result = 0;
213  for (size_t i=terms(); i>0; --i) {
214  result *= x;
215  result += coefficients[i-1];
216  }
217  return result;
218  }
219 
220  double Polynomial::scalar_product(const Polynomial &_rhs, const std::function<double (double)> &_weight, double _minX, double _maxX) const {
221  return xerus::misc::integrate([&](double x){
222  return (*this)(x) * _rhs(x) * _weight(x);
223  }, _minX, _maxX, 1e-10);
224  }
225 
226  double Polynomial::norm(const std::function<double (double)> &_weight, double _minX, double _maxX) const {
227  return std::sqrt(scalar_product(*this, _weight, _minX, _maxX));
228  }
229 
231  Polynomial& Polynomial::orthogonolize(const std::vector<Polynomial> &_orthoBase, const std::function<double (double)> &_weight, double _minX, double _maxX) {
232  for (size_t i=0; i<_orthoBase.size(); ++i) {
233  (*this)-=_orthoBase[i]*scalar_product(_orthoBase[i], _weight, _minX, _maxX);
234  REQUIRE(std::abs(scalar_product(_orthoBase[i], _weight, _minX, _maxX)) < 1e-12, i << " " << std::abs(scalar_product(_orthoBase[i], _weight, _minX, _maxX)));
235  }
236  (*this)/=norm(_weight, _minX, _maxX);
237  return *this;
238  }
239 
241  std::vector<Polynomial> Polynomial::build_orthogonal_base(size_t _N, const std::function<double (double)> &_weight, double _minX, double _maxX) {
242  std::vector<Polynomial> base;
243  while (base.size() < _N) {
244  Polynomial next;
245  next.coefficients.resize(base.size()+1);
246  next.coefficients.back() = 1; // next = x^(base.size)
247  next.orthogonolize(base, _weight, _minX, _maxX);
248  base.push_back(next);
249  LOG(debug, "next basis polynomial " << next.coefficients);
250  }
251  return base;
252  }
253 
254 
255  // - - - - - - - - - - - - - - - ShanksTransformation - - - - - - - - - - - - - - -
256 
257  template<class ft_type>
258  ft_type ShanksTransformation<ft_type>::shanks(ft_type x1, ft_type x2, ft_type x3) {
259  ft_type denominator = x1 - 2*x2 + x3;
260  if (std::abs(denominator) < epsilon * std::max(x1, std::max(x2, x3))) {
261  return x2;
262  }
263  return (x1*x3 - x2*x2) / denominator;
264  }
265 
266  template<class ft_type>
268  values.push_back(_val);
269  for (size_t i=values.size()-1; i>=2; i-=2) {
270  values[i-2] = shanks(values[i-2], values[i-1], values[i]);
271  }
272  }
273 
274  template<class ft_type>
276  if (values.empty()) {
277  XERUS_THROW(generic_error() << "tried to extract limit of empty sequence");
278  } else {
279  return values[(values.size()-1) % 2];
280  }
281  }
282 
283  template<class ft_type>
285  size_t i = (values.size()-1) % 2;
286  if (i+1 >= values.size()) {
287  return 1;
288  }
289  return std::abs(values[i] - values[i+1]);
290  }
291 
292  template<class ft_type>
294  values.clear();
295  }
296 
297  template class ShanksTransformation<float>;
298  template class ShanksTransformation<double>;
299 
300 
301 
302 
303  // - - - - - - - - - - - - - - - ShanksTransformation - - - - - - - - - - - - - - -
304 
305  template<class ft_type>
306  ft_type RichardsonExtrapolation<ft_type>::richard(size_t n, ft_type x1, ft_type x2) {
307  return ft_type(n+1)*x2 - ft_type(n)*x1;
308  }
309 
310  template<class ft_type>
312  values.push_back(_val);
313  for (size_t i=values.size()-1; i>=1; i-=1) {
314  values[i-1] = richard(i-1, values[i-1], values[i]);
315  }
316  }
317 
318  template<class ft_type>
320  if (values.empty()) {
321  XERUS_THROW(generic_error() << "tried to extract limit of empty sequence");
322  } else {
323  return values.front();
324  }
325  }
326 
327  template<class ft_type>
329  if (values.size() < 2) {
330  return 1;
331  }
332  return std::abs(values[0] - values[1]);
333  }
334 
335  template<class ft_type>
337  values.clear();
338  }
339 
340 
341  template class RichardsonExtrapolation<float>;
342  template class RichardsonExtrapolation<double>;
343 
344  } // namespace misc
345 } // namespace xerus
Header file for some additional math functions.
Header file for CHECK and REQUIRE macros.
ft_type error_approximate() const override
ft_type error_approximate() const override
void push_back(ft_type _val) override
double integrate_segmented(const std::function< double(double)> &_f, double _a, double _b, double _segmentation, double _eps=1e-8, uint _minIter=4, uint _maxIter=6, uint _branchFactor=8, uint _maxRecursion=10)
Header file for the standard contaienr support functions.
#define REQUIRE
Definition: internal.h:33
#define LOG
Definition: internal.h:38
ft_type best_estimate() const override
The main namespace of xerus.
Definition: basic.h:37
static ft_type shanks(ft_type x1, ft_type x2, ft_type x3)
static ft_type richard(size_t n, ft_type x1, ft_type x2)
Header file for several elementary numerical methods: polynomials, romberg integration and limit extr...
void push_back(ft_type _val) override
The xerus exception class.
Definition: exceptions.h:37
Header file for xerus::misc::generic_error exception class.
Header file for all logging macros and log-buffer functionality.
unsigned int uint
Definition: standard.h:51
double find_root_bisection(const std::function< double(double)> &_f, double _min, double _max, double _epsilon=1e-14)
Header file for comfort functions and macros that should not be exported in the library.
#define XERUS_THROW(...)
Helper macro to throw a generic_error (or derived exception) with some additional information include...
Definition: exceptions.h:73
Limit extraction using the shanks transformation aka Aitken process.
Limit extraction using the richardson extrapolation.
ft_type best_estimate() const override
double integrate(const std::function< double(double)> &_f, double _a, double _b, double _eps=std::numeric_limits< double >::epsilon(), uint _minIter=4, uint _maxIter=6, uint _branchFactor=7, uint _maxRecursion=10, bool _relativeError=true)
Performs a Romberg Integration (richardson extrapolation of regular riemannian sum) + adaptive refine...
IndexedTensorMoveable< tensor_type > operator*(const value_t _factor, IndexedTensorMoveable< tensor_type > &&_tensor)
constexpr T pow(const T &_base, const uint32 _exp) noexcept
: Calculates _base^_exp by binary exponentiation
Definition: math.h:50