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
51#include "ROL_CDFObjective.hpp"
53#include "ROL_SROMVector.hpp"
54#include "ROL_StdVector.hpp"
56#include "ROL_Bounds.hpp"
57
58namespace ROL {
59
60template<class Real>
61class SROMGenerator : public SampleGenerator<Real> {
62private:
63 // Parameterlist for optimization
64 ROL::ParameterList parlist_;
65 // Vector of distributions (size = dimension of space)
66 std::vector<Ptr<Distribution<Real>>> dist_;
67
68 const int dimension_;
73
74 Real ptol_;
75 Real atol_;
76
78 const AtomVector<Real> &atom) {
79 // Remove points with zero weight
80 std::vector<std::vector<Real>> pts;
81 std::vector<Real> wts;
82 for (int i = 0; i < numMySamples_; i++) {
83 if ( prob.getProbability(i) > ptol_ ) {
84 pts.push_back(*(atom.getAtom(i)));
85 wts.push_back(prob.getProbability(i));
86 }
87 }
88 numMySamples_ = wts.size();
89 // Remove atoms that are within atol of each other
90 Real err = 0.0;
91 std::vector<Real> pt;
92 std::vector<int> ind;
93 for (int i = 0; i < numMySamples_; i++) {
94 pt = pts[i]; ind.clear();
95 for (int j = i+1; j < numMySamples_; j++) {
96 err = 0.0;
97 for (int d = 0; d < dimension_; d++) {
98 err += std::pow(pt[d] - pts[j][d],2);
99 }
100 err = std::sqrt(err);
101 if ( err < atol_ ) {
102 ind.push_back(j);
103 for (int d = 0; d < dimension_; d++) {
104 pts[i][d] += pts[j][d];
105 wts[i] += wts[j];
106 }
107 }
108 }
109 if ( ind.size() > 0 ) {
110 for (int d = 0; d < dimension_; d++) {
111 pts[i][d] /= (Real)(ind.size()+1);
112 }
113 for (int k = ind.size()-1; k >= 0; k--) {
114 pts.erase(pts.begin()+ind[k]);
115 wts.erase(wts.begin()+ind[k]);
116 }
117 }
118 numMySamples_ = wts.size();
119 }
120 // Renormalize weights
121 Real psum = 0.0, sum = 0.0;
122 for (int i = 0; i < numMySamples_; i++) {
123 psum += wts[i];
124 }
125 SampleGenerator<Real>::sumAll(&psum,&sum,1);
126 for (int i = 0; i < numMySamples_; i++) {
127 wts[i] /= sum;
128 }
129 // Set points and weights
132 }
133
134public:
135
136 SROMGenerator(ROL::ParameterList &parlist,
137 const Ptr<BatchManager<Real>> &bman,
138 const std::vector<Ptr<Distribution<Real>>> &dist,
139 std::ostream &outStream = std::cout)
140 : SampleGenerator<Real>(bman), parlist_(parlist), dist_(dist),
141 dimension_(dist.size()) {
142 // Get SROM sublist
143 ROL::ParameterList &list = parlist.sublist("SOL").sublist("Sample Generator").sublist("SROM");
144 numSamples_ = list.get("Number of Samples",50);
145 adaptive_ = list.get("Adaptive Sampling",false);
146 numNewSamples_ = list.get("Number of New Samples Per Adaptation",0);
147 ptol_ = list.get("Probability Tolerance",1.e2*std::sqrt(ROL_EPSILON<Real>()));
148 atol_ = list.get("Atom Tolerance",1.e2*std::sqrt(ROL_EPSILON<Real>()));
149 bool presolve = list.get("Presolve for Atom Locations",false);
150 // Compute batch local number of samples
151 int rank = (int)SampleGenerator<Real>::batchID();
152 int nProc = (int)SampleGenerator<Real>::numBatches();
153 int frac = numSamples_ / nProc;
154 int rem = numSamples_ % nProc;
155 numMySamples_ = frac + ((rank < rem) ? 1 : 0);
156 // Initialize vectors
157 Ptr<ProbabilityVector<Real>> prob, prob_lo, prob_hi, prob_eq;
158 Ptr<AtomVector<Real>> atom, atom_lo, atom_hi, atom_eq;
159 Ptr<Vector<Real>> x, x_lo, x_hi, x_eq;
160 initialize_vectors(prob,prob_lo,prob_hi,prob_eq,atom,atom_lo,atom_hi,atom_eq,x,x_lo,x_hi,x_eq,bman);
161 Ptr<Vector<Real>> l = makePtr<SingletonVector<Real>>(0.0);
162 // Initialize constraints
163 Ptr<BoundConstraint<Real>> bnd = makePtr<Bounds<Real>>(x_lo,x_hi);
164 Ptr<Constraint<Real>> con = makePtr<ScalarLinearConstraint<Real>>(x_eq,1.0);
165 if (presolve) { // Optimize over atom locations only
166 ROL::ParameterList pslist(list);
167 pslist.sublist("Step").set("Type","Trust Region");
168 Ptr<Objective<Real>> obj = initialize_objective(dist,bman,false,true,pslist);
169 OptimizationProblem<Real> optProblem(obj,x,bnd);
170 OptimizationSolver<Real> optSolver(optProblem, pslist);
171 optSolver.solve(outStream);
172 }
173 // Optimization over atom locations and probabilities
174 Ptr<Objective<Real>> obj = initialize_objective(dist,bman,true,true,list);
175 OptimizationProblem<Real> optProblem(obj,x,bnd,con,l);
176 optProblem.check(outStream);
177 OptimizationSolver<Real> optSolver(optProblem, list);
178 optSolver.solve(outStream);
179 // Prune samples with zero weight and set samples/weights
180 pruneSamples(*prob,*atom);
181 }
182
183 void refine(void) {}
184
185private:
186
187 void get_scaling_vectors(std::vector<Real> &typw, std::vector<Real> &typx) const {
188 typw.clear(); typx.clear();
189 typw.resize(numMySamples_,(Real)(numSamples_*numSamples_));
190 typx.resize(numMySamples_*dimension_,0);
191 Real mean = 1, var = 1, one(1);
192 for (int j = 0; j < dimension_; j++) {
193 mean = std::abs(dist_[j]->moment(1));
194 var = dist_[j]->moment(2) - mean*mean;
195 mean = ((mean > ROL_EPSILON<Real>()) ? mean : std::sqrt(var));
196 mean = ((mean > ROL_EPSILON<Real>()) ? mean : one);
197 for (int i = 0; i < numMySamples_; i++) {
198 typx[i*dimension_ + j] = one/(mean*mean);
199 }
200 }
201 }
202
204 Ptr<ProbabilityVector<Real>> &prob_lo,
205 Ptr<ProbabilityVector<Real>> &prob_hi,
206 Ptr<ProbabilityVector<Real>> &prob_eq,
207 Ptr<AtomVector<Real>> &atom,
208 Ptr<AtomVector<Real>> &atom_lo,
209 Ptr<AtomVector<Real>> &atom_hi,
210 Ptr<AtomVector<Real>> &atom_eq,
211 Ptr<Vector<Real>> &vec,
212 Ptr<Vector<Real>> &vec_lo,
213 Ptr<Vector<Real>> &vec_hi,
214 Ptr<Vector<Real>> &vec_eq,
215 const Ptr<BatchManager<Real>> &bman) const {
216 // Compute scaling for probability and atom vectors
217 std::vector<Real> typx, typw;
218 get_scaling_vectors(typw,typx);
219 // Compute initial guess and bounds for probability and atom vectors
220 std::vector<Real> pt(dimension_*numMySamples_,0.), wt(numMySamples_,1./(Real)numSamples_);
221 std::vector<Real> pt_lo(dimension_*numMySamples_,0.), pt_hi(dimension_*numMySamples_,0.);
222 std::vector<Real> wt_lo(numMySamples_,0.), wt_hi(numMySamples_,1.);
223 std::vector<Real> pt_eq(dimension_*numMySamples_,0.), wt_eq(numMySamples_,1.);
224 Real lo = 0., hi = 0.;
225 srand(12345*SampleGenerator<Real>::batchID());
226 for ( int j = 0; j < dimension_; j++) {
227 lo = dist_[j]->lowerBound();
228 hi = dist_[j]->upperBound();
229 for (int i = 0; i < numMySamples_; i++) {
230 pt[i*dimension_ + j] = dist_[j]->invertCDF((Real)rand()/(Real)RAND_MAX);
231 //pt[i*dimension_ + j] = dist_[j]->invertCDF(0);
232 pt_lo[i*dimension_ + j] = lo;
233 pt_hi[i*dimension_ + j] = hi;
234 }
235 }
236 // Build probability, atom, and SROM vectors
237 prob = makePtr<PrimalProbabilityVector<Real>>(
238 makePtr<std::vector<Real>>(wt),bman,
239 makePtr<std::vector<Real>>(typw));
240 atom = makePtr<PrimalAtomVector<Real>>(
241 makePtr<std::vector<Real>>(pt),bman,numMySamples_,dimension_,
242 makePtr<std::vector<Real>>(typx));
243 vec = makePtr<SROMVector<Real>>(prob,atom);
244 // Lower and upper bounds on Probability Vector
245 prob_lo = makePtr<PrimalProbabilityVector<Real>>(
246 makePtr<std::vector<Real>>(wt_lo),bman,
247 makePtr<std::vector<Real>>(typw));
248 prob_hi = makePtr<PrimalProbabilityVector<Real>>(
249 makePtr<std::vector<Real>>(wt_hi),bman,
250 makePtr<std::vector<Real>>(typw));
251 // Lower and upper bounds on Atom Vector
252 atom_lo = makePtr<PrimalAtomVector<Real>>(
253 makePtr<std::vector<Real>>(pt_lo),bman,numMySamples_,dimension_,
254 makePtr<std::vector<Real>>(typx));
255 atom_hi = makePtr<PrimalAtomVector<Real>>(
256 makePtr<std::vector<Real>>(pt_hi),bman,numMySamples_,dimension_,
257 makePtr<std::vector<Real>>(typx));
258 // Lower and upper bounds on SROM Vector
259 vec_lo = makePtr<SROMVector<Real>>(prob_lo,atom_lo);
260 vec_hi = makePtr<SROMVector<Real>>(prob_hi,atom_hi);
261 // Constraint vectors
262 prob_eq = makePtr<DualProbabilityVector<Real>>(
263 makePtr<std::vector<Real>>(wt_eq),bman,
264 makePtr<std::vector<Real>>(typw));
265 atom_eq = makePtr<DualAtomVector<Real>>(
266 makePtr<std::vector<Real>>(pt_eq),bman,numMySamples_,dimension_,
267 makePtr<std::vector<Real>>(typx));
268 vec_eq = makePtr<SROMVector<Real>>(prob_eq,atom_eq);
269 }
270
271 Ptr<Objective<Real>> initialize_objective(const std::vector<Ptr<Distribution<Real>>> &dist,
272 const Ptr<BatchManager<Real>> &bman,
273 const bool optProb, const bool optAtom,
274 ROL::ParameterList &list) const {
275 std::vector<Ptr<Objective<Real>>> obj_vec;
276 // Build CDF objective function
277 Real scale = list.get("CDF Smoothing Parameter",1.e-2);
278 obj_vec.push_back(makePtr<CDFObjective<Real>>(dist,bman,scale,optProb,optAtom));
279 // Build moment matching objective function
280 std::vector<int> tmp_order
281 = ROL::getArrayFromStringParameter<int>(list,"Moments");
282 std::vector<int> order(tmp_order.size(),0);
283 for (unsigned int i = 0; i < tmp_order.size(); i++) {
284 order[i] = static_cast<int>(tmp_order[i]);
285 }
286 obj_vec.push_back(makePtr<MomentObjective<Real>>(dist,order,bman,optProb,optAtom));
287 // Build linear combination objective function
288 std::vector<Real> tmp_coeff
289 = ROL::getArrayFromStringParameter<Real>(list,"Coefficients");
290 std::vector<Real> coeff(2,0.);
291 coeff[0] = tmp_coeff[0]; coeff[1] = tmp_coeff[1];
292 return makePtr<LinearCombinationObjective<Real>>(coeff,obj_vec);
293 }
294
295 int numGlobalSamples(void) const {
296 return numSamples_;
297 }
298};
299
300}
301
302#endif
Provides the std::vector implementation of the ROL::Vector interface.
ROL::Ptr< const std::vector< Real > > getAtom(const int i) const
void check(std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Provides a simplified interface for solving a wide range of optimization problems.
int solve(const ROL::Ptr< StatusTest< Real > > &status=ROL::nullPtr, const bool combineStatus=true)
Solve optimization problem with no iteration output.
Provides the std::vector implementation of the ROL::Vector interface.
const Real getProbability(const int i) const
SROMGenerator(ROL::ParameterList &parlist, const Ptr< BatchManager< Real > > &bman, const std::vector< Ptr< Distribution< Real > > > &dist, std::ostream &outStream=std::cout)
void initialize_vectors(Ptr< ProbabilityVector< Real > > &prob, Ptr< ProbabilityVector< Real > > &prob_lo, Ptr< ProbabilityVector< Real > > &prob_hi, Ptr< ProbabilityVector< Real > > &prob_eq, Ptr< AtomVector< Real > > &atom, Ptr< AtomVector< Real > > &atom_lo, Ptr< AtomVector< Real > > &atom_hi, Ptr< AtomVector< Real > > &atom_eq, Ptr< Vector< Real > > &vec, Ptr< Vector< Real > > &vec_lo, Ptr< Vector< Real > > &vec_hi, Ptr< Vector< Real > > &vec_eq, const Ptr< BatchManager< Real > > &bman) const
Ptr< Objective< Real > > initialize_objective(const std::vector< Ptr< Distribution< Real > > > &dist, const Ptr< BatchManager< Real > > &bman, const bool optProb, const bool optAtom, ROL::ParameterList &list) const
void get_scaling_vectors(std::vector< Real > &typw, std::vector< Real > &typx) const
std::vector< Ptr< Distribution< Real > > > dist_
void pruneSamples(const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom)
int numGlobalSamples(void) const
ROL::ParameterList parlist_
void setPoints(std::vector< std::vector< Real > > &p)
void setWeights(std::vector< Real > &w)
void sumAll(Real *input, Real *output, int dim) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84