ROL
ROL_MomentObjective.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_MOMENTOBJECTIVE_H
45 #define ROL_MOMENTOBJECTIVE_H
46 
47 #include "ROL_Objective.hpp"
48 #include "ROL_BatchManager.hpp"
49 #include "ROL_Distribution.hpp"
50 #include "ROL_SROMVector.hpp"
51 #include "ROL_Types.hpp"
52 #include <iostream>
53 
54 namespace ROL {
55 
56 template <class Real>
57 class MomentObjective : public Objective<Real> {
58 private:
59  std::vector<std::vector<std::pair<int, Real> > > moments_;
60  Teuchos::RCP<BatchManager<Real> > bman_;
63  const bool optProb_;
64  const bool optAtom_;
65 
66  Real momentValue(const int dim, const Real power, const Real moment,
67  const ProbabilityVector<Real> &prob,
68  const AtomVector<Real> &atom) const {
69  const int numSamples = prob.getNumMyAtoms();
70  Real val = 0., xpt = 0., xwt = 0., sum = 0.;
71  for (int k = 0; k < numSamples; k++) {
72  xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
73  val += xwt * ((power==1) ? xpt : std::pow(xpt,power));
74  }
75  bman_->sumAll(&val,&sum,1);
76  return 0.5*std::pow((sum-moment)/moment,2);
77  }
78 
79  void momentGradient(std::vector<Real> &gradx, std::vector<Real> &gradp, Real &scale,
80  const int dim, const Real power, const Real moment,
81  const ProbabilityVector<Real> &prob,
82  const AtomVector<Real> &atom) const {
83  const int numSamples = prob.getNumMyAtoms();
84  gradx.resize(numSamples,0.); gradp.resize(numSamples,0.);
85  scale = 0.;
86  Real xpt = 0., xwt = 0., xpow = 0., psum = 0.;
87  for (int k = 0; k < numSamples; k++) {
88  xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
89  xpow = ((power==1) ? 1. : ((power==2) ? xpt : std::pow(xpt,power-1)));
90  psum += xwt * xpow * xpt;
91  gradx[k] = xwt * xpow * power;
92  gradp[k] = xpow * xpt;
93  }
94  bman_->sumAll(&psum,&scale,1);
95  scale -= moment;
96  scale /= std::pow(moment,2);
97  }
98 
99  void momentHessVec(std::vector<Real> &hvx1, std::vector<Real> &hvx2, std::vector<Real> &hvx3,
100  std::vector<Real> &hvp1, std::vector<Real> &hvp2,
101  Real &scale1, Real &scale2, Real &scale3,
102  const int dim, const Real power, const Real moment,
103  const ProbabilityVector<Real> &prob,
104  const AtomVector<Real> &atom,
105  const ProbabilityVector<Real> &vprob,
106  const AtomVector<Real> &vatom) const {
107  const int numSamples = prob.getNumMyAtoms();
108  hvx1.resize(numSamples,0.); hvx2.resize(numSamples,0.); hvx3.resize(numSamples,0.);
109  hvp1.resize(numSamples,0.); hvp2.resize(numSamples,0.);
110  scale1 = 0.; scale2 = 0.; scale3 = 0.;
111  std::vector<Real> psum(3,0.0), scale(3,0.0);
112  Real xpt = 0., xwt = 0., vpt = 0., vwt = 0.;
113  Real xpow0 = 0., xpow1 = 0., xpow2 = 0.;
114  const Real moment2 = std::pow(moment,2);
115  for (int k = 0; k < numSamples; k++) {
116  xpt = (*atom.getAtom(k))[dim]; xwt = prob.getProbability(k);
117  vpt = (*vatom.getAtom(k))[dim]; vwt = vprob.getProbability(k);
118  xpow2 = ((power==1) ? 0. : ((power==2) ? 1. : ((power==3) ? xpt :
119  std::pow(xpt,power-2))));
120  xpow1 = ((power==1) ? 1. : xpow2 * xpt);
121  xpow0 = xpow1 * xpt;
122  psum[0] += xwt * xpow1 * vpt;
123  psum[1] += xwt * xpow0;
124  psum[2] += vwt * xpow0;
125  hvx1[k] = power * xwt * xpow1;
126  hvx2[k] = power * (power-1.) * xwt * xpow2 * vpt;
127  hvx3[k] = power * vwt * xpow1;
128  hvp1[k] = xpow0;
129  hvp2[k] = power * xpow1 * vpt;
130  }
131  bman_->sumAll(&psum[0],&scale[0],3);
132  scale1 = scale[0] * power/moment2;
133  scale2 = (scale[1] - moment)/moment2 ;
134  scale3 = scale[2]/moment2;
135  }
136 
137 public:
138  MomentObjective(const std::vector<std::vector<std::pair<int, Real> > > &moments,
139  const Teuchos::RCP<BatchManager<Real> > &bman,
140  const bool optProb = true, const bool optAtom = true)
141  : Objective<Real>(), moments_(moments), bman_(bman),
142  optProb_(optProb), optAtom_(optAtom) {
143  dimension_ = moments_.size();
144  numMoments_ = moments_[0].size();
145  }
146 
147  MomentObjective(const std::vector<Teuchos::RCP<Distribution<Real> > > &dist,
148  const std::vector<int> &order,
149  const Teuchos::RCP<BatchManager<Real> > &bman,
150  const bool optProb = true, const bool optAtom = true)
151  : Objective<Real>(), bman_(bman), optProb_(optProb), optAtom_(optAtom) {
152  numMoments_ = order.size();
153  dimension_ = dist.size();
154  std::vector<std::pair<int,Real> > data(numMoments_);
155  moments_.clear(); moments_.resize(dimension_);
156  for (int d = 0; d < dimension_; d++) {
157  for (int i = 0; i < numMoments_; i++) {
158  data[i] = std::make_pair(order[i],dist[d]->moment(order[i]));
159  }
160  moments_[d].assign(data.begin(),data.end());
161  }
162  }
163 
164  Real value( const Vector<Real> &x, Real &tol ) {
165  const SROMVector<Real> &ex = Teuchos::dyn_cast<const SROMVector<Real> >(x);
166  const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
167  const AtomVector<Real> &atom = *(ex.getAtomVector());
168  Real val = 0.;
169  std::vector<std::pair<int, Real> > data;
170  for (int d = 0; d < dimension_; d++) {
171  data = moments_[d];
172  for (int m = 0; m < numMoments_; m++) {
173  val += momentValue(d,(Real)data[m].first,data[m].second,prob,atom);
174  }
175  }
176  return val;
177  }
178 
179  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
180  g.zero();
181  const SROMVector<Real> &ex = Teuchos::dyn_cast<const SROMVector<Real> >(x);
182  const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
183  const AtomVector<Real> &atom = *(ex.getAtomVector());
184  int numSamples = prob.getNumMyAtoms();
185  std::vector<Real> gradx(numSamples,0.), gradp(numSamples,0.);
186  Real scale = 0.;
187  std::vector<std::pair<int, Real> > data;
188  std::vector<Real> val_wt(numSamples,0.), tmp(dimension_,0.);
189  std::vector<std::vector<Real> > val_pt(numSamples,tmp);
190  for (int d = 0; d < dimension_; d++) {
191  data = moments_[d];
192  for (int m = 0; m < numMoments_; m++) {
193  momentGradient(gradx,gradp,scale,d,(Real)data[m].first,data[m].second,prob,atom);
194  for (int k = 0; k < numSamples; k++) {
195  (val_pt[k])[d] += scale*gradx[k];
196  val_wt[k] += scale*gradp[k];
197  }
198  }
199  }
200  SROMVector<Real> &eg = Teuchos::dyn_cast<SROMVector<Real> >(g);
202  AtomVector<Real> &gatom = *(eg.getAtomVector());
203  for (int k = 0; k < numSamples; k++) {
204  if ( optProb_ ) {
205  gprob.setProbability(k,val_wt[k]);
206  }
207  if ( optAtom_ ) {
208  gatom.setAtom(k,val_pt[k]);
209  }
210  }
211  }
212 
213  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
214  hv.zero();
215  const SROMVector<Real> &ev = Teuchos::dyn_cast<const SROMVector<Real> >(v);
216  const ProbabilityVector<Real> &vprob = *(ev.getProbabilityVector());
217  const AtomVector<Real> &vatom = *(ev.getAtomVector());
218  const SROMVector<Real> &ex = Teuchos::dyn_cast<const SROMVector<Real> >(x);
219  const ProbabilityVector<Real> &prob = *(ex.getProbabilityVector());
220  const AtomVector<Real> &atom = *(ex.getAtomVector());
221  const int numSamples = prob.getNumMyAtoms();
222  std::vector<Real> hvx1(numSamples,0.), hvx2(numSamples,0.), hvx3(numSamples,0.);
223  std::vector<Real> hvp1(numSamples,0.), hvp2(numSamples,0.);
224  Real scale1 = 0., scale2 = 0., scale3 = 0.;
225  std::vector<std::pair<int, Real> > data;
226  std::vector<Real> val_wt(numSamples,0.), tmp(dimension_,0.);
227  std::vector<std::vector<Real> > val_pt(numSamples,tmp);
228  for (int d = 0; d < dimension_; d++) {
229  data = moments_[d];
230  for (int m = 0; m < numMoments_; m++) {
231  momentHessVec(hvx1,hvx2,hvx3,hvp1,hvp2,scale1,scale2,scale3,
232  d,(Real)data[m].first,data[m].second,prob,atom,vprob,vatom);
233  for (int k = 0; k < numSamples; k++) {
234  (val_pt[k])[d] += (scale1+scale3)*hvx1[k] + scale2*(hvx2[k]+hvx3[k]);
235  val_wt[k] += (scale1+scale3)*hvp1[k] + scale2*hvp2[k];
236  }
237  }
238  }
239  SROMVector<Real> &ehv = Teuchos::dyn_cast<SROMVector<Real> >(hv);
241  AtomVector<Real> &hatom = *(ehv.getAtomVector());
242  for (int k = 0; k < numSamples; k++) {
243  if ( optProb_ ) {
244  hprob.setProbability(k,val_wt[k]);
245  }
246  if ( optAtom_ ) {
247  hatom.setAtom(k,val_pt[k]);
248  }
249  }
250  }
251 }; // class SROMObjective
252 
253 } // namespace ROL
254 
255 #endif
Provides the interface to evaluate objective functions.
void momentGradient(std::vector< Real > &gradx, std::vector< Real > &gradp, Real &scale, const int dim, const Real power, const Real moment, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom) const
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Provides the std::vector implementation of the ROL::Vector interface.
void momentHessVec(std::vector< Real > &hvx1, std::vector< Real > &hvx2, std::vector< Real > &hvx3, std::vector< Real > &hvp1, std::vector< Real > &hvp2, Real &scale1, Real &scale2, Real &scale3, const int dim, const Real power, const Real moment, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom, const ProbabilityVector< Real > &vprob, const AtomVector< Real > &vatom) const
const Real getProbability(const int i) const
void setAtom(const int i, const std::vector< Real > &pt)
Contains definitions of custom data types in ROL.
const Teuchos::RCP< const ProbabilityVector< Real > > getProbabilityVector(void) const
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
std::vector< std::vector< std::pair< int, Real > > > moments_
Teuchos::RCP< const std::vector< Real > > getAtom(const int i) const
Teuchos::RCP< BatchManager< Real > > bman_
const Teuchos::RCP< const AtomVector< Real > > getAtomVector(void) const
Provides the std::vector implementation of the ROL::Vector interface.
Provides the std::vector implementation of the ROL::Vector interface.
Real momentValue(const int dim, const Real power, const Real moment, const ProbabilityVector< Real > &prob, const AtomVector< Real > &atom) const
int getNumMyAtoms(void) const
MomentObjective(const std::vector< Teuchos::RCP< Distribution< Real > > > &dist, const std::vector< int > &order, const Teuchos::RCP< BatchManager< Real > > &bman, const bool optProb=true, const bool optAtom=true)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
MomentObjective(const std::vector< std::vector< std::pair< int, Real > > > &moments, const Teuchos::RCP< BatchManager< Real > > &bman, const bool optProb=true, const bool optAtom=true)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.