1 // @HEADER
2 // ************************************************************************
3 //
4 //               Rapid Optimization Library (ROL) Package
5 //                 Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 //              Drew Kouri   (dpkouri@sandia.gov) and
39 //              Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #ifndef ROL_QUADRATUREDEF_HPP
45 #define ROL_QUADRATUREDEF_HPP
46 
47 /** \file   ROL_QuadratureDef.hpp
48     \brief  Definition file for the ROL::Quadrature class.
49     \author Created by D. Kouri and D. Ridzal.
50 */
51 
52 #include "ROL_QuadratureHelpers.hpp"
53 
54 namespace ROL {
55 
56 template <class Real>
Quadrature(int dimension)57 Quadrature<Real>::Quadrature(int dimension) : dimension_(dimension) {
58   points_.clear(); weights_.clear(); accuracy_.clear();
59   accuracy_.resize(dimension);
60 }
61 
62 template <class Real>
addPointAndWeight(std::vector<Real> point,Real weight,int loc)63 void Quadrature<Real>::addPointAndWeight(std::vector<Real> point, Real weight, int loc) {
64   this->points_.insert(std::pair<std::vector<Real>,int>(point,loc));
65   this->weights_.push_back(weight);
66 }
67 
68 template <class Real>
getNumPoints() const69 int Quadrature<Real>::getNumPoints() const {
70   return this->weights_.size();
71 }
72 
73 template <class Real>
getAccuracy(std::vector<int> & accuracy) const74 void Quadrature<Real>::getAccuracy(std::vector<int> & accuracy) const {
75   accuracy = this->accuracy_;
76 }
77 
78 template <class Real>
getCubature(std::vector<std::vector<Real>> & points,std::vector<Real> & weights) const79 void Quadrature<Real>::getCubature(std::vector<std::vector<Real> >& points,
80                                    std::vector<Real>& weights) const {
81   points.resize(this->weights_.size());
82   weights.resize(this->weights_.size());
83   typename std::map<std::vector<Real>,int>::const_iterator it;
84   for (it=(this->points_).begin(); it!=(this->points_).end();it++) {
85     points[it->second] = it->first;
86     weights[it->second] = (this->weights_)[it->second];
87   }
88 }
89 
90 template <class Real>
getDimension() const91 int Quadrature<Real>::getDimension() const {
92   return this->dimension_;
93 }
94 
95 template <class Real>
begin()96 typename std::map<std::vector<Real>,int>::iterator Quadrature<Real>::begin() {
97   return (this->points_).begin();
98 }
99 
100 template <class Real>
end()101 typename std::map<std::vector<Real>,int>::iterator Quadrature<Real>::end() {
102   return (this->points_).end();
103 }
104 
105 template <class Real>
insert(typename std::map<std::vector<Real>,int>::iterator it,std::vector<Real> point,Real weight)106 void Quadrature<Real>::insert(typename std::map<std::vector<Real>,int>::iterator it,
107                               std::vector<Real> point, Real weight) {
108   (this->points_).insert(it,std::pair<std::vector<Real>,int>(point,(int)(this->points_).size()));
109   (this->weights_).push_back(weight);
110 }
111 
112 template <class Real>
getNode(typename std::map<std::vector<Real>,int>::iterator it)113 std::vector<Real> Quadrature<Real>::getNode(typename std::map<std::vector<Real>,int>::iterator it) {
114   return it->first;
115 }
116 
117 template <class Real>
getWeight(int node)118 Real Quadrature<Real>::getWeight(int node) {
119   return (this->weights_)[node];
120 }
121 
122 template <class Real>
getWeight(std::vector<Real> point)123 Real Quadrature<Real>::getWeight(std::vector<Real> point) {
124   return (this->weights_)[(this->points_)[point]];
125 }
126 
127 template <class Real>
update(Real alpha,Quadrature<Real> & rule)128 void Quadrature<Real>::update(Real alpha, Quadrature<Real> &rule) {
129   if ( rule.begin() != rule.end() ) {
130     // Initialize an iterator on std::map<std::vector<Real>,Real>
131     typename std::map<std::vector<Real>,int>::iterator it;
132     typename std::map<std::vector<Real>,int>::iterator it2;
133 
134     // Get location of last point and weight
135     int loc = points_.size();
136 
137     // Intersection of rule1 and rule2 and set difference rule2 \ rule1
138     for ( it2=rule.begin(); it2!=rule.end(); ++it2 ) {
139       it = points_.find(it2->first);
140       if ( it != points_.end() ) {
141         weights_[it->second] += alpha*rule.getWeight(it2->second);
142       }
143       else {
144         points_.insert(std::pair<std::vector<Real>, int>(it2->first,loc));
145         weights_.push_back(alpha*rule.getWeight(it2->second));
146         loc++;
147       }
148     }
149   }
150 }
151 
152 template <class Real>
normalize()153 void Quadrature<Real>::normalize(){
154   Real sum = 0.0;
155   typename std::vector<Real>::iterator it;
156   for (it=weights_.begin(); it!=weights_.end(); it++) {
157     sum += *it;
158   }
159   for (it=weights_.begin(); it!=weights_.end(); it++) {
160     *it /= sum;
161   }
162 }
163 
164 // PRIVATE MEMBER FUNCTION DEFINITIONS
165 template<class Real>
buildInitial(const int dimension,const int maxNumPoints,const std::vector<EQuadrature> & rule1D,const std::vector<EGrowth> & growth1D,const bool isNormalized)166 void Quadrature<Real>::buildInitial(const int dimension,
167                                     const int maxNumPoints,
168                                     const std::vector<EQuadrature> &rule1D,
169                                     const std::vector<EGrowth> &growth1D,
170                                     const bool isNormalized) {
171   accuracy_.clear();
172   accuracy_.resize(dimension);
173   std::vector<int> degree(1);
174   Quadrature<Real> newRule(1);
175   for (int i=0; i<dimension; i++) {
176     // Compute 1D rules
177     int numPoints = growthRule1D(maxNumPoints,growth1D[i],rule1D[i]);
178     Quadrature<Real> rule1(rule1D[i],numPoints,isNormalized);
179     rule1.getAccuracy(degree);
180     accuracy_.push_back(degree[0]);
181     // Build Tensor Rule
182     newRule = kron_prod<Real>(newRule,rule1);
183   }
184   typename std::map<std::vector<Real>,int>::iterator it;
185   int loc = 0;
186   std::vector<Real> node;
187   for (it=newRule.begin(); it!=newRule.end(); it++) {
188     node = it->first;
189     addPointAndWeight(node,newRule.getWeight(node),loc);
190     loc++;
191   }
192 }
193 
194 template<class Real>
buildSGMGA(const int dimension,const int maxLevel,const std::vector<EQuadrature> & rule1D,const std::vector<EGrowth> & growth1D)195 void Quadrature<Real>::buildSGMGA(const int dimension,
196                                   const int maxLevel,
197                                   const std::vector<EQuadrature> &rule1D,
198                                   const std::vector<EGrowth> &growth1D) {
199   std::vector<Real> sparse_point;
200   std::vector<Real> sparse_weight;
201   int point_num = 0;
202 
203   Real tol = std::sqrt(ROL_EPSILON<Real>());
204 
205   GWPointer *gw_compute_points;
206   gw_compute_points  = new GWPointer[dimension];
207   GWPointer *gw_compute_weights;
208   gw_compute_weights = new GWPointer[dimension];
209   GWPointer2 *gw_compute_order;
210   gw_compute_order   = new GWPointer2[dimension];
211 
212   int growth = 0;
213 
214   for (int i=0; i<dimension; ++i) {
215     if(rule1D[i] == QUAD_CLENSHAWCURTIS) {
216       gw_compute_points[i]  = &Quadrature<Real>::ClenshawCurtisPoints;
217       gw_compute_weights[i] = &Quadrature<Real>::ClenshawCurtisWeights;
218       gw_compute_order[i]   = &webbur::level_to_order_exp_cc;
219     }
220     else if (rule1D[i] == QUAD_FEJER2) {
221       gw_compute_points[i]  = &Quadrature<Real>::Fejer2Points;
222       gw_compute_weights[i] = &Quadrature<Real>::Fejer2Weights;
223       gw_compute_order[i]   = &webbur::level_to_order_exp_f2;
224     }
225     else if (rule1D[i] == QUAD_LEGENDRE) {
226       gw_compute_points[i]  = &Quadrature<Real>::LegendrePoints;
227       gw_compute_weights[i] = &Quadrature<Real>::LegendreWeights;
228       gw_compute_order[i]   = &webbur::level_to_order_exp_gauss;
229     }
230     else if (rule1D[i] == QUAD_PATTERSON) {
231       gw_compute_points[i]  = &Quadrature<Real>::PattersonPoints;
232       gw_compute_weights[i] = &Quadrature<Real>::PattersonWeights;
233       gw_compute_order[i]   = &webbur::level_to_order_exp_gp;
234     }
235     else if (rule1D[i] == QUAD_HERMITE) {
236       gw_compute_points[i]  = &Quadrature<Real>::HermitePoints;
237       gw_compute_weights[i] = &Quadrature<Real>::HermiteWeights;
238       gw_compute_order[i]   = &webbur::level_to_order_exp_gauss;
239     }
240     else if (rule1D[i] == QUAD_GENZKEISTER) {
241       gw_compute_points[i]  = &Quadrature<Real>::GenzKeisterPoints;
242       gw_compute_weights[i] = &Quadrature<Real>::GenzKeisterWeights;
243       gw_compute_order[i]   = &webbur::level_to_order_exp_hgk;
244     }
245     else if (rule1D[i] == QUAD_LAGUERRE) {
246       gw_compute_points[i]  = &Quadrature<Real>::LaguerrePoints;
247       gw_compute_weights[i] = &Quadrature<Real>::LaguerreWeights;
248       gw_compute_order[i]   = &webbur::level_to_order_exp_gauss;
249     }
250 
251     if( growth1D[i] == GROWTH_DEFAULT ||
252         growth1D[i] == GROWTH_FULLEXP    ) {
253       growth = 2;
254     }
255     else if ( growth1D[i] == GROWTH_SLOWLIN    ||
256               growth1D[i] == GROWTH_SLOWLINODD ||
257               growth1D[i] == GROWTH_SLOWEXP       ) {
258       growth = 0;
259     }
260     else if ( growth1D[i] == GROWTH_MODLIN ||
261               growth1D[i] == GROWTH_MODEXP    ) {
262       growth = 1;
263     }
264   }
265 
266   std::vector<Real> level_weight(dimension,1);
267   std::vector<int> np(dimension,0);
268   int np_sum = webbur::i4vec_sum(dimension,&np[0]);
269   std::vector<Real> p(np_sum,0);
270 
271   //  Compute necessary data.
272   int point_total_num
273     = webbur::sandia_sgmga_size_total(dimension, &level_weight[0], maxLevel,
274                                       growth, gw_compute_order);
275 
276   point_num
277     = webbur::sandia_sgmga_size(dimension, &level_weight[0], maxLevel,
278                                 gw_compute_points, tol, growth,
279                                 gw_compute_order);
280 
281   std::vector<int> sparse_unique_index(point_total_num);
282   webbur::sandia_sgmga_unique_index(dimension, &level_weight[0], maxLevel,
283                                     gw_compute_points, tol, point_num,
284                                     point_total_num, growth, gw_compute_order,
285                                     &sparse_unique_index[0]);
286 
287   std::vector<int> sparse_order(dimension*point_num,0);
288   std::vector<int> sparse_index(dimension*point_num,0);
289   webbur::sandia_sgmga_index(dimension, &level_weight[0], maxLevel, point_num,
290                              point_total_num, &sparse_unique_index[0], growth,
291                              gw_compute_order, &sparse_order[0],
292                              &sparse_index[0]);
293 
294   //  Compute points and weights.
295   sparse_point.resize(dimension*point_num,0.0);
296   webbur::sandia_sgmga_point(dimension, &level_weight[0], maxLevel,
297                              gw_compute_points, point_num, &sparse_order[0],
298                              &sparse_index[0], growth, gw_compute_order,
299                              &sparse_point[0]);
300 
301   sparse_weight.resize(point_num,0.0);
302   webbur::sandia_sgmga_weight(dimension, &level_weight[0], maxLevel,
303                               gw_compute_weights, point_num, point_total_num,
304                               &sparse_unique_index[0], growth,
305                               gw_compute_order, &sparse_weight[0]);
306 
307   delete [] gw_compute_points;
308   delete [] gw_compute_weights;
309   delete [] gw_compute_order;
310   /***********************************************************************/
311   /******** END BUILD SPARSE GRID USING SANDIA_SGMGA *********************/
312   /***********************************************************************/
313   typename std::map<std::vector<Real>,int>::iterator it;
314   Real weight(0);
315   std::vector<Real> node(dimension);
316   for (int i = 0; i < point_num; ++i) {
317     weight = sparse_weight[i];
318     for (int j = 0; j < dimension; ++j) {
319       node[j] = sparse_point[j+i*dimension];
320     }
321     addPointAndWeight(node,weight,i);
322   }
323 }
324 
325 template<class Real>
ClenshawCurtisPoints(int n,int dim,double x[])326 void Quadrature<Real>::ClenshawCurtisPoints(int n, int dim, double x[]) {
327   webbur::clenshaw_curtis_compute_points(n,x);
328 }
329 
330 template<class Real>
ClenshawCurtisWeights(int n,int dim,double w[])331 void Quadrature<Real>::ClenshawCurtisWeights(int n, int dim, double w[]) {
332   webbur::clenshaw_curtis_compute_weights(n,w);
333 }
334 
335 template<class Real>
Fejer2Points(int n,int dim,double x[])336 void Quadrature<Real>::Fejer2Points(int n, int dim, double x[]) {
337   webbur::fejer2_compute_points(n,x);
338 }
339 
340 template<class Real>
Fejer2Weights(int n,int dim,double w[])341 void Quadrature<Real>::Fejer2Weights(int n, int dim, double w[]) {
342   webbur::fejer2_compute_weights(n,w);
343 }
344 
345 template<class Real>
LegendrePoints(int n,int dim,double x[])346 void Quadrature<Real>::LegendrePoints(int n, int dim, double x[]) {
347   webbur::legendre_compute_points(n,x);
348 }
349 
350 template<class Real>
LegendreWeights(int n,int dim,double w[])351 void Quadrature<Real>::LegendreWeights(int n, int dim, double w[]) {
352   webbur::legendre_compute_weights(n,w);
353 }
354 
355 template<class Real>
PattersonPoints(int n,int dim,double x[])356 void Quadrature<Real>::PattersonPoints(int n, int dim, double x[]) {
357   webbur::patterson_lookup_points(n,x);
358 }
359 
360 template<class Real>
PattersonWeights(int n,int dim,double w[])361 void Quadrature<Real>::PattersonWeights(int n, int dim, double w[]) {
362   webbur::patterson_lookup_weights(n,w);
363 }
364 
365 template<class Real>
HermitePoints(int n,int dim,double x[])366 void Quadrature<Real>::HermitePoints(int n, int dim, double x[]) {
367   webbur::hermite_compute_points(n,x);
368 }
369 
370 template<class Real>
HermiteWeights(int n,int dim,double w[])371 void Quadrature<Real>::HermiteWeights(int n, int dim, double w[]) {
372   webbur::hermite_compute_weights(n,w);
373 }
374 
375 template<class Real>
GenzKeisterPoints(int n,int dim,double x[])376 void Quadrature<Real>::GenzKeisterPoints(int n, int dim, double x[]) {
377   webbur::hermite_genz_keister_lookup_points(n,x);
378 }
379 
380 template<class Real>
GenzKeisterWeights(int n,int dim,double w[])381 void Quadrature<Real>::GenzKeisterWeights(int n, int dim, double w[]) {
382   webbur::hermite_genz_keister_lookup_weights(n,w);
383 }
384 
385 template<class Real>
LaguerrePoints(int n,int dim,double x[])386 void Quadrature<Real>::LaguerrePoints(int n, int dim, double x[]) {
387   webbur::laguerre_compute_points(n,x);
388 }
389 
390 template<class Real>
LaguerreWeights(int n,int dim,double w[])391 void Quadrature<Real>::LaguerreWeights(int n, int dim, double w[]) {
392   webbur::laguerre_compute_weights(n,w);
393 }
394 
395 } // end ROL namespace
396 
397 #endif
398