1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA  02110-1301  USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 #ifndef UQ_1D_1D_FUNCTION_H
26 #define UQ_1D_1D_FUNCTION_H
27 
28 #include <queso/ScalarSequence.h>
29 #include <queso/1DQuadrature.h>
30 #include <queso/Environment.h>
31 #include <queso/Defines.h>
32 #include <vector>
33 #include <math.h>
34 #include <fstream>
35 
36 namespace QUESO {
37 
38 /*! \file 1D1DFunction.h
39     \brief One-dimension function class.
40 */
41 
42 
43 
44 //*****************************************************
45 // Base 1D->1D class
46 //*****************************************************
47 
48 /*! \class Base1D1DFunction
49     \brief Class for one-dimensional functions.
50 
51     Base class for one-dimensional functions.
52 */
53 class
54 Base1D1DFunction {
55 public:
56   //! @name Constructor/Destructor methods
57   //@{
58   //! Default constructor.
59   Base1D1DFunction(double minDomainValue,
60 			  double maxDomainValue);
61 
62   //! Destructor.
63   virtual ~Base1D1DFunction();
64   //@}
65 
66   //! @name Mathematical  methods
67   //@{
68   //! Returns the minimum value of the domain of the (one-dimensional) function.
69   double minDomainValue() const;
70 
71   //! Returns the maximum value of the domain of the (one-dimensional) function.
72   double maxDomainValue() const;
73 
74   //! Returns the value of the (one-dimensional) function. See template specialization.
75   virtual  double value         (double domainValue) const = 0;
76 
77   //! Returns the value of the derivative of the function. See template specialization.
78   virtual  double deriv         (double domainValue) const = 0;
79 
80   //! TODO: Multiplies \c this function with \c function, and integrates it numerically.  See template specialization.
81   /*! \todo: Please, implement me!*/
82   virtual  double multiplyAndIntegrate(const Base1D1DFunction& func, unsigned int quadratureOrder, double* resultWithMultiplicationByTAsWell) const;
83   //@}
84 protected:
85   double m_minDomainValue;
86   double m_maxDomainValue;
87 };
88 
89 //*****************************************************
90 // Generic 1D->1D class
91 //*****************************************************
92 /*! \class Generic1D1DFunction
93     \brief Class for generic one-dimensional functions.
94 
95     This class implements generic one-dimensional functions.
96 */
97 class Generic1D1DFunction : public Base1D1DFunction {
98 public:
99   //! @name Constructor/Destructor methods
100   //@{
101   //! Default constructor.
102   Generic1D1DFunction(double minDomainValue,
103                              double maxDomainValue,
104                              double (*valueRoutinePtr)(double domainValue, const void* routinesDataPtr),
105                              double (*derivRoutinePtr)(double domainValue, const void* routinesDataPtr),
106                              const void* routinesDataPtr);
107   //! Destructor.
108   ~Generic1D1DFunction();
109   //@}
110 
111   //! @name Mathematical  methods
112   //@{
113   //! Returns the value of the (one-dimensional) function at point \c domainValue.
114   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
115    * and in affirmative case, it evaluates the function at such point. */
116   double value(double domainValue) const;
117 
118   //! Returns the value of the derivative of the function at point \c domainValue.
119   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
120    * and in affirmative case, it returns the value of the derivative at such point. */
121   double deriv(double domainValue) const;
122   //@}
123 protected:
124   using Base1D1DFunction::m_minDomainValue;
125   using Base1D1DFunction::m_maxDomainValue;
126 
127   double (*m_valueRoutinePtr)(double domainValue, const void* routinesDataPtr);
128   double (*m_derivRoutinePtr)(double domainValue, const void* routinesDataPtr);
129   const void* m_routinesDataPtr;
130 };
131 
132 //*****************************************************
133 // Constant 1D->1D class
134 //*****************************************************
135 /*! \class Constant1D1DFunction
136     \brief Class for constant one-dimensional functions.
137 
138     This class implements constant one-dimensional functions.
139 */
140 class Constant1D1DFunction : public Base1D1DFunction {
141 public:
142   //! @name Constructor/Destructor methods
143   //@{
144   //! Default constructor.
145   Constant1D1DFunction(double minDomainValue,
146                               double maxDomainValue,
147                               double constantValue);
148   //! Destructor.
149   ~Constant1D1DFunction();
150   //@}
151 
152   //! @name Mathematical  methods
153   //@{
154   //! Returns the value of the constant function at point \c domainValue.
155   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
156    * and in affirmative case, it evaluates the function at such point, which is the constant
157    * value \c constantValue passed to the constructor of this class. */
158   double value(double domainValue) const;
159 
160   //! Returns the value of the derivative of the constant function at point \c domainValue.
161   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
162    * and in affirmative case, it returns 0 (derivative of a constant function. */
163   double deriv(double domainValue) const;
164 //@}
165 protected:
166   using Base1D1DFunction::m_minDomainValue;
167   using Base1D1DFunction::m_maxDomainValue;
168 
169   double m_constantValue;
170 };
171 
172 //*****************************************************
173 // Linear 1D->1D class
174 //*****************************************************
175 /*! \class Linear1D1DFunction
176     \brief Class for linear one-dimensional functions.
177 
178     This class implements linear one-dimensional functions.
179     A common linear function is the the equation of a line: \f[f(x) = y_1 + m (x - x_1), \f]
180     which is a linear function with slope (or rate) \f$ m \f$ and passing through the point
181     \f$(x_1,y_1)\f$.
182 
183     This class implements a linear function such as the one describe above where the rate is
184     given by \c rateValue, and the point is <c>(referenceDomainValue,referenceImageValue)</c>.
185 */
186 class Linear1D1DFunction : public Base1D1DFunction {
187 public:
188   //! @name Constructor/Destructor methods
189   //@{
190   //! Default constructor.
191   /*! Constructs a linear function of constant rate of change \c rateValue which passes through
192    *  the point (x,y)=<c>(referenceDomainValue,referenceImageValue)</c>.*/
193   Linear1D1DFunction(double minDomainValue,
194                             double maxDomainValue,
195                             double referenceDomainValue,
196                             double referenceImageValue,
197                             double rateValue);
198   //! Destructor.
199   ~Linear1D1DFunction();
200   //@}
201 
202   //! @name Mathematical  methods
203   //@{
204   //! Returns the value of the linear function at point \c domainValue.
205   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
206    * and in affirmative case, it evaluates the function at such point. */
207   double value(double domainValue) const;
208 
209   //! Returns the value of the derivative of the linear function at point \c domainValue.
210   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
211    * and in affirmative case, it returns the value of the derivative at such point (namely
212    * \c m_rateValue). */
213   double deriv(double domainValue) const;
214   //@}
215 protected:
216   using Base1D1DFunction::m_minDomainValue;
217   using Base1D1DFunction::m_maxDomainValue;
218 
219   //! Reference value in the function domain; \f$ x_1 \f$ in  \f$ f(x) = y_1 + m (x - x_1)\f$.
220   double m_referenceDomainValue;
221 
222   //! Reference value in the function image; \f$ y_1 \f$ in  \f$ f(x) = y_1 + m (x - x_1)\f$.
223   double m_referenceImageValue;
224 
225   //! Rate value; \f$ m\f$ in  \f$f(x) = y_1 + m (x - x_1)\f$.
226   double m_rateValue;
227 };
228 
229 //*****************************************************
230 // PiecewiseLinear 1D->1D class
231 //*****************************************************
232 /*! \class PiecewiseLinear1D1DFunction
233     \brief Class for piecewise-linear one-dimensional functions.
234 
235     This class implements piecewise-linear one-dimensional functions.
236 */
237 class PiecewiseLinear1D1DFunction : public Base1D1DFunction {
238 public:
239   //! @name Constructor/Destructor methods
240   //@{
241   //! Default constructor.
242   PiecewiseLinear1D1DFunction(double                     minDomainValue,
243                                      double                     maxDomainValue,
244                                      const std::vector<double>& referenceDomainValues,
245                                      double                     referenceImageValue0,
246                                      const std::vector<double>& rateValues);
247   //! Destructor.
248   ~PiecewiseLinear1D1DFunction();
249   //@}
250 
251   //! @name Mathematical  methods
252   //@{
253   //! Returns the value of the piecewise-linear function at point \c domainValue.
254   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
255    * and in affirmative case, it evaluates the function at such point. */
256   double value(double domainValue) const;
257 
258   //! Returns the value of the derivative of the piecewise-linear function at point \c domainValue.
259   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
260    * and in affirmative case, it returns the value of the derivative at such point. */
261   double deriv(double domainValue) const;
262   //@}
263 protected:
264   using Base1D1DFunction::m_minDomainValue;
265   using Base1D1DFunction::m_maxDomainValue;
266 
267   //! Number of points which will be evaluated.
268   unsigned int        m_numRefValues;
269 
270   //! Reference values in the piecewise-linear function domain; \f$ x_1i \f$ in   \f$ f_i(x) = y_{1i} + m_i (x - x_{1i})\f$, for each \f$ i \f$=1 .. \c m_numRefValues.
271   std::vector<double> m_referenceDomainValues;
272 
273   //! Reference values in the piecewise-linear function image; \f$ y_{1i} \f$ in  \f$ f_i(x) = y_{1i} + m_i (x - x_{1i})\f$, for each \f$ i \f$=1 .. \c m_numRefValues.
274   std::vector<double> m_referenceImageValues;
275 
276   //! Rate value; \f$ m_i \f$ in \f$ f_i(x) = y_{1i} + m_i (x - x_{1i})\f$, for each \f$ i \f$=1 .. \c m_numRefValues.
277   std::vector<double> m_rateValues;
278 };
279 
280 //*****************************************************
281 // Quadratic 1D->1D class
282 //*****************************************************
283 /*! \class Quadratic1D1DFunction
284     \brief Class for one-dimensional quadratic functions.
285 
286     This class implements quadratic one-dimensional functions.
287     A quadratic function, in mathematics, is a polynomial function of the form
288     \f[ f(x)=ax^2+bx+c,\quad a \ne 0.\f]
289     The graph of a quadratic function is a parabola whose axis of symmetry is parallel
290     to the y-axis.
291 */
292 class Quadratic1D1DFunction : public Base1D1DFunction {
293 public:
294     //! @name Constructor/Destructor methods
295   //@{
296   //! Default constructor.
297   /*! Constructs a quadratic function in the interval <c>[minDomainValue,maxDomainValue]</c>,
298    * of the type: \f$ f(x)=ax^2+bx+c \f$.*/
299   Quadratic1D1DFunction(double minDomainValue,
300                                double maxDomainValue,
301                                double a,
302                                double b,
303                                double c);
304   //! Destructor.
305   ~Quadratic1D1DFunction();
306 
307     //! @name Mathematical  methods
308   //@{
309   //! Returns the value of the quadratic function at point \c domainValue.
310   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
311    * and in affirmative case, it evaluates the function at such point, namely:
312    * <c>imageValue = a*domainValue^2 + b*domainValue + c </c>.*/
313   double value(double domainValue) const;
314 
315   //! Returns the value of the derivative of the quadratic function at point \c domainValue.
316   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
317    * and in affirmative case, it returns the value of the derivative at such point (namely
318    * <c>2*a*domainValue + b </c>. */
319   double deriv(double domainValue) const;
320   //@}
321 protected:
322   using Base1D1DFunction::m_minDomainValue;
323   using Base1D1DFunction::m_maxDomainValue;
324 
325   //! 'Quadratic' coefficient of the quadratic function; \f$ a \f$ in \f$ f(x)=ax^2+bx+c \f$.
326   double m_a;
327 
328   //! 'Linear' coefficient of the quadratic function; \f$ b \f$ in \f$ f(x)=ax^2+bx+c \f$.
329   double m_b;
330 
331   //! 'Free' coefficient of the quadratic function; \f$ c \f$ in \f$ f(x)=ax^2+bx+c \f$.
332   double m_c;
333 };
334 
335 //*****************************************************
336 // Sampled 1D->1D class
337 //*****************************************************
338 /*! \class Sampled1D1DFunction
339     \brief Class for one-dimensional sampled functions.
340 
341     This class implements sampled one-dimensional functions.
342     A sample function is one whose values are known only at discrete values of the independent
343     variable; i.e., its values are only known at grid points. There are several QUESO classes
344     which handle grids, all derived from BaseOneDGrid.
345     Sampled functions are usually defined by arrays. */
346 
347 class Sampled1D1DFunction : public Base1D1DFunction {
348 public:
349   //! @name Constructor/Destructor methods
350   //@{
351   //! Default constructor. It should not be called by the user.
352   Sampled1D1DFunction();
353 
354   //! Constructor.
355   /*! When calling this constructor, the user provides the values of the independent variable
356    *(\c domainValues) and their respective values in the image set (independent variable,
357    \c imageValue)*/
358   Sampled1D1DFunction(const std::vector<double>& domainValues,
359                              const std::vector<double>& imageValues);
360 
361   //! Destructor.
362   virtual ~Sampled1D1DFunction();
363   //@}
364 
365   //! @name Mathematical  methods
366   //@{
367   //! Returns the value of the sampled function at point \c domainValue.
368   /*! This function checks if point \c domainValue belongs to the domain of \c this function,
369    * and in affirmative case,  it looks for the image value associated to the domain value
370    * passed to the function. If there isn't any, it calculates a linear approximation for the
371    * image value of \c domainValue, considering its neighbors points in the domain.*/
372   virtual double       value(double domainValue) const;
373 
374   //! <b>Bogus</b>: Derivative of the function.
375   /*! Derivatives are not defined over sampled functions! Thus, this function simply checks if
376    * point \c domainValue belongs to the domain of \c this function, and in affirmative case,
377    * it returns 0.*/
378   double               deriv(double domainValue) const;
379 
380   //! Array of the domain values (values of the independent variable).
381   const   std::vector<double>& domainValues() const;
382 
383   //! Array of the image values (values of the dependent variable).
384   const   std::vector<double>& imageValues () const;
385 
386   //! Checks whether the domain value \c domainValue matches exactly one of the values in the function domain \c domainValues().
387   bool    domainValueMatchesExactly(double domainValue) const;
388 
389   //! Sets the values of the independent (\c domainValues) and dependent (\c imageValues) variables of this sampled function.
390   void    set(const std::vector<double>& domainValues,
391 	      const std::vector<double>& imageValues);
392   //@}
393 
394   //! @name I/O methods
395   //@{
396   //! Prints the values of the function in Matlab/Octave format.
397   virtual void                 printForMatlab(const BaseEnvironment& env, std::ofstream& ofsvar, const std::string& prefixName) const;
398   //@}
399 protected:
400   using Base1D1DFunction::m_minDomainValue;
401   using Base1D1DFunction::m_maxDomainValue;
402 
403   //! Array of the values in the domain of the function (values of the independent variable).
404   std::vector<double> m_domainValues;
405 
406   //! Array of the values in the image of the function (values of the dependent variable).
407   std::vector<double> m_imageValues;
408 };
409 
410 //*****************************************************
411 // 'ScalarTimesFunc' 1D->1D class
412 //*****************************************************
413 /*! \class ScalarTimesFunc1D1DFunction
414     \brief Class for multiplication of a one-dimensional function by a scalar.
415 
416     This class implements the multiplication of a generic one-dimensional function
417     (implemented through any class derived from Base1D1DFunction) by a given
418     scalar.
419 */
420 
421 class ScalarTimesFunc1D1DFunction : public Base1D1DFunction {
422 public:
423   //! @name Constructor/Destructor methods
424   //@{
425   //! Default constructor.
426   ScalarTimesFunc1D1DFunction(double scalar,
427                                      const Base1D1DFunction& func);
428 
429   //! Destructor.
430   ~ScalarTimesFunc1D1DFunction();
431   //@}
432 
433   //! @name Mathematical  methods
434   //@{
435   //! Returns the value of the (one-dimensional) function at point \c domainValue, already multiplied by the scalar (given during construction).
436   /*! This function calls the value() method of the given function and multiplies its return
437    * value by the scalar passed to the constructor. */
438   double value(double domainValue) const;
439 
440   //! TODO: Returns the value of the derivative of the function multiplied by the given scalar at point \c domainValue.
441   /*! \todo Please, implement me! */
442   double deriv(double domainValue) const;
443   //@}
444 protected:
445   using Base1D1DFunction::m_minDomainValue;
446   using Base1D1DFunction::m_maxDomainValue;
447 
448   double m_scalar;
449   const Base1D1DFunction& m_func;
450 };
451 
452 //*****************************************************
453 // 'FuncTimesFunc' 1D->1D class
454 //*****************************************************
455 /*! \class FuncTimesFunc1D1DFunction
456     \brief Class for multiplication of a one-dimensional function by another.
457 
458     This class implements the multiplication of a generic one-dimensional function
459     by another (both functions may be of any type of functions derived from
460     Base1D1DFunction).
461 */
462 class FuncTimesFunc1D1DFunction : public Base1D1DFunction {
463 public:
464   //! @name Constructor/Destructor methods
465   //@{
466   //! Default constructor.
467   FuncTimesFunc1D1DFunction(const Base1D1DFunction& func1,
468                                    const Base1D1DFunction& func2);
469 
470   //! Destructor.
471   ~FuncTimesFunc1D1DFunction();
472   //@}
473 
474   //! @name Mathematical  methods
475   //@{
476   //! Returns the value of the multiplication of a function \c func1 by another function \c func2 at the point \c domainValue, i.e. returnValue = <c>func1(domainValue)*func2(domainValue) </c>.
477   /*! This function calls the value() method of the both functions and multiplies their return
478    * value with one another.*/
479   double value(double domainValue) const;
480 
481   //! TODO: Returns the value of the derivative of the function \c func1 by another function \c func2 at the point \c domainValue.
482   /*! \todo Please, implement me! \n
483    * This method objective requires clarification: are we calculating either (func1*func2)' or func1'*func2'?   */
484   double deriv(double domainValue) const;
485   //@}
486 
487 protected:
488   using Base1D1DFunction::m_minDomainValue;
489   using Base1D1DFunction::m_maxDomainValue;
490 
491   const Base1D1DFunction& m_func1;
492   const Base1D1DFunction& m_func2;
493 };
494 
495 //*****************************************************
496 // 'FuncPlusFunc' 1D->1D class
497 //*****************************************************
498 /*! \class FuncPlusFunc1D1DFunction
499     \brief Class for addition of a one-dimensional function with another.
500 
501     This class implements the addition of two generic one-dimensional functions
502     (both functions may be of any type of functions derived from Base1D1DFunction).
503 */
504 class FuncPlusFunc1D1DFunction : public Base1D1DFunction {
505 public:
506   //! @name Constructor/Destructor methods
507   //@{
508   //! Default constructor.
509   FuncPlusFunc1D1DFunction(const Base1D1DFunction& func1,
510                                   const Base1D1DFunction& func2);
511 
512   //! Destructor.
513   ~FuncPlusFunc1D1DFunction();
514 //@}
515 
516   //! @name Mathematical  methods
517   //@{
518   //! Returns the value of the addition of function \c func1 and \c func2 evaluated at the point \c domainValue, i.e. returnValue = <c>func1(domainValue)+func2(domainValue) </c>.
519   /*! This function calls the value() method of the both functions and adds their return
520    * value with one another.*/
521   double value(double domainValue) const;
522 
523   //! TODO: Returns the value of the derivative of the addition of two functions.
524   /*! \todo Please, implement me! \n*/
525   double deriv(double domainValue) const;
526   //@}
527 protected:
528   using Base1D1DFunction::m_minDomainValue;
529   using Base1D1DFunction::m_maxDomainValue;
530 
531   const Base1D1DFunction& m_func1;
532   const Base1D1DFunction& m_func2;
533 };
534 //*****************************************************
535 // Lagrange Polynomial 1D->1D class
536 //*****************************************************
537 /*! \class LagrangePolynomial1D1DFunction
538     \brief Class for one-dimensional Lagrange polynomials.
539 */
540 
541 /*! The Lagrange interpolating polynomial of a one-dimensional function \f$ f(x) = y \f$
542  * is the polynomial \f$ P(x) \f$ of degree \f$ \leq n-1 \f$ that passes through the
543  * \f$ n \f$ points \f$ (x_1,y_1=f(x_1)),  (x_2,y_2=f(x_2)), ..., (x_n,y_n=f(x_n))\f$,
544  * and is given by:
545  * \f[ P(x)=\sum_{j=1}^n \prod_{k=1; k\not=j}^n \frac{x-x_k}{x_j-x_k}.\f]
546  *
547  * Written explicitly,
548  * \f[ P(x)=\frac{(x-x_2)(x-x_3)...(x-x_n)}{(x_1-x_2)(x_1-x_3)...(x_1-x_n)}y_1+
549             \frac{(x-x_1)(x-x_3)...(x-x_n)}{(x_2-x_1)(x_2-x_3)...(x_2-x_n)}y_2+...+
550             \frac{(x-x_1)(x-x_2)...(x-x_(n-1))}{(x_n-x_1)(x_n-x_2)...(x_n-x_(n-1))}y_n.\f]
551 
552  * In this class, the array <c> std::vector<double>& positionValues </c> stores the points
553  * \f$ x_1, x_2, ... x_n \f$  and the array <c> std::vector<double>* functionValues </c>
554  * stores the points \f$ y_1, y_2, ... y_n \f$ of the Lagrange polynomial.
555  *
556  * \see Archer, Branden and Weisstein, Eric W. "Lagrange Interpolating Polynomial." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html.*/
557 
558 class LagrangePolynomial1D1DFunction : public Base1D1DFunction {
559 public:
560   //! @name Constructor/Destructor methods
561   //@{
562   //! Default constructor.
563   LagrangePolynomial1D1DFunction(const std::vector<double>& positionValues,
564                                         const std::vector<double>* functionValues);
565 
566   //! Destructor.
567   ~LagrangePolynomial1D1DFunction();
568   //@}
569 
570   //! @name Mathematical  methods
571   //@{
572   //! Returns the value of the Lagrange polynomial at point \c domainValue.
573   double value(double domainValue) const;
574 
575   //! TODO: Returns the value of the derivative of the Lagrange polynomial at point \c domainValue.
576   /*! \todo This function checks if point \c domainValue belongs to the domain of \c this function,
577    * and in affirmative case, it returns the value of the derivative at such point. */
578   double deriv(double domainValue) const;
579   //@}
580 protected:
581   using Base1D1DFunction::m_minDomainValue;
582   using Base1D1DFunction::m_maxDomainValue;
583 
584   std::vector<double> m_positionValues;
585   std::vector<double> m_functionValues;
586 };
587 
588 //*****************************************************
589 // Lagrange Basis 1D->1D class
590 //*****************************************************
591 /*! \class LagrangeBasis1D1DFunction
592  *  \brief Class for Lagrange polynomial basis.
593  *
594  * Given a set of \f$ k+1 \f$ data points \f$(x_0, y_0),\ldots,(x_j, y_j),\ldots,(x_k, y_k)\f$
595  * where no two \f$ x_j \f$ are the same, the interpolation polynomial in the Lagrange form is
596  * a linear combination
597  * \f[ L(x) = \sum_{j=0}^{k} y_j \ell_j(x) \f]
598  * of Lagrange basis polynomials
599  * \f[ \ell_j(x) = \prod_{0\le m\le k;\, m\neq j}
600        \frac{x-x_m}{x_j-x_m} = \frac{(x-x_0)}{(x_j-x_0)} \cdots
601        \frac{(x-x_{j-1})}{(x_j-x_{j-1})} \frac{(x-x_{j+1})}{(x_j-x_{j+1})}
602        \cdots \frac{(x-x_k)}{(x_j-x_k)},\f]
603   * where \f$ 0\le j\le k \f$.\n
604   *
605   * This class implements the one-dimensional function (Lagrange basis) \f$ \ell_j(x) \f$.
606   * In this class, the array <c> std::vector<double>& positionValues </c> stores the points
607   * \f$ x_1, x_2, ... x_n \f$  and the  index \f$ j \f$ is stored in \c basisIndex. */
608 
609 class LagrangeBasis1D1DFunction : public Base1D1DFunction {
610 public:
611   //! @name Constructor/Destructor methods
612   //@{
613   //! Default constructor.
614   LagrangeBasis1D1DFunction(const std::vector<double>& positionValues,
615                                    unsigned int basisIndex);
616 
617   //! Destructor.
618   ~LagrangeBasis1D1DFunction();
619   //@}
620 
621   //! @name Mathematical  methods
622   //@{
623   //! Returns the value of the Lagrange basis at point \c domainValue.
624   double value(double domainValue) const;
625 
626   //! TODO: Returns the value of the derivative of the Lagrange basis at point \c domainValue.
627   /*! \todo This function checks if point \c domainValue belongs to the domain of \c this function,
628    * and in affirmative case, it returns the value of the derivative at such point. */
629   double deriv(double domainValue) const;
630   //@}
631 protected:
632   using Base1D1DFunction::m_minDomainValue;
633   using Base1D1DFunction::m_maxDomainValue;
634 
635   std::vector<double> m_positionValues;
636   unsigned int        m_basisIndex;
637 };
638 
639 //----------------------------------------------------------------------
640 //! Calculates the integral of a 2D Gaussian KDE.
641 /*! This function performs numerical integration (via Hermite-Gauss quadrature,
642  * see GaussianHermite1DQuadrature), of a bi-variate Gaussian KDE (refer
643  * to ScalarSequence::subGaussian1dKde(), ArrayOfSequences::gaussianKDE(),
644  * or SequenceOfVectors::subGaussian1dKde(). */
645 
646 template <class T>
647 double SubF1F2Gaussian2dKdeIntegral(const ScalarSequence<T>& scalarSeq1,
648     const ScalarSequence<T>& scalarSeq2, unsigned int initialPos,
649     double scaleValue1, double scaleValue2, const Base1D1DFunction& func1,
650     const Base1D1DFunction& func2, unsigned int quadratureOrder);
651 
652 }  // End namespace QUESO
653 
654 #endif // UQ_1D_1D_FUNCTION_H
655