1 #pragma once
2
3 /// \file Functional.h
4 /// \author Basic routines for integrating arbitrary functions in N dimensions, finding roots, etc.
5 /// \author Pavel Sevecek (sevecek at sirrah.troja.mff.cuni.cz)
6 /// \date 2016-2021
7
8 #include "math/rng/Rng.h"
9 #include "objects/geometry/Domain.h"
10 #include "objects/wrappers/AutoPtr.h"
11
12 NAMESPACE_SPH_BEGIN
13
14 /// \brief Integrates a one-dimensional function using Simpson's rule.
15 ///
16 /// The precision of the integral is currently fixed.
17 template <typename TFunctor>
integrate(const Interval range,const TFunctor & functor)18 INLINE Float integrate(const Interval range, const TFunctor& functor) {
19 const Size N = 1000;
20 const Float h = range.size() / N;
21 double result = functor((Float)range.lower()) + functor((Float)range.upper());
22 for (Size j = 1; j < N; ++j) {
23 const Float x = (Float)range.lower() + j * h;
24 if (j % 2 == 0) {
25 result += 2._f * functor(x);
26 } else {
27 result += 4._f * functor(x);
28 }
29 }
30 return Float(result * h / 3._f);
31 }
32
33
34 /// \brief Object for integrating a generic three-dimensional scalar function.
35 template <typename TRng = UniformRng, typename TInternal = double>
36 class Integrator : public Noncopyable {
37 private:
38 const IDomain& domain;
39 TRng rng;
40
41 static constexpr Size chunk = 100;
42
43 public:
44 /// \brief Constructs an integrator given the domain of integration.
Integrator(const IDomain & domain)45 Integrator(const IDomain& domain)
46 : domain(domain) {}
47
48 Integrator(IDomain&& domain) = delete;
49
50 /// \brief Integrates a function.
51 ///
52 /// \param f Functor with a Vector parameter, returning the integral as a scalar value
53 template <typename TFunctor>
54 Float integrate(TFunctor&& f, const Float targetError = 0.001_f) {
55 double sum = 0.;
56 double sumSqr = 0.;
57 double errorSqr = sqr(targetError);
58 Size n = 0;
59 StaticArray<Vector, chunk> buffer;
60 Array<Size> inside;
61 Box box = domain.getBoundingBox();
62 while (true) {
63 for (Size i = 0; i < chunk; ++i) {
64 // dimensionless vector
65 const Vector q(rng(0), rng(1), rng(2));
66 // vector properly scaled and moved
67 const Vector v = box.lower() + q * box.size();
68 buffer[i] = v;
69 }
70 inside.clear();
71 domain.getSubset(buffer, inside, SubsetType::INSIDE);
72 for (Size i : inside) {
73 const double x = (double)f(buffer[i]);
74 sum += x;
75 sumSqr += x * x;
76 n++;
77 }
78 const double m = double(n);
79 if (m * sumSqr - sum * sum < m * m * errorSqr * sumSqr) {
80 return Float(sum / m) * domain.getVolume();
81 }
82 }
83 }
84 };
85
86 /// \brief Returns a root of a R->R function on given range.
87 ///
88 /// If there is no root or the root cannot be found, returns NOTHING. For functions with multiple roots,
89 /// returns one of them; the selection of such root is not specified.
90 template <typename TFunctor>
getRoot(const Interval & range,const Float eps,const TFunctor & functor)91 INLINE Optional<Float> getRoot(const Interval& range, const Float eps, const TFunctor& functor) {
92 Interval r = range;
93 if (functor(r.lower()) * functor(r.upper()) > 0._f) { // same sign
94 return NOTHING;
95 }
96 while (r.size() > eps * range.size()) {
97 Float x = r.center();
98 if (functor(x) * functor(r.upper()) > 0._f) {
99 r = Interval(r.lower(), x);
100 } else {
101 r = Interval(x, r.upper());
102 }
103 }
104 return r.center();
105 }
106
107 /// \brief Checks if the given R->R function is continuous in given interval.
108 ///
109 /// The eps and delta values are the same for all values and must be specified by the caller.
110 template <typename TFunctor>
isContinuous(const Interval & range,const Float delta,const Float eps,const TFunctor & functor)111 INLINE bool isContinuous(const Interval& range, const Float delta, const Float eps, const TFunctor& functor) {
112 Float y1 = functor(range.lower());
113 for (Float x = range.lower(); x <= range.upper(); x += delta) {
114 const Float y2 = functor(x);
115 if (abs(y1 - y2) > eps) {
116 return false;
117 }
118 y1 = y2;
119 }
120 return true;
121 }
122
123 NAMESPACE_SPH_END
124