ROL
ROL_HMCR.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_HMCR_HPP
11 #define ROL_HMCR_HPP
12 
14 #include "ROL_PlusFunction.hpp"
15 #include "ROL_RiskVector.hpp"
16 
38 namespace ROL {
39 
40 template<class Real>
41 class HMCR : public RandVarFunctional<Real> {
42 private:
43  // User inputs
44  ROL::Ptr<PlusFunction<Real> > plusFunction_;
45  Real prob_;
46  Real lambda_;
47  unsigned order_;
48 
49  // 1/(1-prob)
50  Real coeff_;
51 
52  // Temporary vector storage
53  ROL::Ptr<Vector<Real> > mDualVector_;
54 
55  // Temporary scalar storage
56  Real pnorm_;
57  Real coeff0_;
58  Real coeff1_;
59  Real coeff2_;
60 
61  // Flag to initialized vector storage
63 
69 
72 
77 
78  void checkInputs(void) const {
79  const Real zero(0), one(1);
80  ROL_TEST_FOR_EXCEPTION((prob_ <= zero) || (prob_ >= one), std::invalid_argument,
81  ">>> ERROR (ROL::HMCR): Confidence level must be between 0 and 1!");
82  ROL_TEST_FOR_EXCEPTION((lambda_ < zero) || (lambda_ > one), std::invalid_argument,
83  ">>> ERROR (ROL::HMCR): Convex combination parameter must be positive!");
84  ROL_TEST_FOR_EXCEPTION((order_ < 2), std::invalid_argument,
85  ">>> ERROR (ROL::HMCR): Norm order is less than 2!");
86  ROL_TEST_FOR_EXCEPTION(plusFunction_ == ROL::nullPtr, std::invalid_argument,
87  ">>> ERROR (ROL::HMCR): PlusFunction pointer is null!");
88  }
89 
90 public:
100  HMCR( const Real prob, const Real lambda, const unsigned order,
101  const ROL::Ptr<PlusFunction<Real> > &pf )
102  : RandVarFunctional<Real>(),
103  plusFunction_(pf), prob_(prob), lambda_(lambda), order_(order),
104  pnorm_(0), coeff0_(0), coeff1_(0), coeff2_(0),
105  HMCR_firstReset_(true) {
106  checkInputs();
107  const Real one(1);
108  coeff_ = one/(one-prob_);
109  }
110 
122  HMCR( ROL::ParameterList &parlist )
123  : RandVarFunctional<Real>(),
124  pnorm_(0), coeff0_(0), coeff1_(0), coeff2_(0),
125  HMCR_firstReset_(true) {
126  ROL::ParameterList &list
127  = parlist.sublist("SOL").sublist("Risk Measure").sublist("HMCR");
128  // Check HMCR inputs
129  prob_ = list.get<Real>("Confidence Level");
130  lambda_ = list.get<Real>("Convex Combination Parameter");
131  order_ = (unsigned)list.get<int>("Order",2);
132  // Build (approximate) plus function
133  plusFunction_ = ROL::makePtr<PlusFunction<Real> >(list);
134  // Check inputs
135  checkInputs();
136  const Real one(1);
137  coeff_ = one/(one-prob_);
138  }
139 
140  void initialize(const Vector<Real> &x) {
142  // Initialize additional vector storage
143  if ( HMCR_firstReset_ ) {
144  mDualVector_ = x.dual().clone();
145  HMCR_firstReset_ = false;
146  }
147  // Zero temporary storage
148  const Real zero(0);
149  mDualVector_->zero();
150  pnorm_ = zero; coeff0_ = zero;
151  coeff1_ = zero; coeff2_ = zero;
152  }
153 
155  const Vector<Real> &x,
156  const std::vector<Real> &xstat,
157  Real &tol) {
158  const Real rorder = static_cast<Real>(order_);
159  // Expected value
160  Real val = computeValue(obj,x,tol);
161  val_ += weight_*val;
162  // Higher moment
163  Real pf = plusFunction_->evaluate(val-xstat[0],0);
164  pnorm_ += weight_*std::pow(pf,rorder);
165  }
166 
167  Real getValue(const Vector<Real> &x,
168  const std::vector<Real> &xstat,
169  SampleGenerator<Real> &sampler) {
170  const Real one(1);
171  const Real power = one/static_cast<Real>(order_);
172  std::vector<Real> val_in(2), val_out(2);
173  val_in[0] = val_;
174  val_in[1] = pnorm_;
175  sampler.sumAll(&val_in[0],&val_out[0],2);
176  return (one-lambda_)*val_out[0]
177  + lambda_*(xstat[0] + coeff_*std::pow(val_out[1],power));
178  }
179 
181  const Vector<Real> &x,
182  const std::vector<Real> &xstat,
183  Real &tol) {
184  const Real one(1);
185  const Real rorder0 = static_cast<Real>(order_);
186  const Real rorder1 = rorder0 - one;
187  // Expected value
188  computeGradient(*dualVector_,obj,x,tol);
189  g_->axpy(weight_,*dualVector_);
190  // Higher moment
191  Real val = computeValue(obj,x,tol);
192  Real pf0 = plusFunction_->evaluate(val-xstat[0],0);
193  Real pf1 = plusFunction_->evaluate(val-xstat[0],1);
194 
195  Real pf0p0 = std::pow(pf0,rorder0);
196  Real pf0p1 = std::pow(pf0,rorder1);
197 
198  pnorm_ += weight_*pf0p0;
199  coeff0_ += weight_*pf0p1*pf1;
200 
201  hv_->axpy(weight_*pf0p1*pf1,*dualVector_);
202  }
203 
205  std::vector<Real> &gstat,
206  const Vector<Real> &x,
207  const std::vector<Real> &xstat,
208  SampleGenerator<Real> &sampler) {
209  const Real zero(0), one(1);
210  std::vector<Real> val_in(2), val_out(2);
211  val_in[0] = pnorm_; val_in[1] = coeff0_;
212 
213  sampler.sumAll(&val_in[0],&val_out[0],2);
214  sampler.sumAll(*g_,g);
215  g.scale(one-lambda_);
216  Real var = lambda_;
217  // If the higher moment term is positive then compute gradient
218  if ( val_out[0] > zero ) {
219  const Real rorder0 = static_cast<Real>(order_);
220  const Real rorder1 = rorder0 - one;
221  Real denom = std::pow(val_out[0],rorder1/rorder0);
222  // Sum higher moment contribution
223  dualVector_->zero();
224  sampler.sumAll(*hv_,*dualVector_);
225  g.axpy(lambda_*coeff_/denom,*dualVector_);
226  // Compute statistic gradient
227  var -= lambda_*coeff_*((denom > zero) ? val_out[1]/denom : zero);
228  }
229  gstat[0] = var;
230  }
231 
233  const Vector<Real> &v,
234  const std::vector<Real> &vstat,
235  const Vector<Real> &x,
236  const std::vector<Real> &xstat,
237  Real &tol) {
238  const Real one(1);
239  const Real rorder0 = static_cast<Real>(order_);
240  const Real rorder1 = rorder0-one;
241  const Real rorder2 = rorder1-one;
242  // Expected value
243  computeHessVec(*dualVector_,obj,v,x,tol);
244  hv_->axpy(weight_,*dualVector_);
245  // Higher moment
246  Real val = computeValue(obj,x,tol);
247  Real pf0 = plusFunction_->evaluate(val-xstat[0],0);
248  Real pf1 = plusFunction_->evaluate(val-xstat[0],1);
249  Real pf2 = plusFunction_->evaluate(val-xstat[0],2);
250 
251  Real pf0p0 = std::pow(pf0,rorder0);
252  Real pf0p1 = std::pow(pf0,rorder1);
253  Real pf0p2 = std::pow(pf0,rorder2);
254 
255  Real scale1 = pf0p1*pf1;
256  g_->axpy(weight_*scale1,*dualVector_);
257 
258  Real gv = computeGradVec(*dualVector_,obj,v,x,tol);
259  Real scale0 = (rorder1*pf0p2*pf1*pf1 + pf0p1*pf2)*(gv-vstat[0]);
260 
261  pnorm_ += weight_*pf0p0;
262  coeff0_ += weight_*scale0;
263  coeff1_ += weight_*scale1;
264  coeff2_ += weight_*rorder1*scale1*(vstat[0]-gv);
265 
266  g_->axpy(weight_*scale0,*dualVector_);
267  mDualVector_->axpy(weight_*scale1,*dualVector_);
268  }
269 
271  std::vector<Real> &hvstat,
272  const Vector<Real> &v,
273  const std::vector<Real> &vstat,
274  const Vector<Real> &x,
275  const std::vector<Real> &xstat,
276  SampleGenerator<Real> &sampler) {
277  const Real zero(0), one(1);
278  std::vector<Real> val_in(4), val_out(4);
279  val_in[0] = pnorm_; val_in[1] = coeff0_;
280  val_in[2] = coeff1_; val_in[3] = coeff2_;
281 
282  sampler.sumAll(&val_in[0],&val_out[0],4);
283  sampler.sumAll(*hv_,hv);
284 
285  Real var = zero;
286  hv.scale(one-lambda_);
287 
288  if ( val_out[0] > zero ) {
289  const Real rorder0 = static_cast<Real>(order_);
290  const Real rorder1 = rorder0-one;
291  const Real rorder2 = rorder0 + rorder1;
292  const Real coeff = lambda_*coeff_;
293 
294  Real denom1 = std::pow(val_out[0],rorder1/rorder0);
295  Real denom2 = std::pow(val_out[0],rorder2/rorder0);
296 
297  dualVector_->zero();
298  sampler.sumAll(*g_,*dualVector_);
299  hv.axpy(coeff/denom1,*dualVector_);
300 
301  dualVector_->zero();
302  sampler.sumAll(*mDualVector_,*dualVector_);
303  hv.axpy(coeff*val_out[3]/denom2,*dualVector_);
304 
305  var = -coeff*(val_out[1]/denom1 + val_out[3]*val_out[2]/denom2);
306  }
307  hvstat[0] = var;
308  }
309 };
310 
311 }
312 
313 #endif
Provides the interface to evaluate objective functions.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
void computeHessVec(Vector< Real > &hv, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
virtual void scale(const Real alpha)=0
Compute where .
Ptr< Vector< Real > > g_
Real computeValue(Objective< Real > &obj, const Vector< Real > &x, Real &tol)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
ROL::Ptr< PlusFunction< Real > > plusFunction_
Definition: ROL_HMCR.hpp:44
Ptr< Vector< Real > > hv_
Real getValue(const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure value.
Definition: ROL_HMCR.hpp:167
bool HMCR_firstReset_
Definition: ROL_HMCR.hpp:62
Ptr< Vector< Real > > dualVector_
ROL::Ptr< Vector< Real > > mDualVector_
Definition: ROL_HMCR.hpp:53
Real coeff1_
Definition: ROL_HMCR.hpp:58
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Real prob_
Definition: ROL_HMCR.hpp:45
void sumAll(Real *input, Real *output, int dim) const
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Real coeff_
Definition: ROL_HMCR.hpp:50
void initialize(const Vector< Real > &x)
Initialize temporary variables.
Definition: ROL_HMCR.hpp:140
void checkInputs(void) const
Definition: ROL_HMCR.hpp:78
void updateGradient(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for gradient computation.
Definition: ROL_HMCR.hpp:180
Real coeff2_
Definition: ROL_HMCR.hpp:59
void getHessVec(Vector< Real > &hv, std::vector< Real > &hvstat, const Vector< Real > &v, const std::vector< Real > &vstat, const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure Hessian-times-a-vector.
Definition: ROL_HMCR.hpp:270
Real pnorm_
Definition: ROL_HMCR.hpp:56
void updateHessVec(Objective< Real > &obj, const Vector< Real > &v, const std::vector< Real > &vstat, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal risk measure storage for Hessian-time-a-vector computation.
Definition: ROL_HMCR.hpp:232
Real lambda_
Definition: ROL_HMCR.hpp:46
HMCR(const Real prob, const Real lambda, const unsigned order, const ROL::Ptr< PlusFunction< Real > > &pf)
Constructor.
Definition: ROL_HMCR.hpp:100
unsigned order_
Definition: ROL_HMCR.hpp:47
Real coeff0_
Definition: ROL_HMCR.hpp:57
void computeGradient(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &x, Real &tol)
HMCR(ROL::ParameterList &parlist)
Constructor.
Definition: ROL_HMCR.hpp:122
void updateValue(Objective< Real > &obj, const Vector< Real > &x, const std::vector< Real > &xstat, Real &tol)
Update internal storage for value computation.
Definition: ROL_HMCR.hpp:154
Real computeGradVec(Vector< Real > &g, Objective< Real > &obj, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Provides the interface to implement any functional that maps a random variable to a (extended) real n...
Provides an interface for a convex combination of the expected value and the higher moment coherent r...
Definition: ROL_HMCR.hpp:41
virtual void initialize(const Vector< Real > &x)
Initialize temporary variables.
void getGradient(Vector< Real > &g, std::vector< Real > &gstat, const Vector< Real > &x, const std::vector< Real > &xstat, SampleGenerator< Real > &sampler)
Return risk measure (sub)gradient.
Definition: ROL_HMCR.hpp:204