Medusa  1.1
Coordinate Free Mehless Method implementation
RBFBasis.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_APPROXIMATIONS_RBFBASIS_HPP_
2 #define MEDUSA_BITS_APPROXIMATIONS_RBFBASIS_HPP_
3 
9 #include <medusa/bits/utils/assert.hpp>
10 #include "RBFBasis_fwd.hpp"
11 
12 namespace mm {
13 
15 template <class vec_t, class RBFType>
16 typename RBFBasis<vec_t, RBFType>::scalar_t RBFBasis<vec_t, RBFType>::eval(
17  int index, const vector_t& point, const std::vector<vector_t>& support) const {
18  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
19  assert_msg(static_cast<int>(support.size()) >= size_,
20  "Not enough support points given for this basis size, got %d support points for "
21  "basis size %d.", support.size(), size_);
22  return rbf_((point-support[index]).squaredNorm());
23 }
24 
25 template <class RBFType, class vec_t>
26 template <typename operator_t>
27 typename vec_t::scalar_t RBFBasis<RBFType, vec_t>::evalOp(int index, const vector_t& point,
28  operator_t op, const std::vector<vector_t>& support, scalar_t scale) const {
29  return op.apply(*this, index, point, support, scale);
30 }
32 
33 template <class RBFType, class vec_t>
34 typename vec_t::scalar_t RBFBasis<RBFType, vec_t>::evalOp(
35  int index, const vector_t& point, Lap<dim> op, const std::vector<vector_t>& support,
36  scalar_t scale) const {
37  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
38  assert_msg(static_cast<int>(support.size()) >= size_,
39  "Not enough support points given for this basis size, got %d support points for "
40  "basis size %d.", support.size(), size_);
41  return rbf_((point - support[index]).squaredNorm(), op)/scale/scale;
42 }
43 
44 template <class RBFType, class vec_t>
45 typename vec_t::scalar_t RBFBasis<RBFType, vec_t>::evalOp(
46  int index, const vector_t& point, Der1<dim> op, const std::vector<vector_t>& support,
47  scalar_t scale) const {
48  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
49  assert_msg(static_cast<int>(support.size()) >= size_,
50  "Not enough support points given for this basis size, got %d support points for "
51  "basis size %d.", support.size(), size_);
52  vec_t p = point - support[index];
53  return 2*p[op.var]*rbf_(p.squaredNorm(), 1)/scale;
54 }
55 
56 template <class RBFType, class vec_t>
57 typename vec_t::scalar_t RBFBasis<RBFType, vec_t>::evalOp(
58  int index, const vector_t& point, Der2<dim> op, const std::vector<vector_t>& support,
59  scalar_t scale) const {
60  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
61  assert_msg(static_cast<int>(support.size()) >= size_,
62  "Not enough support points given for this basis size, got %d support points for "
63  "basis size %d.", support.size(), size_);
64  vec_t p = point - support[index];
65  scalar_t r2 = p.squaredNorm();
66  double res = 4*p[op.var1]*p[op.var2]*rbf_(r2, 2);
67  if (op.var1 == op.var2) {
68  res += 2*rbf_(r2, 1);
69  }
70  return res/scale/scale;
71 }
72 
73 template <class vec_t, class RBFType>
74 typename RBFBasis<vec_t, RBFType>::scalar_t RBFBasis<vec_t, RBFType>::evalAt0(
75  int index, const std::vector<vector_t>& support) const {
76  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
77  assert_msg(static_cast<int>(support.size()) >= size_,
78  "Not enough support points given for this basis size, got %d support points for "
79  "basis size %d.", support.size(), size_);
80  return rbf_(support[index].squaredNorm());
81 }
82 
83 template <class RBFType, class vec_t>
84 typename vec_t::scalar_t
85 RBFBasis<RBFType, vec_t>::evalOpAt0(int index, const Lap<dim>& lap,
86  const std::vector<vector_t>& support, scalar_t scale) const {
87  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
88  assert_msg(static_cast<int>(support.size()) >= size_,
89  "Not enough support points given for this basis size, got %d support points for "
90  "basis size %d.", support.size(), size_);
91  return rbf_(support[index].squaredNorm(), lap)/scale/scale;
92 }
93 
94 template <class RBFType, class vec_t>
95 typename vec_t::scalar_t RBFBasis<RBFType, vec_t>::evalOpAt0(int index, const Der1<dim>& der1,
96  const std::vector<vector_t>& support, scalar_t scale) const {
97  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
98  assert_msg(static_cast<int>(support.size()) >= size_,
99  "Not enough support points given for this basis size, got %d support points for "
100  "basis size %d.", support.size(), size_);
101  assert_msg(der1.var >= 0, "Index of derived variable %d should be higher or equal to 0.",
102  der1.var);
103  assert_msg(der1.var < dim, "Index of derived variable %d should be lower"
104  " than the number of dimensions %d", der1.var, dim);
105 
106  const auto& s = support[index];
107  scalar_t r2 = s.squaredNorm();
108  scalar_t v = rbf_(r2, 1);
109  v *= -2*s[der1.var];
110  return v/scale;
111 }
112 
113 template <class RBFType, class vec_t>
114 typename vec_t::scalar_t RBFBasis<RBFType, vec_t>::evalOpAt0(int index, const Der2<dim>& der2,
115  const std::vector<vector_t>& support, scalar_t scale) const {
116  assert_msg(0 <= index && index < size_, "Basis index %d out of range [0, %d)", index, size_);
117  assert_msg(static_cast<int>(support.size()) >= size_,
118  "Not enough support points given for this basis size, got %d support points for "
119  "basis size %d.", support.size(), size_);
120  assert_msg(der2.var1 >= 0 && der2.var2 >=0, "Indexes of derived variables %d and %d should"
121  " be both higher or equal to 0.",
122  der2.var1, der2.var2);
123  assert_msg(der2.var1 < dim && der2.var2 < dim , "Index of derived variable %d and %d should"
124  " both be lower than the number of dimensions %d", der2.var1, der2.var2, dim);
125 
126  auto& s = support[index];
127  scalar_t r2 = s.squaredNorm();
128  int d1 = der2.var1;
129  int d2 = der2.var2;
130  scalar_t v = 4*s[d1]*s[d2]*rbf_(r2, 2);
131  if (d1 == d2) {
132  v += 2*rbf_(r2, 1);
133  }
134  return v/scale/scale;
135 }
136 
138 template <typename V, typename R>
139 std::ostream& operator<<(std::ostream& os, const RBFBasis<V, R>& m) {
140  return os << "RBF basis " << m.dim << "D, " << m.size() << "x " << m.rbf();
141 }
142 
143 } // namespace mm
144 
145 #endif // MEDUSA_BITS_APPROXIMATIONS_RBFBASIS_HPP_
mm::sh::d1
static const shape_flags d1
Indicates to calculate d1 shapes.
Definition: shape_flags.hpp:23
RBFBasis_fwd.hpp
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
scalar_t
Scalar scalar_t
Type of the elements, alias of Scalar.
Definition: MatrixBaseAddons.hpp:16
mm::RBFBasis::vector_t
vec_t vector_t
Vector type.
Definition: RBFBasis_fwd.hpp:35
mm::sh::lap
static const shape_flags lap
Indicates to calculate laplace shapes.
Definition: shape_flags.hpp:24
dim
@ dim
Number of elements of this matrix.
Definition: MatrixBaseAddons.hpp:14
mm::RBFBasis::dim
@ dim
Dimensionality of the function domain.
Definition: RBFBasis_fwd.hpp:38
mm::operator<<
std::ostream & operator<<(std::ostream &os, const Gaussian< S > &b)
Output basic information about given Gaussian RBF.
Definition: Gaussian.hpp:37
mm::RBFBasis::eval
scalar_t eval(int index, const vector_t &point, const std::vector< vector_t > &support) const
Evaluates index-th RBF's at point.
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
mm::Der1::var
int var
Index representing derived variable (x = 0, y = 1, z = 2)
Definition: Operators_fwd.hpp:144
mm::RBFBasis::evalAt0
scalar_t evalAt0(int index, const std::vector< vector_t > &support) const
Evaluate index-th function at zero.
Definition: RBFBasis.hpp:74
mm::sh::d2
static const shape_flags d2
Indicates to calculate d2 shapes.
Definition: shape_flags.hpp:25
mm::RBFBasis::scalar_t
vector_t::scalar_t scalar_t
Scalar type.
Definition: RBFBasis_fwd.hpp:36
mm::RBFBasis::evalOp
scalar_t evalOp(int index, const vector_t &point, operator_t op, const std::vector< vector_t > &support, scalar_t scale=1) const
Apply an operator at a given point.
mm::Der2::var2
int var2
Higher index representing derived variable (x = 0, y = 1, z = 2, ...)
Definition: Operators_fwd.hpp:168
mm::Der2
Represents a second derivative wrt. var1 and var2.
Definition: Monomials_fwd.hpp:24
assert.hpp
mm::RBFBasis
Represents a basis of Radial Basis Functions over a local neighbourhood.
Definition: RBFBasis_fwd.hpp:30
mm::Lap
Represents the Laplacian operator.
Definition: Monomials_fwd.hpp:20
mm::Der1
Represents a first derivative wrt. var.
Definition: Monomials_fwd.hpp:22
mm::RBFBasis::size
int size() const
Returns basis size.
Definition: RBFBasis_fwd.hpp:51
mm::RBFBasis::rbf
const rbf_t & rbf() const
Returns underlying RBF object.
Definition: RBFBasis_fwd.hpp:54
mm::RBFBasis::evalOpAt0
scalar_t evalOpAt0(int index, const operator_t&op, const std::vector< vector_t > &support, scalar_t scale=1) const
Evaluate general operator op at zero.
Definition: RBFBasis_fwd.hpp:115
mm::Der2::var1
int var1
Lower index representing derived variable (x = 0, y = 1, z = 2, ...)
Definition: Operators_fwd.hpp:167