Flow123d  master-f44eb46
interpolant_impl.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file interpolant_impl.hh
15  * @brief
16  */
17 
18 #ifndef INTERPOLANT_IMPL_H
19 #define INTERPOLANT_IMPL_H
20 
21 #include <cmath>
22 
23 #include "interpolant.hh"
24 
25 /********************************** InterpolantBase ********************************/
26 
27 inline double InterpolantBase::error()
28 {
29  return error_;
30 }
31 
33 {
34  return stats;
35 }
36 
37 inline double InterpolantBase::bound_a() const
38 {
39  return bound_a_;
40 }
41 
42 inline double InterpolantBase::bound_b() const
43 {
44  return bound_b_;
45 }
46 
47 inline unsigned int InterpolantBase::size() const
48 {
49  return size_;
50 }
51 
52 
53 /********************************** Interpolant ********************************/
54 
55 template<template<class> class Func, class Type >
56 Interpolant::Interpolant(Func<Type>* func, bool interpolate_derivative)
57  : func(func),
58  interpolate_derivative(interpolate_derivative)
59 {
60  func_diff = new Func<B<Type> >();
61  func_diffn = new Func<T<Type> >();
62 
65 
66  checks[Check::functor] = true;
67 }
68 
69 
70 template<template<class> class Func, class Type >
71 void Interpolant::set_functor(Func<Type>* func, bool interpolate_derivative)
72 {
73  this->interpolate_derivative = interpolate_derivative;
74  this->func = func;
75  func_diff = new Func<B<Type> >();
76  func_diffn = new Func<T<Type> >();
79 
80  checks[Check::functor] = true;
81 }
82 
83 
84 inline double Interpolant::val(double x)
85 {
86  //increase calls
88 
89  //uncomment if we want to do the check automatically
90  //every STATISTIC_CHECK calls check the statistics
91  //if((stats.total_calls % STATISTIC_CHECK == 0) && use_statistics)
92  // this->check_and_reinterpolate();
93 
94  //return value
95  if(x < bound_a_) //left miss
96  {
97  //DebugOut() << "test\n";
98  //std::cout << "test a" << std::endl;
100  stats.min = std::min(stats.min, x);
101 
102  switch(extrapolation)
103  {
105  return f_vec[0];
107  return (f_vec[1] - f_vec[0])*(x/step-a_div_step) + f_vec[0]; //f_vec[0] + p1_vec[0]*(x-x_vec[0]);
109  default:
110  return f_val(x); //otherwise compute original functor
111  }
112  }
113  else if(x > bound_b_) //right miss
114  {
115  //DebugOut() << "test\n";
116  //std::cout << "test b" << std::endl;
118  stats.max = std::max(stats.max, x);
119 
120  switch(extrapolation)
121  {
123  return f_vec[size_];
125  return (f_vec[size_] - f_vec[size_-1])*(x-bound_b_+step)/step + f_vec[size_-1];//f_vec[size_-1] + p1_vec[size_-1]*(x-x_vec[size_-1]);
127  default:
128  return f_val(x); //otherwise compute original functor
129  }
130  }
131  else //hit the interval
132  {
133  return val_p1(x);
134  }
135 }
136 
137 inline double Interpolant::val_test(double x)
138 {
139  /*
140  //increase calls
141  stats.total_calls++;
142 
143  if(x < bound_a_)
144  {
145  //std::cout << "test a" << std::endl;
146  return 0;
147  }
148  else if(x > bound_b_) //right miss
149  {
150  //std::cout << "test b" << std::endl;
151  return 0;
152  }
153  else //*/
154  return val_p1(x);
155 }
156 
157 inline DiffValue Interpolant::diff(double x)
158 {
159  ASSERT_PERMANENT(interpolate_derivative).error("Derivative is not interpolated. Flag must be switched true in constructor (or set_functor).");
160  //increase calls
161  stats.total_calls++;
162 
163  //uncomment if we want to do the check automatically
164  //every STATISTIC_CHECK calls check the statistics
165  //if((stats.total_calls % STATISTIC_CHECK == 0) && use_statistics)
166  // this->check_and_reinterpolate();
167 
168  if(x < bound_a_) //left miss
169  {
171  stats.min = std::min(stats.min, x);
172  switch(extrapolation)
173  {
175  return DiffValue(f_vec[0],
176  df_vec[0]);
178  return DiffValue((f_vec[1] - f_vec[0])*(x/step-a_div_step) + f_vec[0], //f_vec[0] + p1_vec[0]*(x-x_vec[0]),
179  (df_vec[1] - df_vec[0])*(x/step-a_div_step) + df_vec[0]);//df_vec[0] + p1d_vec[0]*(x-x_vec[0]));
181  default:
182  return f_diff(x); //otherwise compute original functor
183  }
184  }
185  else if(x > bound_b_) //right miss
186  {
188  stats.max = std::max(stats.max, x);
189  switch(extrapolation)
190  {
192  return DiffValue(f_vec[size_],
193  df_vec[size_]);
195  return DiffValue((f_vec[size_] - f_vec[size_-1])*(x-bound_b_+step)/step + f_vec[size_-1], //f_vec[size_-1] + p1_vec[size_-1]*(x-x_vec[size_-1]),
196  (df_vec[size_] - df_vec[size_-1])*(x-bound_b_+step)/step + df_vec[size_-1]);//df_vec[size_-1] + p1d_vec[size_-1]*(x-x_vec[size_-1]) );
198  default:
199  return f_diff(x); //otherwise compute original functor
200  }
201  }
202  else //hit the interval
203  {
204  return diff_p1(x);
205  }
206 }
207 
208 inline double Interpolant::f_val(double x)
209 {
210  return func->operator()(x);
211 }
212 
214 {
215  B<double> xx(x); // Initialize arguments
216  //Func func; // Instantiate functor
217  B<double> f(func_diff->operator()(xx)); // Evaluate function and record DAG
218  f.diff(0,1); // Differentiate
219 
220  DiffValue d;
221  d.first = f.x(); // Value of function
222  d.second = xx.d(0); // Value of df/dx
223  return d; // Return function value
224 }
225 
226 /*
227 inline unsigned int Interpolant::find_interval(double x)
228 {
229  //counts in which interval x is (the last node before x)
230  return floor((x - bound_a_) / step);
231 }
232 */
233 
234 /* //CONSTANT INTERPOLATION
235 inline double Interpolant::val_p0(double x)
236 {
237  return f_vec[find_interval(x)];
238 }
239 
240 inline DiffValue Interpolant::diff_p0(double x)
241 {
242  DiffValue result;
243  unsigned int i = find_interval(x);
244  result.f = f_vec[i];
245  result.dfdx = df_vec[i];
246  return result;
247 }
248 */
249 
250 inline double Interpolant::val_p1(double x)
251 {
252  /*
253 [ OK ] InterpolantTest.InterpolantEval_val (184 ms)
254 [ OK ] InterpolantTest.InterpolantEval_val_p1 (192 ms)
255 [ OK ] InterpolantTest.InterpolantEval_val_test (190 ms)
256 
257  When we used precomputed coeficients of P1 interpolation.
258  */
259  //unsigned int i = floor(x / step - a_div_step);
260  //return p1_vec[i]*(x-x_vec[i]) + f_vec[i];
261 
262  /*
263 [ OK ] InterpolantTest.InterpolantEval_val (162 ms)
264 [ OK ] InterpolantTest.InterpolantEval_val_p1 (163 ms)
265 [ OK ] InterpolantTest.InterpolantEval_val_test (162 ms)
266 
267 These results are satisfying:
268  + the function val_p1 itself is by 14% faster than before
269  + we use only one array to keep the function values
270  + all three functions take the same time which means
271  that the collecting statistics and going through several IF commands
272  does not affect the evaluation time
273  */
274  double x_step = x / step;
275  unsigned int i = floor(x_step - a_div_step);
276 
277  return (f_vec[i+1] - f_vec[i])*(x_step-i-a_div_step) + f_vec[i];
278 }
279 
280 
282 {
283  DiffValue result;
284  //unsigned int i = find_interval(x);
285  //result.first = p1_vec[i]*(x-x_vec[i]) + f_vec[i];
286  //result.second = p1d_vec[i]*(x-x_vec[i]) + df_vec[i];
287 
288  double x_step = x / step;
289  unsigned int i = floor(x_step - a_div_step);
290 
291  result.first = (f_vec[i+1] - f_vec[i])*(x_step-i-a_div_step) + f_vec[i];
292  result.second = (df_vec[i+1] - df_vec[i])*(x_step-i-a_div_step) + df_vec[i];
293 
294  return result;
295 }
296 
297 
298 /** Functor class that computes the argument of the integral \f$ (f(x)-i(x))^p + (f'(x)-i'(x))^p \f$ in the norm \f$ \|f-i\|_{W^1_p} \f$.
299  * It is used as input functor to integration.
300  */
302 {
303 public:
305  : interpolant(interpolant), p_(p), tol_(tol) {}
306 
307  inline double p() {return p_;}
308  inline double tol() {return tol_;}
309 
310  virtual double operator()(double x)
311  {
312  double f_val = interpolant->f_val(x);
313  double a = std::abs(f_val - interpolant->val(x)) / (std::abs(f_val) + tol_);
314 
315  return std::pow(a, p_);
316  }
317 
318 private:
320  double p_;
321  double tol_;
322 };
323 
324  /** Functor class that computes the argument of the integral \f$ (f(x)-i(x))^p + (f'(x)-i'(x))^p \f$ in the norm \f$ \|f-i\|_{W^1_p} \f$.
325  * It is used as input functor to integration.
326  */
328 {
329 public:
331  : interpolant(interpolant), p_(p), tol_(tol) {}
332 
333  inline double p() {return p_;}
334  inline double tol() {return tol_;}
335 
336  virtual double operator()(double x)
337  {
338  DiffValue f = interpolant->f_diff(x);
339  DiffValue g = interpolant->diff(x);
340  double a = std::abs(f.first - g.first) / (std::abs(f.first) + tol_)
341  + std::abs(f.second - g.second) / (std::abs(f.second) + tol_);
342 
343  return std::pow(a, p_);
344  }
345 
346 private:
348  double p_;
349  double tol_;
350 };
351 
352 
353 
354 /********************************** InterpolantImplicit ********************************/
355 
356 template<template<class> class Func, class Type >
357 void InterpolantImplicit::set_functor(Func<Type>* func, bool interpolate_derivative)
358 {
359  this->func = func;
360  func_diff = new Func<B<Type> >();
361  func_diffn = new Func<T<Type> >();
362 
365 
366  this->interpolate_derivative = interpolate_derivative;
367 
368  checks[Check::functor] = true;
369 }
370 
371 template<template<class> class Func, class Type >
372 InterpolantImplicit::InterpolantImplicit(Func<Type>* func, bool interpolate_derivative)
373  : func(func),
374  interpolate_derivative(interpolate_derivative),
375  explicit_interpolant(NULL),
376  fix_(IFixVariable::no_fix)
377 {
378  func_diff = new Func<B<Type> >();
379  func_diffn = new Func<T<Type> >();
380 
383 
384  checks[Check::functor] = true;
385 }
386 
387  ///class FuncExplicit.
388  /** This functor transforms implicit functor with two variables into
389  * an explicit functor with only one variable and the other one fixed.
390  */
391  template<class Type>
393  {
394  public:
395  ///Constructor.
397 
398  //constructor from templated implicit functor
399  template<class TType>
401  : func_impl(&func_impl), fix_(fix), fix_val(fix_val){}
402 
403  Type operator()(Type u)
404  {
406  Type ret = 0;
408  ret = func_impl->operator()(fix_val,u);
410  ret = func_impl->operator()(u,fix_val);
411 
412  return ret;
413  }
414 
415  private:
418  double fix_val;
419  };
420 
421 
422 inline double InterpolantImplicit::val(double u)
423 {
424  return explicit_interpolant->val(u);
425 }
426 
428 {
429  return explicit_interpolant->diff(u);
430 }
431 
432 #endif //INTERPOLATION_IMPL_H
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
Abstract templated explicit functor class.
Definition: functors.hh:108
void set_param_from_func(FunctorCommon< TType > *func)
Sets a functor's parameters from another functor.
Abstract templated implicit functor class.
Definition: functors.hh:129
unsigned int size() const
Returns the size of the interpolation table.
double bound_a() const
Returns left boundary of the interval.
double bound_b_
Right interval boundary.
Definition: interpolant.hh:168
double bound_b() const
Returns right boundary of the interval.
double step
Chosen interpolation step.
Definition: interpolant.hh:169
std::vector< bool > checks
Vector of boolean values telling us which parameters are set or not.
Definition: interpolant.hh:206
EvalStatistics statistics() const
Returns structure with evaluation statistics.
double error_
Error of the interpolation.
Definition: interpolant.hh:181
double error()
Returns error of the interpolation.
Extrapolation::Type extrapolation
Definition: interpolant.hh:186
double bound_a_
Left interval boundary.
Definition: interpolant.hh:167
unsigned int size_
Number of dividing intervals.
Definition: interpolant.hh:172
double a_div_step
bound_ divided by step - precomputed value for evaluation
Definition: interpolant.hh:170
EvalStatistics stats
Structure which keeps evaluation statistics. See InterpolantBase::eval_statistics.
Definition: interpolant.hh:214
FuncExplicit(IFunctorBase< TType > &func_impl, IFixVariable::Type fix, double fix_val)
Type operator()(Type u)
Virtual operator () with type Type.
DiffValue diff(double u)
Returns interpolated value of the derivation.
IFunctorBase< B< double > > * func_diff
Definition: interpolant.hh:526
InterpolantImplicit()
constructor
Definition: interpolant.cc:537
Interpolant * explicit_interpolant
Definition: interpolant.hh:531
IFunctorBase< double > * func
Definition: interpolant.hh:525
IFunctorBase< T< double > > * func_diffn
Definition: interpolant.hh:527
bool interpolate_derivative
Is true if we want to interpolate the derivative too.
Definition: interpolant.hh:529
double val(double u)
Returns interpolated value.
void set_functor(Func< Type > *func, bool interpolate_derivative=false)
Sets the implicit functor.
virtual double operator()(double x)
Virtual operator () with type Type.
FuncError_lp(Interpolant *interpolant, double p, double tol)
FuncError_wp1(Interpolant *interpolant, double p, double tol)
virtual double operator()(double x)
Virtual operator () with type Type.
The main class for interpolation of functors.
Definition: interpolant.hh:292
FunctorBase< T< double > > * func_diffn
Pointer to original functor with FADBAD T type.
Definition: interpolant.hh:384
std::vector< double > f_vec
Vector of function values at nodes.
Definition: interpolant.hh:389
FunctorBase< B< double > > * func_diff
Pointer to original functor with FADBAD B type.
Definition: interpolant.hh:383
DiffValue diff_p1(double x)
Finds interval on which x lies.
double f_val(double x)
Returns value of the original functor.
double val_test(double x)
Do NOT use, only for testing purpose.
std::vector< double > df_vec
Vector of function derivatives values at nodes.
Definition: interpolant.hh:390
bool interpolate_derivative
Is true if we want to interpolate the derivative too.
Definition: interpolant.hh:386
double val(double x)
Returns interpolated value.
double val_p1(double x)
Do NOT use, unless you are 100% sure.
Interpolant()
Default constructor.
Definition: interpolant.cc:155
DiffValue f_diff(double x)
Returns 1st derivative of original functor using FADBAD.
DiffValue diff(double x)
Returns interpolated value of the derivation.
FunctorBase< double > * func
Pointer to original functor with double type.
Definition: interpolant.hh:380
void set_functor(Func< Type > *func, bool interpolate_derivative=false)
Sets the functor.
std::pair< double, double > DiffValue
Definition: interpolant.hh:26
@ linear
means that the linear approximation on the first and last subintervals will be used outside the inter...
Definition: interpolant.hh:43
@ functor
means that the original functor will be called outside interpolation interval
Definition: interpolant.hh:44
@ constant
means that the values at the boundary points will be used outside the interval_miss_a
Definition: interpolant.hh:42
@ fix_x
x variable will be fixed with given value
Definition: interpolant.hh:64
@ fix_y
y variable will be fixed with given value
Definition: interpolant.hh:65
@ no_fix
no variable is fixed (used when InterpolantImplicit is created)
Definition: interpolant.hh:66
Structure that keeps statistics of evaluation.
Definition: interpolant.hh:144
double max
maximal x for which the evaluation was called outside the interval ((initially is equal the right bou...
Definition: interpolant.hh:150
unsigned int interval_miss_a
counts left misses of the interval
Definition: interpolant.hh:145
unsigned int interval_miss_b
counts right misses of the interval
Definition: interpolant.hh:146
unsigned int total_calls
counts total calls of evaluation
Definition: interpolant.hh:147
double min
minimal x for which the evaluation was called outside the interval (initially is equal the left bound...
Definition: interpolant.hh:149