ROL
ROL_SerialConstraint.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 #pragma once
45 #ifndef ROL_SERIALCONSTRAINT_HPP
46 #define ROL_SERIALCONSTRAINT_HPP
47 
51 #include "ROL_VectorWorkspace.hpp"
52 
60 namespace ROL {
61 
62 namespace details {
63 
64 //using namespace std;
65 
66 template<typename Real>
68 
71 
72 private:
73 
74  Ptr<DynamicConstraint<Real>> con_; // Constraint over a single time step
75  Ptr<Vector<Real>> u_initial_; // Initial condition for state
76  Ptr<Vector<Real>> u_zero_; // Zero vector the size of u_initial_
77  Ptr<std::vector<TimeStamp<Real>>> timeStamp_; // Set of all time stamps
78  VectorWorkspace<Real> workspace_; // For efficient cloning
79  size_type Nt_; // Number of time steps
80  bool skipInitialCond_; // skip the initial condition application
81 
82  PV& partition ( Vector<Real>& x ) const { return static_cast<PV&>(x); }
83  const PV& partition ( const Vector<Real>& x ) const { return static_cast<const PV&>(x); }
84 
85 public:
86 
88  const Vector<Real>& u_initial,
89  const Ptr<vector<TimeStamp<Real>>>& timeStamp ) :
90  con_(con), u_initial_(u_initial.clone()),
91  u_zero_(u_initial.clone()),
92  timeStamp_(timeStamp), Nt_(timeStamp_->size()),
93  skipInitialCond_(false) {
94  u_initial_->set(u_initial);
95  u_zero_->zero();
96  }
97 
98  size_type numTimeSteps(void) const { return Nt_; }
99 
100  const Vector<Real>& getInitialCondition() const { return *u_initial_; }
101  void setInitialCondition( const Vector<Real>& u_initial ) { u_initial_->set(u_initial); }
102 
103  bool getSkipInitialCondition() const { return skipInitialCond_; }
104  void setSkipInitialCondition(bool skip) { skipInitialCond_ = skip; }
105 
106  Ptr<vector<TimeStamp<Real>>> getTimeStamp() const { return timeStamp_; }
107  void setTimeStamp(const Ptr<vector<TimeStamp<Real>>>& timeStamp )
108  { timeStamp_ = timeStamp; Nt_ = timeStamp_->size(); }
109 
110  virtual void solve( Vector<Real>& c, Vector<Real>& u,
111  const Vector<Real>& z, Real& tol ) override {
112  auto& cp = partition(c);
113  auto& up = partition(u);
114  auto& zp = partition(z);
115 
116  if(!skipInitialCond_)
117  con_->solve( *(cp.get(0)), *u_initial_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
118 
119  for( size_type k=1; k<Nt_; ++k )
120  con_->solve( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
121  } // solve
122 
123 
124  virtual void value( Vector<Real>& c, const Vector<Real>& u,
125  const Vector<Real>& z, Real& tol ) override {
126 
127  auto& cp = partition(c);
128  auto& up = partition(u);
129  auto& zp = partition(z);
130 
131  if(!skipInitialCond_)
132  con_->value( *(cp.get(0)), *u_initial_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
133 
134  for( size_type k=1; k<Nt_; ++k )
135  con_->value( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
136  } // value
137 
138 
139  virtual void applyJacobian_1( Vector<Real>& jv, const Vector<Real>& v,
140  const Vector<Real>& u, const Vector<Real>& z,
141  Real& tol ) override {
142 
143  auto& jvp = partition(jv); auto& vp = partition(v);
144  auto& up = partition(u); auto& zp = partition(z);
145 
146  auto tmp = workspace_.clone(jvp.get(0));
147  auto& x = *tmp;
148 
149  if(!skipInitialCond_)
150  con_->applyJacobian_un( *(jvp.get(0)), *(vp.get(0)), *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
151 
152  for( size_type k=1; k<Nt_; ++k ) {
153  con_->applyJacobian_uo( x, *(vp.get(k-1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
154  con_->applyJacobian_un( *(jvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
155  jvp.get(k)->plus(x);
156  } // end for
157  } // applyJacobian_1
158 
159 
160  virtual void applyInverseJacobian_1( Vector<Real>& ijv, const Vector<Real>& v, const Vector<Real>& u,
161  const Vector<Real>& z, Real& tol) override {
162 
163  auto& ijvp = partition(ijv); auto& vp = partition(v);
164  auto& up = partition(u); auto& zp = partition(z);
165 
166  auto tmp = workspace_.clone(ijvp.get(0));
167  auto& x = *tmp;
168 
169  if(!skipInitialCond_)
170  con_->applyInverseJacobian_un( *(ijvp.get(0)), *(vp.get(0)), *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
171 
172  for( size_type k=1; k<Nt_; ++k ) {
173  con_->applyJacobian_uo( x, *(ijvp.get(k-1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
174  x.scale(-1.0);
175  x.plus( *(vp.get(k)) );
176  con_->applyInverseJacobian_un( *(ijvp.get(k)), x, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
177  } // end for
178  } // applyInverseJacobian_1
179 
180 
181  virtual void applyAdjointJacobian_1( Vector<Real>& ajv, const Vector<Real>& v, const Vector<Real>& u,
182  const Vector<Real>& z, const Vector<Real>& dualv, Real& tol) override {
183 
184  auto& ajvp = partition(ajv); auto& vp = partition(v);
185  auto& up = partition(u); auto& zp = partition(z);
186 
187  auto tmp = workspace_.clone(ajvp.get(0));
188  auto& x = *tmp;
189 
190  if(!skipInitialCond_)
191  con_->applyAdjointJacobian_un( *(ajvp.get(0)), *(vp.get(0)), *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
192 
193  for( size_type k=1; k<Nt_; ++k ) {
194  con_->applyAdjointJacobian_un( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
195  con_->applyAdjointJacobian_uo( x, *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
196  ajvp.get(k-1)->plus(x);
197  } // end for
198  } // applyAdjointJacobian_1
199 
201  const Vector<Real>& z, Real& tol) override {
202 
203  auto& iajvp = partition(iajv); auto& vp = partition(v);
204  auto& up = partition(u); auto& zp = partition(z);
205 
206  auto tmp = workspace_.clone(iajvp.get(0));
207  auto& x = *tmp;
208 
209  size_type k = Nt_-1;
210 
211  con_->applyInverseAdjointJacobian_un( *(iajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
212 
213  for( size_type k=Nt_-2; k>0; --k ) {
214  con_->applyAdjointJacobian_uo( x, *(iajvp.get(k+1)), *(up.get(k)), *(up.get(k+1)), *(zp.get(k+1)), timeStamp_->at(k+1) );
215  x.scale(-1.0);
216  x.plus( *(vp.get(k) ) );
217  con_->applyInverseAdjointJacobian_un( *(iajvp.get(k)), x, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
218  } // end for
219 
220  con_->applyAdjointJacobian_uo( x, *(iajvp.get(1)), *(up.get(0)), *(up.get(1)), *(zp.get(1)), timeStamp_->at(1) );
221  x.scale(-1.0);
222  x.plus( *(vp.get(0) ) );
223  if(!skipInitialCond_) {
224  con_->applyInverseAdjointJacobian_un( *(iajvp.get(0)), x, *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
225  }
226  else {
227  // this weird condition places iajvp in the final vector slot
228  iajvp.get(0)->set(x);
229  }
230 
231  } // applyInverseAdjointJacobian_1
232 
233 
234  virtual void applyJacobian_2( Vector<Real>& jv, const Vector<Real>& v, const Vector<Real>& u,
235  const Vector<Real>& z, Real &tol ) override {
236 
237  auto& jvp = partition(jv); auto& vp = partition(v);
238  auto& up = partition(u); auto& zp = partition(z);
239 
240  if(!skipInitialCond_)
241  con_->applyJacobian_z( *(jvp.get(0)), *(vp.get(0)), *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
242 
243  for( size_type k=1; k<Nt_; ++k )
244  con_->applyJacobian_z( *(jvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
245  } // applyJacobian_2
246 
247 
248  virtual void applyAdjointJacobian_2( Vector<Real>& ajv, const Vector<Real>& v, const Vector<Real>& u,
249  const Vector<Real>& z, Real& tol ) override {
250 
251  auto& ajvp = partition(ajv); auto& vp = partition(v);
252  auto& up = partition(u); auto& zp = partition(z);
253 
254  if(!skipInitialCond_)
255  con_->applyAdjointJacobian_z( *(ajvp.get(0)), *(vp.get(0)), *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
256 
257  for( size_type k=1; k<Nt_; ++k )
258  con_->applyAdjointJacobian_z( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
259  } // applyAdjointJacobian_2
260 
261 
262 
263  virtual void applyAdjointHessian_11( Vector<Real> &ahwv, const Vector<Real> &w, const Vector<Real> &v,
264  const Vector<Real> &u, const Vector<Real> &z, Real &tol) override {
265 
266  auto& ahwvp = partition(ahwv); auto& wp = partition(w);
267  auto& vp = partition(v); auto& up = partition(u);
268  auto& zp = partition(z);
269 
270  auto tmp = workspace_.clone(ahwvp.get(0));
271  auto& x = *tmp;
272 
273  if( !skipInitialCond_ ) {
274  con_->applyAdjointHessian_un_un( *(ahwvp.get(0)), *(wp.get(0)), *(vp.get(0)), *u_zero_, *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
275  con_->applyAdjointHessian_un_uo( x, *(wp.get(1)), *(vp.get(1)), *(up.get(0)), *(up.get(1)), *(zp.get(1)), timeStamp_->at(1) );
276  ahwvp.get(0)->plus(x);
277  con_->applyAdjointHessian_uo_uo( x, *(wp.get(1)), *(vp.get(0)), *(up.get(0)), *(up.get(1)), *(zp.get(1)), timeStamp_->at(1) );
278  ahwvp.get(0)->plus(x);
279 
280  }
281 
282  for( size_type k=1; k<Nt_; ++k ) {
283  con_->applyAdjointHessian_un_un( *(ahwvp.get(k)), *(wp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
284  con_->applyAdjointHessian_uo_un( x, *(wp.get(k)), *(vp.get(k-1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
285  ahwvp.get(k)->plus(x);
286 
287  if( k < Nt_-1 ) {
288  con_->applyAdjointHessian_un_uo( x, *(wp.get(k+1)), *(vp.get(k+1)), *(up.get(k)), *(up.get(k+1)), *(zp.get(k+1)), timeStamp_->at(k+1) );
289  ahwvp.get(k)->plus(x);
290  con_->applyAdjointHessian_uo_uo( x, *(wp.get(k+1)), *(vp.get(k)), *(up.get(k)), *(up.get(k+1)), *(zp.get(k+1)), timeStamp_->at(k+1) );
291  ahwvp.get(k)->plus(x);
292  }
293  }
294  } // applyAdjointHessian_11
295 
296 
298  const Vector<Real> &w,
299  const Vector<Real> &v,
300  const Vector<Real> &u,
301  const Vector<Real> &z,
302  Real &tol) override {}
303 
304 
306  const Vector<Real> &w,
307  const Vector<Real> &v,
308  const Vector<Real> &u,
309  const Vector<Real> &z,
310  Real &tol) override {}
311 
312  virtual void applyAdjointHessian_22(Vector<Real> &ahwv, const Vector<Real> &w, const Vector<Real> &v,
313  const Vector<Real> &u, const Vector<Real> &z, Real &tol) override {
314 
315  auto& ahwvp = partition(ahwv); auto& vp = partition(v); auto& wp = partition(w);
316  auto& up = partition(u); auto& zp = partition(z);
317 
318  con_->applyAdjointHessian_z_z( *(ahwvp.get(0)), *(wp.get(0)), *(vp.get(0)), *u_initial_,
319  *(up.get(0)), *(zp.get(0)), timeStamp_->at(0) );
320 
321  for( size_type k=1; k<Nt_; ++k ) {
322  con_->applyAdjointHessian_z_z( *(ahwvp.get(k)), *(wp.get(k)), *(vp.get(k)), *(up.get(k-1)),
323  *(up.get(k)), *(zp.get(k)), timeStamp_->at(k) );
324  }
325  } // applyAdjointHessian_22
326 
327 
328 }; // details::SerialConstraint
329 
330 } // namespace details
331 
332 template<typename Real>
334 
335 } // namespace ROL
336 
337 
338 #endif // ROL_SERIALCONSTRAINT_HPP
virtual void value(Vector< Real > &c, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Evaluate the constraint operator at .
typename PV< Real >::size_type size_type
void setTimeStamp(const Ptr< vector< TimeStamp< Real >>> &timeStamp)
Defines the time-dependent constraint operator interface for simulation-based optimization.
virtual void applyAdjointJacobian_2(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to vector . This is the primary interface...
Defines the linear algebra of vector space on a generic partitioned vector.
PV & partition(Vector< Real > &x) const
virtual void applyAdjointHessian_12(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the optimization-space derivative of the adjoint of the constraint simulation-space Jacobian at...
Contains local time step information.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
virtual void applyJacobian_2(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
Ptr< DynamicConstraint< Real > > con_
virtual void applyAdjointHessian_11(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the simulation-space derivative of the adjoint of the constraint simulation-space Jacobian at ...
VectorWorkspace< Real > workspace_
const Vector< Real > & getInitialCondition() const
SerialConstraint(const Ptr< DynamicConstraint< Real >> &con, const Vector< Real > &u_initial, const Ptr< vector< TimeStamp< Real >>> &timeStamp)
Ptr< std::vector< TimeStamp< Real > > > timeStamp_
virtual void applyAdjointJacobian_1(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &dualv, Real &tol) override
Apply the adjoint of the partial constraint Jacobian at , , to the vector . This is the secondary int...
virtual void applyJacobian_1(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the partial constraint Jacobian at , , to the vector .
Ptr< vector< TimeStamp< Real > > > getTimeStamp() const
void setInitialCondition(const Vector< Real > &u_initial)
virtual void solve(Vector< Real > &c, Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Given , solve for .
const PV & partition(const Vector< Real > &x) const
Evaluates ROL::DynamicConstraint over a sequential set of time intervals.
virtual void applyAdjointHessian_22(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the optimization-space derivative of the adjoint of the constraint optimization-space Jacobian ...
virtual void applyAdjointHessian_21(Vector< Real > &ahwv, const Vector< Real > &w, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the simulation-space derivative of the adjoint of the constraint optimization-space Jacobian at...
typename std::vector< Real >::size_type size_type
void applyInverseAdjointJacobian_1(Vector< Real > &iajv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the inverse of the adjoint of the partial constraint Jacobian at , , to the vector ...
virtual void applyInverseJacobian_1(Vector< Real > &ijv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol) override
Apply the inverse partial constraint Jacobian at , , to the vector .