queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
LinearLagrangeInterpolationSurrogate.C
Go to the documentation of this file.
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 // This class
26 #include <queso/LinearLagrangeInterpolationSurrogate.h>
27 
28 // QUESO
29 #include <queso/GslVector.h>
30 #include <queso/GslMatrix.h>
31 #include <queso/MultiDimensionalIndexing.h>
32 #include <queso/InterpolationSurrogateData.h>
33 
34 namespace QUESO
35 {
36  template<class V, class M>
38  : InterpolationSurrogateBase<V,M>(data)
39  {}
40 
41  template<class V, class M>
42  double LinearLagrangeInterpolationSurrogate<V,M>::evaluate(const V & domainVector) const
43  {
44  /* Populate indices. These are the lower bound global indices for the
45  "element" containing the domainVector */
46  std::vector<unsigned int> indices(this->m_data.dim());
47  this->compute_interval_indices(domainVector, indices);
48 
49  // lower bound coordinates
50  std::vector<double> x_min(this->m_data.dim());
51 
52  // upper bound coordinates
53  std::vector<double> x_max(this->m_data.dim());
54 
55  // Function value at each of the "nodes" for the "element"
56  std::vector<double> values(this->n_coeffs());
57 
58  // Use the indices to populate the x_min, etc. structures
59  this->compute_interval_values( indices, x_min, x_max, values );
60 
61  // Now evaluate the interpolant
62  return this->eval_interpolant( x_min, x_max, values, domainVector );
63  }
64 
65  template<class V, class M>
67  std::vector<unsigned int>& indices) const
68  {
69  queso_assert_equal_to( domainVector.sizeGlobal(), this->m_data.dim() );
70  queso_assert_equal_to( indices.size(), this->m_data.dim() );
71 
72  for( unsigned int d = 0; d < this->m_data.dim(); d++ )
73  {
74  double spacing = this->m_data.spacing(d);
75  indices[d] = std::floor( (domainVector[d] - this->m_data.x_min(d))/spacing );
76 
77  // Index should be less than the number of point along this dimension
78  queso_assert_less( indices[d], this->m_data.get_n_points()[d] );
79  }
80  }
81 
82  template<class V, class M>
83  void LinearLagrangeInterpolationSurrogate<V,M>::compute_interval_values( const std::vector<unsigned int>& indices,
84  std::vector<double>& x_min,
85  std::vector<double>& x_max,
86  std::vector<double>& values ) const
87  {
88  queso_assert_equal_to( x_min.size(), this->m_data.dim() );
89  queso_assert_equal_to( x_max.size(), this->m_data.dim() );
90  queso_assert_equal_to( values.size(), this->n_coeffs() );
91  queso_assert_equal_to( indices.size(), this->m_data.dim() );
92 
93  // First, use the lower bound (global) indices to populate x_min, x_max
94  for( unsigned int d = 0; d < this->m_data.dim(); d++ )
95  {
96  x_min[d] = this->m_data.get_x( d, indices[d] );
97  x_max[d] = x_min[d] + this->m_data.spacing( d );
98  }
99 
100  // Now populate the values.
101  std::vector<unsigned int> local_indices(this->m_data.dim());
102  std::vector<unsigned int> global_indices(this->m_data.dim());
103  for( unsigned int n = 0 ; n < this->n_coeffs(); n++ )
104  {
105  // Figure out local indices (each coordinate is 0 or 1)
106  this->singleToCoords(n,local_indices);
107 
108  /* For each dimension, if the local index is 0, use the lower
109  bound of the global index. If the local index is 1, use the
110  "upper" global index.
111 
112  In 2-D:
113  The lower left corner of the element will have local_indices = [0,0];
114  Then, global_indices = [ indices[0], indices[1] ];
115  The upper right corner will have local_indices = [1,1];
116  Then, global_indices = [ indices[0]+1, indices[1]+1 ]; */
117  for( unsigned int d = 0; d < this->m_data.dim(); d++ )
118  {
119  if( local_indices[d] == 0 )
120  global_indices[d] = indices[d];
121 
122  else if( local_indices[d] == 1 )
123  global_indices[d] = indices[d]+1;
124 
125  // This shouldn't happen
126  else
127  queso_error();
128  }
129 
130  /* Now that we have the global indices for each coordinate,
131  we get the "global" index. This is the index into the global
132  values array */
133  unsigned int global = MultiDimensionalIndexing::coordToGlobal( global_indices, this->m_data.get_n_points() );
134  values[n] = this->m_data.get_value(global);
135  }
136  }
137 
138  template<class V, class M>
139  double LinearLagrangeInterpolationSurrogate<V,M>::eval_interpolant( const std::vector<double>& x_min,
140  const std::vector<double>& x_max,
141  const std::vector<double>& values,
142  const V & domainVector ) const
143  {
144  queso_assert_equal_to( x_min.size(), this->m_data.dim() );
145  queso_assert_equal_to( x_max.size(), this->m_data.dim() );
146  queso_assert_equal_to( values.size(), this->n_coeffs() );
147  queso_assert_equal_to( domainVector.sizeGlobal(), this->m_data.dim() );
148 
149  double interp_value = 0.0;
150 
151  std::vector<unsigned int> indices( this->m_data.dim() );
152 
153  for( unsigned int n = 0; n < this->n_coeffs(); n++ )
154  {
155  this->singleToCoords( n, indices );
156 
157  double shape_fn = this->tensor_product_lagrange( x_min, x_max, indices, domainVector );
158 
159  interp_value += values[n]*shape_fn;
160  }
161 
162  return interp_value;
163  }
164 
165  template<class V, class M>
166  unsigned int LinearLagrangeInterpolationSurrogate<V,M>::coordsToSingle( const std::vector<unsigned int>& indices ) const
167  {
168  unsigned int local_dim = 2;
169  std::vector<unsigned int> n_points( this->m_data.dim(), local_dim );
170 
171  /* We're abusing this function as it does what we need to with local_dim = 2
172  for all entries of n_points. */
173  return MultiDimensionalIndexing::coordToGlobal( indices, n_points );
174  }
175 
176  template<class V, class M>
177  void LinearLagrangeInterpolationSurrogate<V,M>::singleToCoords( unsigned int global, std::vector<unsigned int>& indices ) const
178  {
179  unsigned int local_dim = 2;
180  std::vector<unsigned int> n_points( this->m_data.dim(), local_dim );
181 
182  /* We're abusing this function as it does what we need to with local_dim = 2
183  for all entries of n_points. */
184  return MultiDimensionalIndexing::globalToCoord( global, n_points, indices );
185  }
186 
187  template<class V, class M>
189  const std::vector<double>& x_max,
190  const std::vector<unsigned int>& indices,
191  const V & domainVector ) const
192  {
193  queso_assert_equal_to( x_min.size(), this->m_data.dim() );
194  queso_assert_equal_to( x_max.size(), this->m_data.dim() );
195  queso_assert_equal_to( indices.size(), this->m_data.dim() );
196  queso_assert_equal_to( domainVector.sizeGlobal(), this->m_data.dim() );
197 
198  double value = 1.0;
199 
200  for( unsigned int d = 0; d < this->m_data.dim(); d++ )
201  {
202  value *= this->lagrange_poly( x_min[d], x_max[d], domainVector[d], indices[d] );
203  }
204 
205  return value;
206  }
207 
208  template<class V, class M>
209  double LinearLagrangeInterpolationSurrogate<V,M>::lagrange_poly( double x0, double x1, double x, unsigned int index ) const
210  {
211  // Make sure we're trying to interpolate between the given points
212  queso_assert_less_equal( x, x1 );
213  queso_assert_greater_equal( x, x0 );
214 
215  // Make sure index is either 0 or 1
216  queso_assert( (index == 0) || (index == 1) );
217 
218  double value = 0.0;
219 
220  if( index == 0 )
221  value = (x-x1)/(x0-x1);
222  else
223  value = (x-x0)/(x1-x0);
224 
225  return value;
226  }
227 
228 } // end namespace QUESO
229 
230 // Instantiate
virtual double evaluate(const V &domainVector) const
Evaluates value of the interpolant for the given domainVector.
void compute_interval_values(const std::vector< unsigned int > &indices, std::vector< double > &x_min, std::vector< double > &x_max, std::vector< double > &values) const
Helper function to populate bounding values for the intervals in each dimension.
Base class for interpolation-based surrogates.
double lagrange_poly(double x0, double x1, double x, unsigned int index) const
Evaluate a 1-D Lagrange polynomial at point x.
double eval_interpolant(const std::vector< double > &x_min, const std::vector< double > &x_max, const std::vector< double > &values, const V &domainVector) const
Evaluate multidimensional linear Lagrange interpolant.
void singleToCoords(unsigned int global, std::vector< unsigned int > &indices) const
Inverse of map computed in coordsToSingle.
static void globalToCoord(unsigned int global, const std::vector< unsigned int > &n_points, std::vector< unsigned int > &coord_indices)
Inverse of coordToGlobal map.
static unsigned int coordToGlobal(const std::vector< unsigned int > &coord_indices, const std::vector< unsigned int > &n_points)
Map coordinate indices to a singal global index.
double tensor_product_lagrange(const std::vector< double > &x_min, const std::vector< double > &x_max, const std::vector< unsigned int > &indices, const V &domainVector) const
Compute multidimensional Lagrange polynomial as tensor product of 1D polynomials. ...
void compute_interval_indices(const V &domainVector, std::vector< unsigned int > &indices) const
Helper function to get lower bound indices for each dimension.
unsigned int coordsToSingle(const std::vector< unsigned int > &indices) const
Convert local indices to single &quot;global&quot; index.

Generated on Tue Jun 5 2018 19:48:56 for queso-0.57.1 by  doxygen 1.8.5