ROL_Reduced_EqualityConstraint_SimOpt.hpp

Go to the documentation of this file.

// @HEADER
// ************************************************************************
//
//               Rapid Optimization Library (ROL) Package
//                 Copyright (2014) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact lead developers:
//              Drew Kouri   (dpkouri@sandia.gov) and
//              Denis Ridzal (dridzal@sandia.gov)
//
// ************************************************************************
// @HEADER


#ifndef ROL_REDUCED_EQUALITYCONSTRAINT_SIMOPT_H
#define ROL_REDUCED_EQUALITYCONSTRAINT_SIMOPT_H

#include "ROL_EqualityConstraint_SimOpt.hpp"
#include "ROL_SimController.hpp"

namespace ROL {

template <class Real>
class Reduced_EqualityConstraint_SimOpt : public EqualityConstraint<Real> {
private:
  const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > conVal_;          
  const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > conRed_; 
  const Teuchos::RCP<SimController<Real> > stateStore_;
  Teuchos::RCP<SimController<Real> > adjointStore_;

  // Primal vectors
  Teuchos::RCP<Vector<Real> > state_;                              
  Teuchos::RCP<Vector<Real> > adjoint_;                            
  Teuchos::RCP<Vector<Real> > residual_;                            
  Teuchos::RCP<Vector<Real> > state_sens_;                              
  Teuchos::RCP<Vector<Real> > adjoint_sens_;                            

  // Dual vectors
  Teuchos::RCP<Vector<Real> > dualstate_;
  Teuchos::RCP<Vector<Real> > dualstate1_;
  Teuchos::RCP<Vector<Real> > dualadjoint_;
  Teuchos::RCP<Vector<Real> > dualcontrol_;
  Teuchos::RCP<Vector<Real> > dualresidual_;                            

  const bool storage_;             
  const bool useFDhessVec_;

  bool updateFlag_;
  int  updateIter_;

  void solve_state_equation(const Vector<Real> &z, Real &tol) { 
    // Check if state has been computed.
    bool isComputed = false;
    if (storage_) {
      isComputed = stateStore_->get(*state_,EqualityConstraint<Real>::getParameter());
    }
    // Solve state equation if not done already.
    if (!isComputed || !storage_) {
      // Update equality constraint with new Opt variable.
      conRed_->update_2(z,updateFlag_,updateIter_);
      // Solve state equation.
      conRed_->solve(*dualadjoint_,*state_,z,tol);
      // Update equality constraint with new Sim variable.
      conRed_->update_1(*state_,updateFlag_,updateIter_);
      // Update full objective function.
      conVal_->update(*state_,z,updateFlag_,updateIter_);
      // Store state.
      if (storage_) {
        stateStore_->set(*state_,EqualityConstraint<Real>::getParameter());
      }
    }
  }

  void solve_adjoint_equation(const Vector<Real> &w, const Vector<Real> &z, Real &tol) { 
    // Check if adjoint has been computed.
    bool isComputed = false;
    if (storage_) {
      isComputed = adjointStore_->get(*adjoint_,EqualityConstraint<Real>::getParameter());
    }
    // Solve adjoint equation if not done already.
    if (!isComputed || !storage_) {
      // Evaluate the full gradient wrt u
      conVal_->applyAdjointJacobian_1(*dualstate_,w,*state_,z,tol);
      // Solve adjoint equation
      conRed_->applyInverseAdjointJacobian_1(*adjoint_,*dualstate_,*state_,z,tol);
      adjoint_->scale(static_cast<Real>(-1));
      // Store adjoint
      if (storage_) {
        adjointStore_->set(*adjoint_,EqualityConstraint<Real>::getParameter());
      }
    }
  }

  void solve_state_sensitivity(const Vector<Real> &v, const Vector<Real> &z, Real &tol) {
    // Solve state sensitivity equation
    conRed_->applyJacobian_2(*dualadjoint_,v,*state_,z,tol);
    dualadjoint_->scale(static_cast<Real>(-1));
    conRed_->applyInverseJacobian_1(*state_sens_,*dualadjoint_,*state_,z,tol);
  }

  void solve_adjoint_sensitivity(const Vector<Real> &w, const Vector<Real> &v, const Vector<Real> &z, Real &tol) {
    // Evaluate full hessVec in the direction (s,v)
    conVal_->applyAdjointHessian_11(*dualstate_,w,*state_sens_,*state_,z,tol);
    conVal_->applyAdjointHessian_12(*dualstate1_,w,v,*state_,z,tol);
    dualstate_->plus(*dualstate1_);
    // Apply adjoint Hessian of constraint
    conRed_->applyAdjointHessian_11(*dualstate1_,*adjoint_,*state_sens_,*state_,z,tol);
    dualstate_->plus(*dualstate1_);
    conRed_->applyAdjointHessian_21(*dualstate1_,*adjoint_,v,*state_,z,tol);
    dualstate_->plus(*dualstate1_);
    // Solve adjoint sensitivity equation
    dualstate_->scale(static_cast<Real>(-1));
    conRed_->applyInverseAdjointJacobian_1(*adjoint_sens_,*dualstate_,*state_,z,tol);
  }

public:
  Reduced_EqualityConstraint_SimOpt(
      const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > &conVal, 
      const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > &conRed, 
      const Teuchos::RCP<SimController<Real> > &stateStore,
      const Teuchos::RCP<Vector<Real> > &state, 
      const Teuchos::RCP<Vector<Real> > &control, 
      const Teuchos::RCP<Vector<Real> > &adjoint,
      const Teuchos::RCP<Vector<Real> > &residual,
      const bool storage = true,
      const bool useFDhessVec = false) 
    : conVal_(conVal), conRed_(conRed), stateStore_(stateStore),
      storage_(storage), useFDhessVec_(useFDhessVec),
      updateFlag_(true), updateIter_(0) {
    adjointStore_ = Teuchos::rcp(new SimController<Real>());
    state_        = state->clone();
    adjoint_      = adjoint->clone();
    residual_     = residual->clone();
    state_sens_   = state->clone();
    adjoint_sens_ = adjoint->clone();
    dualstate_    = state->dual().clone();
    dualstate1_   = state->dual().clone();
    dualadjoint_  = adjoint->dual().clone();
    dualcontrol_  = control->dual().clone();
    dualresidual_ = residual->dual().clone();
  }

  Reduced_EqualityConstraint_SimOpt(
      const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > &conVal, 
      const Teuchos::RCP<EqualityConstraint_SimOpt<Real> > &conRed,
      const Teuchos::RCP<SimController<Real> > &stateStore, 
      const Teuchos::RCP<Vector<Real> > &state,
      const Teuchos::RCP<Vector<Real> > &control, 
      const Teuchos::RCP<Vector<Real> > &adjoint,
      const Teuchos::RCP<Vector<Real> > &residual,
      const Teuchos::RCP<Vector<Real> > &dualstate,
      const Teuchos::RCP<Vector<Real> > &dualcontrol, 
      const Teuchos::RCP<Vector<Real> > &dualadjoint,
      const Teuchos::RCP<Vector<Real> > &dualresidual,
      const bool storage = true,
      const bool useFDhessVec = false)
    : conVal_(conVal), conRed_(conRed), stateStore_(stateStore),
      storage_(storage), useFDhessVec_(useFDhessVec),
      updateFlag_(true), updateIter_(0) {
    adjointStore_ = Teuchos::rcp(new SimController<Real>());
    state_        = state->clone();
    adjoint_      = adjoint->clone();
    residual_     = residual->clone();
    state_sens_   = state->clone();
    adjoint_sens_ = adjoint->clone();
    dualstate_    = dualstate->clone();
    dualstate1_   = dualstate->clone();
    dualadjoint_  = dualadjoint->clone();
    dualcontrol_  = dualcontrol->clone();
    dualresidual_ = dualresidual->clone();
  }

  void update( const Vector<Real> &z, bool flag = true, int iter = -1 ) {
    updateFlag_ = flag;
    updateIter_ = iter;
    stateStore_->equalityConstraintUpdate(true);
    adjointStore_->equalityConstraintUpdate(flag);
  }

  void value( Vector<Real> &c, const Vector<Real> &z, Real &tol ) {
    // Solve state equation
    solve_state_equation(z,tol);
    // Get constraint value
    conVal_->value(c,*state_,z,tol);
  }

  void applyJacobian( Vector<Real> &jv, const Vector<Real> &v,
                      const Vector<Real> &z, Real &tol ) {
    // Solve state equation.
    solve_state_equation(z,tol);
    // Solve state sensitivity equation.
    solve_state_sensitivity(v,z,tol);
    // Apply Sim Jacobian to state sensitivity.
    conVal_->applyJacobian_1(*residual_,*state_sens_,*state_,z,tol);
    // Apply Opt Jacobian to vector.
    conVal_->applyJacobian_2(jv,v,*state_,z,tol);
    jv.plus(*residual_);
  }

  void applyAdjointJacobian( Vector<Real> &ajw, const Vector<Real> &w,
                             const Vector<Real> &z, Real &tol ) {
    // Solve state equation
    solve_state_equation(z,tol);
    // Solve adjoint equation
    solve_adjoint_equation(w,z,tol);
    // Evaluate the full gradient wrt z
    conVal_->applyAdjointJacobian_2(*dualcontrol_,w,*state_,z,tol);
    // Build gradient
    conRed_->applyAdjointJacobian_2(ajw,*adjoint_,*state_,z,tol);
    ajw.plus(*dualcontrol_);
  }

  void applyAdjointHessian( Vector<Real> &ahwv, const Vector<Real> &w,
                            const Vector<Real> &v, const Vector<Real> &z,
                            Real &tol ) {
    if ( useFDhessVec_ ) {
      EqualityConstraint<Real>::applyAdjointHessian(ahwv,w,v,z,tol);
    }
    else {
      // Solve state equation
      solve_state_equation(z,tol);
      // Solve adjoint equation
      solve_adjoint_equation(w,z,tol);
      // Solve state sensitivity equation
      solve_state_sensitivity(v,z,tol);
      // Solve adjoint sensitivity equation
      solve_adjoint_sensitivity(w,v,z,tol);
      // Build hessVec
      conRed_->applyAdjointJacobian_2(ahwv,*adjoint_sens_,*state_,z,tol);
      conVal_->applyAdjointHessian_21(*dualcontrol_,w,*state_sens_,*state_,z,tol);
      ahwv.plus(*dualcontrol_);
      conVal_->applyAdjointHessian_22(*dualcontrol_,w,v,*state_,z,tol);
      ahwv.plus(*dualcontrol_);
      conRed_->applyAdjointHessian_12(*dualcontrol_,*adjoint_,*state_sens_,*state_,z,tol);
      ahwv.plus(*dualcontrol_);
      conRed_->applyAdjointHessian_22(*dualcontrol_,*adjoint_,v,*state_,z,tol);
      ahwv.plus(*dualcontrol_);
    }
  }

// For parametrized (stochastic) objective functions and constraints
public:
  void setParameter(const std::vector<Real> &param) {
    EqualityConstraint<Real>::setParameter(param);
    conVal_->setParameter(param);
    conRed_->setParameter(param);
  }
}; // class Reduced_EqualityConstraint_SimOpt

} // namespace ROL

#endif