ROL
ROL_SROMGenerator.hpp
Go to the documentation of this file.
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_SROMGENERATOR_HPP
45 #define ROL_SROMGENERATOR_HPP
46 
47 #include "ROL_SampleGenerator.hpp"
48 
49 #include "ROL_Objective.hpp"
50 #include "ROL_BoundConstraint.hpp"
52 
53 #include "ROL_Algorithm.hpp"
54 #include "ROL_BoundConstraint.hpp"
55 
56 #include "ROL_MomentObjective.hpp"
57 #include "ROL_CDFObjective.hpp"
59 #include "ROL_SROMVector.hpp"
60 
61 #include "ROL_StdVector.hpp"
62 
63 namespace ROL {
64 
65 template<class Real>
66 class SROMGenerator : public SampleGenerator<Real> {
67 private:
68  // Parameterlist for optimization
69  Teuchos::ParameterList parlist_;
70  // Vector of distributions (size = dimension of space)
71  std::vector<Teuchos::RCP<Distribution<Real> > > dist_;
72 
73  const int dimension_;
77  bool adaptive_;
78  bool print_;
79 
80  Real ptol_;
81  Real atol_;
82 
84  const AtomVector<Real> &atom) {
85  // Remove points with zero weight
86  std::vector<std::vector<Real> > pts;
87  std::vector<Real> wts;
88  for (int i = 0; i < numMySamples_; i++) {
89  if ( prob.getProbability(i) > ptol_ ) {
90  pts.push_back(*(atom.getAtom(i)));
91  wts.push_back(prob.getProbability(i));
92  }
93  }
94  numMySamples_ = wts.size();
95  // Remove atoms that are within atol of each other
96  Real err = 0.0;
97  std::vector<Real> pt;
98  std::vector<int> ind;
99  for (int i = 0; i < numMySamples_; i++) {
100  pt = pts[i]; ind.clear();
101  for (int j = i+1; j < numMySamples_; j++) {
102  err = 0.0;
103  for (int d = 0; d < dimension_; d++) {
104  err += std::pow(pt[d] - pts[j][d],2);
105  }
106  err = std::sqrt(err);
107  if ( err < atol_ ) {
108  ind.push_back(j);
109  for (int d = 0; d < dimension_; d++) {
110  pts[i][d] += pts[j][d];
111  wts[i] += wts[j];
112  }
113  }
114  }
115  if ( ind.size() > 0 ) {
116  for (int d = 0; d < dimension_; d++) {
117  pts[i][d] /= (Real)(ind.size()+1);
118  }
119  for (int k = ind.size()-1; k >= 0; k--) {
120  pts.erase(pts.begin()+ind[k]);
121  wts.erase(wts.begin()+ind[k]);
122  }
123  }
124  numMySamples_ = wts.size();
125  }
126  // Renormalize weights
127  Real psum = 0.0, sum = 0.0;
128  for (int i = 0; i < numMySamples_; i++) {
129  psum += wts[i];
130  }
131  SampleGenerator<Real>::sumAll(&psum,&sum,1);
132  for (int i = 0; i < numMySamples_; i++) {
133  wts[i] /= sum;
134  }
135  // Set points and weights
138  }
139 
140 public:
141 
142  SROMGenerator(Teuchos::ParameterList &parlist,
143  const Teuchos::RCP<BatchManager<Real> > &bman,
144  const std::vector<Teuchos::RCP<Distribution<Real> > > &dist)
145  : SampleGenerator<Real>(bman), parlist_(parlist), dist_(dist),
146  dimension_(dist.size()) {
147  // Get SROM sublist
148  Teuchos::ParameterList list = parlist.sublist("SOL").sublist("Sample Generator").sublist("SROM");
149  numSamples_ = list.get("Number of Samples",50);
150  adaptive_ = list.get("Adaptive Sampling",false);
151  numNewSamples_ = list.get("Number of New Samples Per Adaptation",0);
152  print_ = list.get("Output to Screen",false);
153  ptol_ = list.get("Probability Tolerance",1.e2*std::sqrt(ROL_EPSILON));
154  atol_ = list.get("Atom Tolerance",1.e2*std::sqrt(ROL_EPSILON));
155  print_ *= !SampleGenerator<Real>::batchID();
156  // Compute batch local number of samples
157  int rank = (int)SampleGenerator<Real>::batchID();
158  int nProc = (int)SampleGenerator<Real>::numBatches();
159  int frac = numSamples_ / nProc;
160  int rem = numSamples_ % nProc;
161  numMySamples_ = frac + ((rank < rem) ? 1 : 0);
162  // Initialize vectors
163  Teuchos::RCP<ProbabilityVector<Real> > prob, prob_lo, prob_hi, prob_eq;
164  Teuchos::RCP<AtomVector<Real> > atom, atom_lo, atom_hi, atom_eq;
165  Teuchos::RCP<Vector<Real> > x, x_lo, x_hi, x_eq;
166  initialize_vectors(prob,prob_lo,prob_hi,prob_eq,atom,atom_lo,atom_hi,atom_eq,x,x_lo,x_hi,x_eq,bman);
167  Teuchos::RCP<Vector<Real> > l
168  = Teuchos::rcp(new StdVector<Real>(Teuchos::rcp(new std::vector<Real>(1,0.))));
169  bool optProb = false, optAtom = true;
170  for ( int i = 0; i < 2; i++ ) {
171  if ( i == 0 ) { optProb = false; optAtom = true; }
172  if ( i == 1 ) { optProb = true; optAtom = true; }
173  // Initialize objective function
174  std::vector<Teuchos::RCP<Objective<Real> > > obj_vec;
175  Teuchos::RCP<Objective<Real> > obj;
176  initialize_objective(obj_vec,obj,dist,bman,optProb,optAtom,list);
177  // Initialize constraints
178  Teuchos::RCP<BoundConstraint<Real> > bnd
179  = Teuchos::rcp(new BoundConstraint<Real>(x_lo,x_hi));
180  Teuchos::RCP<EqualityConstraint<Real> > con
181  = Teuchos::rcp(new ScalarLinearEqualityConstraint<Real>(x_eq,1.0));
182  // Test objective and constraints
183  if ( print_ ) { std::cout << "\nCheck derivatives of CDFObjective\n"; }
184  check_objective(*x,obj_vec[0],bman,optProb,optAtom);
185  if ( print_ ) { std::cout << "\nCheck derivatives of MomentObjective\n"; }
186  check_objective(*x,obj_vec[1],bman,optProb,optAtom);
187  if ( print_ ) { std::cout << "\nCheck derivatives of LinearCombinationObjective\n"; }
188  check_objective(*x,obj,bman,optProb,optAtom);
189  if ( print_ && optProb ) { std::cout << "\nCheck ScalarLinearEqualityConstraint\n"; }
190  check_constraint(*x,con,bman,optProb);
191  // Solve optimization problems to sample
192  Teuchos::RCP<Algorithm<Real> > algo;
193  initialize_optimizer(algo,list,optProb);
194  if ( optProb ) {
195  Teuchos::RCP<Teuchos::ParameterList> plist = Teuchos::rcp(&list,false);
196  ROL::OptimizationProblem<Real> optProblem(obj,x,bnd,con,l,plist);
197  algo->run(optProblem,print_);
198  }
199  else {
200  algo->run(*x,*obj,*bnd,print_);
201  }
202  }
203  // Prune samples with zero weight and set samples/weights
204  pruneSamples(*prob,*atom);
205  }
206 
207  void refine(void) {}
208 
209 private:
210 
211  void get_scaling_vectors(std::vector<Real> &typw, std::vector<Real> &typx) const {
212  typw.clear(); typx.clear();
213  typw.resize(numMySamples_,(Real)(numSamples_*numSamples_));
214  typx.resize(numMySamples_*dimension_,0.);
215  Real mean = 1., var = 1.;
216  for (int j = 0; j < dimension_; j++) {
217  mean = std::abs(dist_[j]->moment(1));
218  var = dist_[j]->moment(2) - mean*mean;
219  mean = ((mean > ROL_EPSILON) ? mean : std::sqrt(var));
220  for (int i = 0; i < numMySamples_; i++) {
221  typx[i*dimension_ + j] = 1./(mean*mean);
222  }
223  }
224  }
225 
226  void initialize_vectors(Teuchos::RCP<ProbabilityVector<Real> > &prob,
227  Teuchos::RCP<ProbabilityVector<Real> > &prob_lo,
228  Teuchos::RCP<ProbabilityVector<Real> > &prob_hi,
229  Teuchos::RCP<ProbabilityVector<Real> > &prob_eq,
230  Teuchos::RCP<AtomVector<Real> > &atom,
231  Teuchos::RCP<AtomVector<Real> > &atom_lo,
232  Teuchos::RCP<AtomVector<Real> > &atom_hi,
233  Teuchos::RCP<AtomVector<Real> > &atom_eq,
234  Teuchos::RCP<Vector<Real> > &vec,
235  Teuchos::RCP<Vector<Real> > &vec_lo,
236  Teuchos::RCP<Vector<Real> > &vec_hi,
237  Teuchos::RCP<Vector<Real> > &vec_eq,
238  const Teuchos::RCP<BatchManager<Real> > &bman) const {
239  // Compute scaling for probability and atom vectors
240  std::vector<Real> typx, typw;
241  get_scaling_vectors(typw,typx);
242  // Compute initial guess and bounds for probability and atom vectors
243  std::vector<Real> pt(dimension_*numMySamples_,0.), wt(numMySamples_,1./(Real)numSamples_);
244  std::vector<Real> pt_lo(dimension_*numMySamples_,0.), pt_hi(dimension_*numMySamples_,0.);
245  std::vector<Real> wt_lo(numMySamples_,0.), wt_hi(numMySamples_,1.);
246  std::vector<Real> pt_eq(dimension_*numMySamples_,0.), wt_eq(numMySamples_,1.);
247  Real lo = 0., hi = 0.;
248  srand(12345*SampleGenerator<Real>::batchID());
249  for ( int j = 0; j < dimension_; j++) {
250  lo = dist_[j]->lowerBound();
251  hi = dist_[j]->upperBound();
252  for (int i = 0; i < numMySamples_; i++) {
253  pt[i*dimension_ + j] = dist_[j]->invertCDF((Real)rand()/(Real)RAND_MAX);
254  //pt[i*dimension_ + j] = dist_[j]->invertCDF(0);;
255  pt_lo[i*dimension_ + j] = lo;
256  pt_hi[i*dimension_ + j] = hi;
257  }
258  }
259  // Build probability, atom, and SROM vectors
260  prob = Teuchos::rcp(new PrimalProbabilityVector<Real>(
261  Teuchos::rcp(new std::vector<Real>(wt)),
262  Teuchos::rcp(new std::vector<Real>(typw)),bman));
263  atom = Teuchos::rcp(new PrimalAtomVector<Real>(
264  Teuchos::rcp(new std::vector<Real>(pt)),numMySamples_,dimension_,
265  Teuchos::rcp(new std::vector<Real>(typx)),bman));
266  vec = Teuchos::rcp(new SROMVector<Real>(prob,atom));
267  // Lower and upper bounds on Probability Vector
268  prob_lo = Teuchos::rcp(new PrimalProbabilityVector<Real>(
269  Teuchos::rcp(new std::vector<Real>(wt_lo)),
270  Teuchos::rcp(new std::vector<Real>(typw)),bman));
271  prob_hi = Teuchos::rcp(new PrimalProbabilityVector<Real>(
272  Teuchos::rcp(new std::vector<Real>(wt_hi)),
273  Teuchos::rcp(new std::vector<Real>(typw)),bman));
274  // Lower and upper bounds on Atom Vector
275  atom_lo = Teuchos::rcp(new PrimalAtomVector<Real>(
276  Teuchos::rcp(new std::vector<Real>(pt_lo)),numMySamples_,dimension_,
277  Teuchos::rcp(new std::vector<Real>(typx)),bman));
278  atom_hi = Teuchos::rcp(new PrimalAtomVector<Real>(
279  Teuchos::rcp(new std::vector<Real>(pt_hi)),numMySamples_,dimension_,
280  Teuchos::rcp(new std::vector<Real>(typx)),bman));
281  // Lower and upper bounds on SROM Vector
282  vec_lo = Teuchos::rcp(new SROMVector<Real>(prob_lo,atom_lo));
283  vec_hi = Teuchos::rcp(new SROMVector<Real>(prob_hi,atom_hi));
284  // Equality constraint vectors
285  prob_eq = Teuchos::rcp(new DualProbabilityVector<Real>(
286  Teuchos::rcp(new std::vector<Real>(wt_eq)),
287  Teuchos::rcp(new std::vector<Real>(typw)), bman));
288  atom_eq = Teuchos::rcp(new DualAtomVector<Real>(
289  Teuchos::rcp(new std::vector<Real>(pt_eq)),numMySamples_,dimension_,
290  Teuchos::rcp(new std::vector<Real>(typx)), bman));
291  vec_eq = Teuchos::rcp(new SROMVector<Real>(prob_eq,atom_eq));
292  }
293 
294  void initialize_objective(std::vector<Teuchos::RCP<Objective<Real> > > &obj_vec,
295  Teuchos::RCP<Objective<Real> > &obj,
296  const std::vector<Teuchos::RCP<Distribution<Real> > > &dist,
297  const Teuchos::RCP<BatchManager<Real> > &bman,
298  const bool optProb, const bool optAtom,
299  Teuchos::ParameterList &list) const {
300  // Build CDF objective function
301  Real scale = list.get("CDF Smoothing Parameter",1.e-2);
302  obj_vec.push_back(Teuchos::rcp(new CDFObjective<Real>(dist,bman,scale,optProb,optAtom)));
303  // Build moment matching objective function
304  Teuchos::Array<int> tmp_order
305  = Teuchos::getArrayFromStringParameter<int>(list,"Moments");
306  std::vector<int> order(tmp_order.size(),0);
307  for ( int i = 0; i < tmp_order.size(); i++) {
308  order[i] = static_cast<int>(tmp_order[i]);
309  }
310  obj_vec.push_back(Teuchos::rcp(new MomentObjective<Real>(dist,order,bman,optProb,optAtom)));
311  // Build linear combination objective function
312  Teuchos::Array<Real> tmp_coeff
313  = Teuchos::getArrayFromStringParameter<Real>(list,"Coefficients");
314  std::vector<Real> coeff(2,0.);
315  coeff[0] = tmp_coeff[0]; coeff[1] = tmp_coeff[1];
316  obj = Teuchos::rcp(new LinearCombinationObjective<Real>(coeff,obj_vec));
317  }
318 
319  void initialize_optimizer(Teuchos::RCP<Algorithm<Real> > &algo,
320  Teuchos::ParameterList &parlist,
321  const bool optProb) const {
322  std::string type = parlist.sublist("Step").get("Type","Trust Region");
323  if ( optProb ) {
324  if ( type == "Moreau Yosida" ) {
325  algo = Teuchos::rcp(new Algorithm<Real>("Moreau-Yosida Penalty",parlist,false));
326  }
327  else if ( type == "Augmented Lagrangian" ) {
328  algo = Teuchos::rcp(new Algorithm<Real>("Augmented Lagrangian",parlist,false));
329  }
330  else {
331  algo = Teuchos::rcp(new Algorithm<Real>("Interior Point",parlist,false));
332  }
333  }
334  else {
335  algo = Teuchos::rcp(new Algorithm<Real>("Trust Region",parlist,false));
336  }
337  }
338 
340  const Teuchos::RCP<Objective<Real> > &obj,
341  const Teuchos::RCP<BatchManager<Real> > &bman,
342  const bool optProb, const bool optAtom) {
343  // Get scaling for probability and atom vectors
344  std::vector<Real> typx, typw;
345  get_scaling_vectors(typw,typx);
346  // Set random direction
347  std::vector<Real> pt(dimension_*numMySamples_,0.), wt(numMySamples_,0.);
348  for (int i = 0; i < numMySamples_; i++) {
349  wt[i] = (optProb ? (Real)rand()/(Real)RAND_MAX : 0);
350  for ( int j = 0; j < dimension_; j++) {
351  pt[i*dimension_ + j] = (optAtom ? dist_[j]->invertCDF((Real)rand()/(Real)RAND_MAX) : 0);
352  }
353  }
354  Teuchos::RCP<ProbabilityVector<Real> > dprob
355  = Teuchos::rcp(new PrimalProbabilityVector<Real>(
356  Teuchos::rcp(new std::vector<Real>(wt)),
357  Teuchos::rcp(new std::vector<Real>(typw)),bman));
358  Teuchos::RCP<AtomVector<Real> > datom
359  = Teuchos::rcp(new PrimalAtomVector<Real>(
360  Teuchos::rcp(new std::vector<Real>(pt)),numMySamples_,dimension_,
361  Teuchos::rcp(new std::vector<Real>(typx)),bman));
362  SROMVector<Real> d = SROMVector<Real>(dprob,datom);
363  // Check derivatives
364  obj->checkGradient(x,d,print_);
365  obj->checkHessVec(x,d,print_);
366  }
367 
369  const Teuchos::RCP<EqualityConstraint<Real> > &con,
370  const Teuchos::RCP<BatchManager<Real> > &bman,
371  const bool optProb) {
372  if ( optProb ) {
373  StdVector<Real> c(Teuchos::rcp(new std::vector<Real>(1,1.0)));
374  // Get scaling for probability and atom vectors
375  std::vector<Real> typx, typw;
376  get_scaling_vectors(typw,typx);
377  // Set random direction
378  std::vector<Real> pt(dimension_*numMySamples_,0.), wt(numMySamples_,0.);
379  for (int i = 0; i < numMySamples_; i++) {
380  wt[i] = (Real)rand()/(Real)RAND_MAX;
381  for ( int j = 0; j < dimension_; j++) {
382  pt[i*dimension_ + j] = dist_[j]->invertCDF((Real)rand()/(Real)RAND_MAX);
383  }
384  }
385  Teuchos::RCP<ProbabilityVector<Real> > dprob
386  = Teuchos::rcp(new PrimalProbabilityVector<Real>(
387  Teuchos::rcp(new std::vector<Real>(wt)),
388  Teuchos::rcp(new std::vector<Real>(typw)),bman));
389  Teuchos::RCP<AtomVector<Real> > datom
390  = Teuchos::rcp(new PrimalAtomVector<Real>(
391  Teuchos::rcp(new std::vector<Real>(pt)),numMySamples_,dimension_,
392  Teuchos::rcp(new std::vector<Real>(typx)),bman));
393  SROMVector<Real> d = SROMVector<Real>(dprob,datom);
394  // Check derivatives
395  con->checkApplyJacobian(x,d,c,print_);
396  con->checkAdjointConsistencyJacobian(c,d,x,print_);
397  }
398  }
399 };
400 
401 }
402 
403 #endif
Provides the interface to evaluate objective functions.
Provides the std::vector implementation of the ROL::Vector interface.
const Real getProbability(const int i) const
This equality constraint defines an affine hyperplane.
void check_objective(const Vector< Real > &x, const Teuchos::RCP< Objective< Real > > &obj, const Teuchos::RCP< BatchManager< Real > > &bman, const bool optProb, const bool optAtom)
void check_constraint(const Vector< Real > &x, const Teuchos::RCP< EqualityConstraint< Real > > &con, const Teuchos::RCP< BatchManager< Real > > &bman, const bool optProb)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void sumAll(Real *input, Real *output, int dim) const
Defines the equality constraint operator interface.
Teuchos::RCP< const std::vector< Real > > getAtom(const int i) const
void initialize_optimizer(Teuchos::RCP< Algorithm< Real > > &algo, Teuchos::ParameterList &parlist, const bool optProb) const
void get_scaling_vectors(std::vector< Real > &typw, std::vector< Real > &typx) const
Provides an interface to run optimization algorithms.
Teuchos::ParameterList parlist_
Provides the std::vector implementation of the ROL::Vector interface.
Provides the std::vector implementation of the ROL::Vector interface.
void initialize_vectors(Teuchos::RCP< ProbabilityVector< Real > > &prob, Teuchos::RCP< ProbabilityVector< Real > > &prob_lo, Teuchos::RCP< ProbabilityVector< Real > > &prob_hi, Teuchos::RCP< ProbabilityVector< Real > > &prob_eq, Teuchos::RCP< AtomVector< Real > > &atom, Teuchos::RCP< AtomVector< Real > > &atom_lo, Teuchos::RCP< AtomVector< Real > > &atom_hi, Teuchos::RCP< AtomVector< Real > > &atom_eq, Teuchos::RCP< Vector< Real > > &vec, Teuchos::RCP< Vector< Real > > &vec_lo, Teuchos::RCP< Vector< Real > > &vec_hi, Teuchos::RCP< Vector< Real > > &vec_eq, const Teuchos::RCP< BatchManager< Real > > &bman) const
SROMGenerator(Teuchos::ParameterList &parlist, const Teuchos::RCP< BatchManager< Real > > &bman, const std::vector< Teuchos::RCP< Distribution< Real > > > &dist)
Provides the interface to apply upper and lower bound constraints.
void initialize_objective(std::vector< Teuchos::RCP< Objective< Real > > > &obj_vec, Teuchos::RCP< Objective< Real > > &obj, const std::vector< Teuchos::RCP< Distribution< Real > > > &dist, const Teuchos::RCP< BatchManager< Real > > &bman, const bool optProb, const bool optAtom, Teuchos::ParameterList &list) const
std::vector< Teuchos::RCP< Distribution< Real > > > dist_
void setPoints(std::vector< std::vector< Real > > &p)
void setWeights(std::vector< Real > &w)
void pruneSamples(const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom)
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118