ROL
ROL_MeanDeviationFromTarget.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_MEANDEVIATIONFROMTARGET_HPP
45 #define ROL_MEANDEVIATIONFROMTARGET_HPP
46 
47 #include "ROL_RiskMeasure.hpp"
48 #include "ROL_PositiveFunction.hpp"
49 #include "ROL_PlusFunction.hpp"
50 #include "ROL_AbsoluteValue.hpp"
51 
52 #include "Teuchos_ParameterList.hpp"
53 #include "Teuchos_Array.hpp"
54 
55 namespace ROL {
56 
57 template<class Real>
58 class MeanDeviationFromTarget : public RiskMeasure<Real> {
59  typedef typename std::vector<Real>::size_type uint;
60 private:
61  Teuchos::RCP<PositiveFunction<Real> > positiveFunction_;
62 
63  Teuchos::RCP<Vector<Real> > dualVector1_;
64  Teuchos::RCP<Vector<Real> > dualVector2_;
65  Teuchos::RCP<Vector<Real> > dualVector3_;
66  Teuchos::RCP<Vector<Real> > dualVector4_;
67 
68  std::vector<Real> target_;
69  std::vector<Real> order_;
70  std::vector<Real> coeff_;
72 
73  std::vector<Real> pval_;
74  std::vector<Real> pgv_;
75 
76  std::vector<Teuchos::RCP<Vector<Real> > > pg0_;
77  std::vector<Teuchos::RCP<Vector<Real> > > pg_;
78  std::vector<Teuchos::RCP<Vector<Real> > > phv_;
79 
81 
82  void initialize(void) {
83  // Initialize additional storage
84  pg_.clear(); pg0_.clear(); phv_.clear(); pval_.clear(); pgv_.clear();
85  pg_.resize(NumMoments_);
86  pg0_.resize(NumMoments_);
87  phv_.resize(NumMoments_);
88  pval_.resize(NumMoments_,0.0);
89  pgv_.resize(NumMoments_,0.0);
90  }
91 
92 public:
93 
94  MeanDeviationFromTarget( Real target, Real order, Real coeff,
95  Teuchos::RCP<PositiveFunction<Real> > &pf )
96  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
97  // Initialize storage for problem data
98  target_.clear(); order_.clear(); coeff_.clear();
99  target_.push_back(target);
100  order_.push_back((order < 2.0) ? 2.0 : order);
101  coeff_.push_back((coeff < 0.0) ? 1.0 : coeff);
102  NumMoments_ = order_.size();
103  initialize();
104  }
105 
106  MeanDeviationFromTarget( std::vector<Real> &target, std::vector<Real> &order, std::vector<Real> &coeff,
107  Teuchos::RCP<PositiveFunction<Real> > &pf )
108  : RiskMeasure<Real>(), positiveFunction_(pf), firstReset_(true) {
109  // Initialize storage for problem data
110  NumMoments_ = order.size();
111  target_.clear(); order_.clear(); coeff_.clear();
112  if ( NumMoments_ != target.size() ) {
113  target.resize(NumMoments_,0.0);
114  }
115  if ( NumMoments_ != coeff.size() ) {
116  coeff.resize(NumMoments_,1.0);
117  }
118  for ( uint i = 0; i < NumMoments_; i++ ) {
119  target_.push_back(target[i]);
120  order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
121  coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
122  }
123  initialize();
124  }
125 
126  MeanDeviationFromTarget( Teuchos::ParameterList &parlist )
127  : RiskMeasure<Real>(), firstReset_(true) {
128  Teuchos::ParameterList &list
129  = parlist.sublist("SOL").sublist("Risk Measure").sublist("Mean Plus Deviation From Target");
130  // Get data from parameter list
131  Teuchos::Array<Real> target
132  = Teuchos::getArrayFromStringParameter<double>(list,"Targets");
133  Teuchos::Array<Real> order
134  = Teuchos::getArrayFromStringParameter<double>(list,"Orders");
135  Teuchos::Array<Real> coeff
136  = Teuchos::getArrayFromStringParameter<double>(list,"Coefficients");
137  // Check inputs
138  NumMoments_ = order.size();
139  target_.clear(); order_.clear(); coeff_.clear();
140  if ( NumMoments_ != static_cast<uint>(target.size()) ) {
141  target.resize(NumMoments_,0.0);
142  }
143  if ( NumMoments_ != static_cast<uint>(coeff.size()) ) {
144  coeff.resize(NumMoments_,1.0);
145  }
146  for ( uint i = 0; i < NumMoments_; i++ ) {
147  target_.push_back(target[i]);
148  order_.push_back((order[i] < 2.0) ? 2.0 : order[i]);
149  coeff_.push_back((coeff[i] < 0.0) ? 1.0 : coeff[i]);
150  }
151  initialize();
152  // Build (approximate) positive function
153  if ( list.get("Deviation Type","Upper") == "Upper" ) {
154  positiveFunction_ = Teuchos::rcp(new PlusFunction<Real>(list));
155  }
156  else {
157  positiveFunction_ = Teuchos::rcp(new AbsoluteValue<Real>(list));
158  }
159  }
160 
161  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x) {
163  if (firstReset_) {
164  for ( uint p = 0; p < NumMoments_; p++ ) {
165  pg0_[p] = (x0->dual()).clone();
166  pg_[p] = (x0->dual()).clone();
167  phv_[p] = (x0->dual()).clone();
168  }
169  dualVector1_ = (x0->dual()).clone();
170  dualVector2_ = (x0->dual()).clone();
171  dualVector3_ = (x0->dual()).clone();
172  dualVector4_ = (x0->dual()).clone();
173  firstReset_ = false;
174  }
175  for ( uint p = 0; p < NumMoments_; p++ ) {
176  pg0_[p]->zero(); pg_[p]->zero(); phv_[p]->zero();
177  pval_[p] = 0.0; pgv_[p] = 0.0;
178  }
179  dualVector1_->zero(); dualVector2_->zero();
180  dualVector3_->zero(); dualVector4_->zero();
181  }
182 
183  void reset(Teuchos::RCP<Vector<Real> > &x0, const Vector<Real> &x,
184  Teuchos::RCP<Vector<Real> > &v0, const Vector<Real> &v) {
185  reset(x0,x);
186  v0 = Teuchos::rcp_const_cast<Vector<Real> >(Teuchos::dyn_cast<const RiskVector<Real> >(
187  Teuchos::dyn_cast<const Vector<Real> >(v)).getVector());
188  }
189 
190  void update(const Real val, const Real weight) {
191  Real diff = 0.0, pf0 = 0.0;
192  RiskMeasure<Real>::val_ += weight * val;
193  for ( uint p = 0; p < NumMoments_; p++ ) {
194  diff = val-target_[p];
195  pf0 = positiveFunction_->evaluate(diff,0);
196  pval_[p] += weight * std::pow(pf0,order_[p]);
197  }
198  }
199 
200  void update(const Real val, const Vector<Real> &g, const Real weight) {
201  Real diff = 0.0, pf0 = 0.0, pf1 = 0.0, c = 0.0;
202  for ( uint p = 0; p < NumMoments_; p++ ) {
203  diff = val-target_[p];
204  pf0 = positiveFunction_->evaluate(diff,0);
205  pf1 = positiveFunction_->evaluate(diff,1);
206  c = std::pow(pf0,order_[p]-1.0) * pf1;
207  (pg_[p])->axpy(weight * c,g);
208  pval_[p] += weight * std::pow(pf0,order_[p]);
209  }
210  RiskMeasure<Real>::g_->axpy(weight,g);
211  }
212 
213  void update(const Real val, const Vector<Real> &g, const Real gv, const Vector<Real> &hv,
214  const Real weight) {
215  Real diff = 0.0, pf0 = 0.0, pf1 = 0.0, pf2 = 0.0, p0 = 0.0, p1 = 0.0, p2 = 0.0, c = 0.0;
216  for ( uint p = 0; p < NumMoments_; p++ ) {
217  diff = val - target_[p];
218  pf0 = positiveFunction_->evaluate(diff,0);
219  pf1 = positiveFunction_->evaluate(diff,1);
220  pf2 = positiveFunction_->evaluate(diff,2);
221  p0 = std::pow(pf0,order_[p]);
222  p1 = std::pow(pf0,order_[p]-1.0);
223  p2 = std::pow(pf0,order_[p]-2.0);
224  c = -(order_[p]-1.0)*p1*pf1;
225  pg0_[p]->axpy(weight*c,g);
226  c = gv*((order_[p]-1.0)*p2*pf1*pf1 + p1*pf2);
227  pg_[p]->axpy(weight*c,g);
228  c = p1*pf1;
229  phv_[p]->axpy(weight*c,hv);
230  pval_[p] += weight*p0;
231  pgv_[p] += weight*p1*pf1*gv;
232  }
233  RiskMeasure<Real>::hv_->axpy(weight,hv);
234  }
235 
237  Real val = RiskMeasure<Real>::val_;
238  Real dev = 0.0;
239  sampler.sumAll(&val,&dev,1);
240  std::vector<Real> pval_sum(NumMoments_);
241  sampler.sumAll(&(pval_)[0],&pval_sum[0],NumMoments_);
242  for ( uint p = 0; p < NumMoments_; p++ ) {
243  dev += coeff_[p] * std::pow(pval_sum[p],1.0/order_[p]);
244  }
245  return dev;
246  }
247 
249  sampler.sumAll(*(RiskMeasure<Real>::g_),*dualVector1_);
250  std::vector<Real> pval_sum(NumMoments_);
251  sampler.sumAll(&(pval_)[0],&pval_sum[0],NumMoments_);
252  Teuchos::RCP<Vector<Real> > pg;
253  for ( uint p = 0; p < NumMoments_; p++ ) {
254  if ( pval_sum[p] > 0.0 ) {
255  pg = (pg_[p])->clone();
256  sampler.sumAll(*(pg_[p]),*pg);
257  dualVector1_->axpy(coeff_[p]/std::pow(pval_sum[p],1.0-1.0/order_[p]),*pg);
258  }
259  }
260  // Set RiskVector
261  (Teuchos::dyn_cast<RiskVector<Real> >(g)).setVector(*dualVector1_);
262  }
263 
265  sampler.sumAll(*(RiskMeasure<Real>::hv_),*dualVector1_);
266  std::vector<Real> pval_sum(NumMoments_);
267  sampler.sumAll(&(pval_)[0],&pval_sum[0],NumMoments_);
268  std::vector<Real> pgv_sum(NumMoments_);
269  sampler.sumAll(&(pgv_)[0],&pgv_sum[0],NumMoments_);
270  Real c = 0.0;
271  for ( uint p = 0; p < NumMoments_; p++ ) {
272  if ( pval_sum[p] > 0.0 ) {
273  sampler.sumAll(*(pg_[p]),*dualVector2_);
274  sampler.sumAll(*(pg0_[p]),*dualVector3_);
275  sampler.sumAll(*(phv_[p]),*dualVector4_);
276  c = coeff_[p]*(pgv_sum[p]/std::pow(pval_sum[p],2.0-1.0/order_[p]));
277  dualVector1_->axpy(c,*dualVector3_);
278  c = coeff_[p]/std::pow(pval_sum[p],1.0-1.0/order_[p]);
279  dualVector1_->axpy(c,*dualVector2_);
280  dualVector1_->axpy(c,*dualVector4_);
281  }
282  }
283  // Set RiskVector
284  (Teuchos::dyn_cast<RiskVector<Real> >(hv)).setVector(*dualVector1_);
285  }
286 };
287 
288 }
289 
290 #endif
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
MeanDeviationFromTarget(std::vector< Real > &target, std::vector< Real > &order, std::vector< Real > &coeff, Teuchos::RCP< PositiveFunction< Real > > &pf)
void getGradient(Vector< Real > &g, SampleGenerator< Real > &sampler)
void update(const Real val, const Vector< Real > &g, const Real gv, const Vector< Real > &hv, const Real weight)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
void sumAll(Real *input, Real *output, int dim) const
Teuchos::RCP< const Vector< Real > > getVector() const
std::vector< Real >::size_type uint
Teuchos::RCP< Vector< Real > > dualVector1_
MeanDeviationFromTarget(Real target, Real order, Real coeff, Teuchos::RCP< PositiveFunction< Real > > &pf)
Real getValue(SampleGenerator< Real > &sampler)
std::vector< Teuchos::RCP< Vector< Real > > > phv_
void update(const Real val, const Real weight)
Teuchos::RCP< Vector< Real > > dualVector2_
Teuchos::RCP< PositiveFunction< Real > > positiveFunction_
Teuchos::RCP< Vector< Real > > dualVector4_
void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x, Teuchos::RCP< Vector< Real > > &v0, const Vector< Real > &v)
virtual void reset(Teuchos::RCP< Vector< Real > > &x0, const Vector< Real > &x)
MeanDeviationFromTarget(Teuchos::ParameterList &parlist)
Teuchos::RCP< Vector< Real > > dualVector3_
std::vector< Teuchos::RCP< Vector< Real > > > pg_
std::vector< Teuchos::RCP< Vector< Real > > > pg0_
void getHessVec(Vector< Real > &hv, SampleGenerator< Real > &sampler)
void update(const Real val, const Vector< Real > &g, const Real weight)