DGtal  1.4.beta
LagrangeInterpolation.h
1
17 #pragma once
18
31 #if defined(LagrangeInterpolation_RECURSES)
32 #error Recursive header files inclusion detected in LagrangeInterpolation.h
33 #else // defined(LagrangeInterpolation_RECURSES)
35 #define LagrangeInterpolation_RECURSES
36
37 #if !defined LagrangeInterpolation_h
39 #define LagrangeInterpolation_h
40
42 // Inclusions
43 #include <iostream>
44 #include "DGtal/base/Common.h"
45 #include "DGtal/math/MPolynomial.h"
46 #include "DGtal/kernel/NumberTraits.h"
47 #include "DGtal/kernel/CEuclideanRing.h"
48 #include <vector>
50
51 namespace DGtal
52 {
61  template <typename TEuclideanRing>
63  {
64
65  // ----------------------- public types -----------------------------------
66  public:
68  typedef TEuclideanRing Ring;
70
73  typedef std::size_t Size;
74
75  // ----------------------- Standard services ------------------------------
76  public:
77
82
87  : myDeterminant( NumberTraits< Ring >::ZERO )
88  {}
89
94  LagrangeInterpolation( const std::vector< Ring >& xvalues )
95  {
96  init( xvalues );
97  }
98
103  LagrangeInterpolation( const LagrangeInterpolation & other ) = default;
104
110
117
124
126  inline Size size() const
127  {
128  return myX.size();
129  }
130
132  inline Size degree() const
133  {
134  return size() - 1;
135  }
136
138  inline Ring denominator() const
139  {
140  return myDeterminant;
141  }
142
148  void init( const std::vector< Ring >& xvalues )
149  {
150  myX = xvalues;
151  // Compute Vandermonde determinant
153  for ( Dimension i = 0; (i+1) < size(); ++i )
154  for ( Dimension j = i + 1; j < size(); ++j )
155  myDeterminant *= myX[ j ] - myX[ i ];
156  // compute Lagrange polynomial basis.
157  myLagrangeBasis.clear();
158  myLagrangeBasis.reserve( size() );
159  for ( Dimension j = 0; j < size(); ++j )
160  {
161  Ring c = myDeterminant;
162  for ( Dimension m = 0; m < size(); ++m )
163  if ( m != j )
164  c /= myX[ j ] - myX[ m ];
165  Polynomial P = c; // constant
166  for ( Dimension m = 0; m < size(); ++m )
167  if ( m != j )
168  P *= mmonomial<Ring>( 1 ) - myX[ m ] * mmonomial<Ring>( 0 );
169  myLagrangeBasis.push_back( P );
170  }
171  }
172
182  Polynomial polynomial( const std::vector< Ring >& yvalues )
183  {
184  Polynomial P;
185  if ( yvalues.size() != size() ) return P;
186  for ( Dimension j = 0; j < size(); ++j )
187  P += yvalues[ j ] * myLagrangeBasis[ j ];
188  return P;
189  }
190
191  // ----------------------- Accessors ------------------------------
192  public:
193
199  Polynomial basis( Size i ) const
200  {
201  ASSERT( i < size() );
202  return myLagrangeBasis[ i ];
203  }
204
205  // ----------------------- Interface --------------------------------------
206  public:
207
212  void selfDisplay( std::ostream & that_stream ) const
213  {
214  that_stream << "[LagrangeInterpolation det=" << myDeterminant << std::endl;
215  for ( Size i = 0; i < size(); i++ )
216  that_stream << "l_" << i << "=" << basis( i ) << std::endl;
217  that_stream << "]";
218  }
219
224  bool OK() const
225  {
227  }
228
229  // ------------------------- Datas ----------------------------------------
230  protected:
231
233  std::vector< Ring > myX;
237  std::vector<Polynomial> myLagrangeBasis;
238
239  };
240
247  template <typename TEuclideanRing>
248  std::ostream&
249  operator<<( std::ostream & that_stream,
250  const LagrangeInterpolation< TEuclideanRing > & that_object_to_display )
251  {
252  that_object_to_display.selfDisplay( that_stream );
253  return that_stream;
254  }
255
256
257 } // namespace DGtal
258
259
261 // Includes inline functions/methods if necessary.
262
263 // //
265
266 #endif // !defined LagrangeInterpolation_h
267
268 #undef LagrangeInterpolation_RECURSES
269 #endif // else defined(LagrangeInterpolation_RECURSES)
Aim: This class implements Lagrange basis functions and Lagrange interpolation.
Polynomial polynomial(const std::vector< Ring > &yvalues)
LagrangeInterpolation(const std::vector< Ring > &xvalues)
BOOST_CONCEPT_ASSERT((concepts::CEuclideanRing< Ring >))
LagrangeInterpolation(LagrangeInterpolation &&other)=default
LagrangeInterpolation & operator=(const LagrangeInterpolation &other)=default
LagrangeInterpolation< TEuclideanRing > Self
void selfDisplay(std::ostream &that_stream) const
std::vector< Polynomial > myLagrangeBasis
The Lagrange polynomial basis corresponding to X-values.
Polynomial basis(Size i) const
LagrangeInterpolation & operator=(LagrangeInterpolation &&other)=default
DGtal::MPolynomial< 1, Ring > Polynomial
The monovariate polynomial type.
Ring myDeterminant
The determinant of the Vandermonde matrix corresponding to X-values.
std::vector< Ring > myX
The vector of X-values (abscissa)
void init(const std::vector< Ring > &xvalues)
LagrangeInterpolation(const LagrangeInterpolation &other)=default
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
Definition: MPolynomial.h:965
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::uint32_t Dimension
Definition: Common.h:136
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:564
Aim: Defines the mathematical concept equivalent to a unitary commutative ring with a division operat...