ROL
ROL_PH_bPOEObjective.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 PH_BPOEOBJECTIVE_H
11 #define PH_BPOEOBJECTIVE_H
12 
13 #include "ROL_Objective.hpp"
14 
21 namespace ROL {
22 
23 template <class Real>
24 class PH_bPOEObjective : public Objective<Real> {
25 private:
26  const Ptr<Objective<Real>> obj_;
27  Real threshold_;
28  Real order_;
29 
31  Real val_;
32 
35  Ptr<Vector<Real>> g_;
36 
37  void getValue(const Vector<Real> &x, Real &tol) {
38  if (!isValueComputed_) {
39  val_ = obj_->value(x,tol);
40  isValueComputed_ = true;
41  }
42  }
43 
44  void getGradient(const Vector<Real> &x, Real &tol) {
46  g_ = x.dual().clone();
48  }
49  if (!isGradientComputed_) {
50  obj_->gradient(*g_,x,tol);
51  isGradientComputed_ = true;
52  }
53  }
54 
55  Ptr<const Vector<Real>> getConstVector(const Vector<Real> &x) const {
56  const RiskVector<Real> &xrv = dynamic_cast<const RiskVector<Real>&>(x);
57  return xrv.getVector();
58  }
59 
60  Ptr<Vector<Real>> getVector(Vector<Real> &x) const {
61  RiskVector<Real> &xrv = dynamic_cast<RiskVector<Real>&>(x);
62  return xrv.getVector();
63  }
64 
65  Ptr<const std::vector<Real>> getConstStat(const Vector<Real> &x) const {
66  const RiskVector<Real> &xrv = dynamic_cast<const RiskVector<Real>&>(x);
67  Ptr<const std::vector<Real>> xstat = xrv.getStatistic();
68  if (xstat == nullPtr) {
69  xstat = makePtr<const std::vector<Real>>(0);
70  }
71  return xstat;
72  }
73 
74  Ptr<std::vector<Real>> getStat(Vector<Real> &x) const {
75  RiskVector<Real> &xrv = dynamic_cast<RiskVector<Real>&>(x);
76  Ptr<std::vector<Real>> xstat = xrv.getStatistic();
77  if (xstat == nullPtr) {
78  xstat = makePtr<std::vector<Real>>(0);
79  }
80  return xstat;
81  }
82 
83  // pth power of the positive part function
84  Real pplus(const Real x, const int deriv = 0) const {
85  const Real zero(0), one(1), two(2), three(3);
86  Real val(0);
87  if (x > zero) {
88  if (deriv==0) {
89  val = (order_==one ? x
90  : std::pow(x,order_));
91  }
92  else if (deriv==1) {
93  val = order_*(order_==one ? one
94  : (order_==two ? x
95  : std::pow(x,order_-one)));
96  }
97  else if (deriv==2) {
98  val = order_*(order_-one)*(order_==one ? zero
99  : (order_==two ? one
100  : (order_==three ? x
101  : std::pow(x,order_-two))));
102  }
103  }
104  return val;
105  }
106 
107  Real bPOEobjective(const Real t, const Real x, const int deriv = 0) const {
108  const Real one(1);
109  Real arg = t*(x-threshold_)+one;
110  return pplus(arg,deriv);
111  }
112 
113 public:
114 
116  ParameterList &parlist)
117  : obj_(obj),
118  isValueComputed_(false),
119  isGradientInitialized_(false),
120  isGradientComputed_(false) {
121  ParameterList &list = parlist.sublist("SOL").sublist("Probability").sublist("bPOE");
122  threshold_ = list.get<Real>("Threshold");
123  order_ = list.get<Real>("Moment Order");
124  }
125 
126  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
127  Ptr<const Vector<Real>> xvec = getConstVector(x);
128  obj_->update(*xvec,flag,iter);
129  isValueComputed_ = false;
130  isGradientComputed_ = false;
131  }
132 
133  Real value( const Vector<Real> &x, Real &tol ) {
134  Ptr<const Vector<Real>> xvec = getConstVector(x);
135  Ptr<const std::vector<Real>> xstat = getConstStat(x);
136  getValue(*xvec,tol);
137  Real xt = (*xstat)[0];
138  Real prob = bPOEobjective(xt,val_,0);
139  return prob;
140  }
141 
142  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
143  Ptr<Vector<Real>> gvec = getVector(g);
144  Ptr<std::vector<Real>> gstat = getStat(g);
145  Ptr<const Vector<Real>> xvec = getConstVector(x);
146  Ptr<const std::vector<Real>> xstat = getConstStat(x);
147  getValue(*xvec,tol);
148  Real xt = (*xstat)[0], diff = val_-threshold_;
149  Real prob = bPOEobjective(xt,val_,1);
150  getGradient(*xvec,tol);
151  gvec->set(*g_);
152  gvec->scale(prob*xt);
153  (*gstat)[0] = prob*diff;
154  }
155 
156  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
157  Ptr<Vector<Real>> hvec = getVector(hv);
158  Ptr<std::vector<Real>> hstat = getStat(hv);
159  Ptr<const Vector<Real>> vvec = getConstVector(v);
160  Ptr<const std::vector<Real>> vstat = getConstStat(v);
161  Ptr<const Vector<Real>> xvec = getConstVector(x);
162  Ptr<const std::vector<Real>> xstat = getConstStat(x);
163  getValue(*xvec,tol);
164  Real xt = (*xstat)[0], vt = (*vstat)[0], diff = val_-threshold_;
165  Real prob1 = bPOEobjective(xt,val_,1);
166  Real prob2 = bPOEobjective(xt,val_,2);
167  getGradient(*xvec,tol);
168  //Real gv = vvec->dot(g_->dual());
169  Real gv = vvec->apply(*g_);
170  obj_->hessVec(*hvec,*vvec,*xvec,tol);
171  hvec->scale(prob1*xt);
172  hvec->axpy(prob2*xt*(vt*diff+xt*gv)+vt*prob1,*g_);
173  (*hstat)[0] = prob2*std::pow(diff,2)*vt+(prob2*diff*xt+prob1)*gv;
174  }
175 
176  void setParameter(const std::vector<Real> &param) {
177  obj_->setParameter(param);
179  }
180 
181 };
182 
183 }
184 #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
Ptr< Vector< Real > > g_
void getValue(const Vector< Real > &x, Real &tol)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Real pplus(const Real x, const int deriv=0) const
void setParameter(const std::vector< Real > &param)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Provides the interface for the progressive hedging probability objective.
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Ptr< Vector< Real > > getVector(Vector< Real > &x) const
Ptr< std::vector< Real > > getStatistic(const int comp=0, const int index=0)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Ptr< std::vector< Real > > getStat(Vector< Real > &x) const
PH_bPOEObjective(const Ptr< Objective< Real >> &obj, ParameterList &parlist)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
void getGradient(const Vector< Real > &x, Real &tol)
Real bPOEobjective(const Real t, const Real x, const int deriv=0) const
virtual void setParameter(const std::vector< Real > &param)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Ptr< const Vector< Real > > getVector(void) const
Ptr< const Vector< Real > > getConstVector(const Vector< Real > &x) const
const Ptr< Objective< Real > > obj_
Ptr< const std::vector< Real > > getConstStat(const Vector< Real > &x) const