ROL_Constraint_SerialSimOpt.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


#pragma once
#ifndef ROL_CONSTRAINT_SERIALSIMOPT_HPP
#define ROL_CONSTRAINT_SERIALSIMOPT_HPP

namespace ROL {

#include "ROL_Constraint_TimeSimOpt.hpp"
#include "ROL_PartitionedVector.hpp"
#include "ROL_VectorWorkspace.hpp"


template<typename Real> 
class Constraint_SerialSimOpt : public Constraint_SimOpt<Real> {

  using V      = Vector<Real>;
  using PV     = PartitionedVector<Real>;
  using Con = Constraint_TimeSimOpt<Real>;
 
  using size_type = typename PV<Real>::size_type;

private:

  const Ptr<Con> con_;          // Constraint for single time step
  const Ptr<V>   ui_;           // Initial condition
  Ptr<V> u0_;                   // Zero sim vector on time step
  Pre<V> z0_;                   // Zero opt vector on time step
 
  VectorWorkspace<Real> workspace_; // Memory management

protected:

  PV& partition( V& x ) const { return static_cast<PV&>(x); }

  const V& partition( const V& x ) const { return static_cast<const PV&>(x); }


public:

  Constraint_SerialSimOpt( const Ptr<Con>& con, const V& ui, const V& zi ) :
    Constraint_SimOpt<Real>(), con_(con), 
    ui_( workspace_.copy( ui ) ),
    u0_( workspace_.clone( ui ) ),
    z0_( workspace_.clone( ui ) ) {
    u0_->zero(); z0_->zero();
  }

//  virtual void update_1( const V& u, const V& z, bool flag = true, int iter = -1 ) override {
//    
//  }
//  
//  virtual void update_2( const V& u, const V& z, bool flag = true, int iter = -1 ) override {
//
//  }
  
  virtual void value( V& c, const V& u, const V& z, Real& tol ) override {
  
    auto& cp = partition(c);
    auto& up = partition(u);
    auto& zp = partition(z);
    
    // First interval value uses initial condition
    con_->value( *(cp.get(0)), *ui_, *(zp.get(0)), tol );

    for( size_type k=1, k<cp.numVectors(), ++k ) {
      con_->value( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
    }
  }
  
  virtual void solve( V& c, V& u, const V& z, Real& tol ) override {

    auto& cp = partition(c);
    auto& up = partition(u);
    auto& zp = partition(z);
    
    // First interval value uses initial condition
    con_->solve( *(cp.get(0)), *ui_, *(zp.get(0)), tol );

    for( size_type k=1, k<cp.numVectors(), ++k ) {
      con_->solve( *(cp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
    }
  }
  
  virtual void applyJacobian_1( V& jv, const V& v, const V& u, const V& z, Real& tol ) override {
    
    auto& jvp = partition(jv);      auto& vp = partition(v);
    auto& up  = partition(u);       auto& zp = partition(z);

    auto tmp = workspace_.clone( up.get(0) ); 

    con_->applyJacobian_1_new( *(jvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)), tol );
    
    for( size_type k=1; k<jvp.numVectors(); ++jvp ) {
      con_->applyJacobian_1_new( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
      con_->applyJacobian_1_old( *tmp, *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
      jvp.plus(tmp);
    }
  } // applyJacobian_1
  
  virtual void applyJacobian_2( V& jv, const V& v, const V& u, const V& z, Real& tol)  override { 

    auto& jvp = partition(jv);      auto& vp = partition(v);
    auto& up  = partition(u);       auto& zp = partition(z);

    for( size_type k=0; k<jvp.numVectors(); ++k ) {
      con_->applyJacobian_2( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
    }
  } // applyJacobian_2
  

  virtual void applyInverseJacobian_1( V& ijv, const V& v, const V& u, 
                                       const V& z, Real& tol ) override {

    auto& ijvp = partition(ijv);     auto& vp = partition(v);
    auto& up   = partition(u);       auto& zp = partition(z);

    // First step 
    con_->applyInverseJacobian_1_new( *(ijvp.get(0)), *(vp.get(0), *u0_, *(up.get(0)), *(zp.get(0)), tol );
    
    auto tmp = workspace.clone( vp.get(0) );

    for( size_type k=1; k<ijvp.numVectors(); ++k ) { 

      con_->applyJacobian_1_old( *tmp, *(ijvp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
      tmp->scale(-1.0);
      tmp->plus( *(vp.get(k) );
   
      con_->applyInverseJacobian_1_new( *(ijvp.get(k)), *tmp, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
    }
  } // applyInverseJacobian_1


  
  virtual void applyAdjointJacobian_1( V& ajv, const V& v, const V& u, 
                                       const V& z, Real& tol) override {
 
    auto& ajvp = partition(ajv);     auto& vp = partition(v);
    auto& up   = partition(u);       auto& zp = partition(z);

    auto tmp = workspace.clone( ajvp.get(0) );

    size_type N = ajvp.numVectors()-1;

    // First step
    con_->applyAdjointJacobian_1_new( *(ajvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)) );
    con_->applyAdjointJacobian_1_old( *(ajvp.get(0)), *(vp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)) );

    for( size_type k=1; k<N; ++k ) {
      con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
      con_->applyAdjointJacobian_1_old( *tmp, *(vp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
      ajvp.plus(*tmp);
    }

    // Last step
    con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );

  } // applyAdjointJacobian_1
  


  virtual void applyAdjointJacobian_1( V& ajv, const V& v, const V& u, const V &z, 
                                       const V& dualv, Real& tol) override {

    auto& ajvp = partition(ajv);     auto& vp = partition(dualv);
    auto& up   = partition(u);       auto& zp = partition(z);
    auto tmp = workspace.clone( ajvp.get(0) );

    size_type N = ajvp.numVectors()-1;

    // First step
    con_->applyAdjointJacobian_1_new( *(ajvp.get(0)), *(vp.get(0)), *u0_, *(up.get(0)), *(zp.get(0)) );
    con_->applyAdjointJacobian_1_old( *(ajvp.get(0)), *(vp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)) );

    for( size_type k=1; k<N; ++k ) {
      con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
      con_->applyAdjointJacobian_1_old( *tmp, *(vp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
      ajvp.plus(*tmp);
    }

    // Last step
    con_->applyAdjointJacobian_1_new( *(ajvp.get(k)), *(vp.get(k)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)) );
 
  } // applyAdjointJacobian_1


  
  virtual void applyAdjointJacobian_2( V& ajv, const V& v, const V& u, 
                                       const V& z, Real& tol) override {
    auto& ajvp = partition(ajv);    auto& vp = partition(v);
    auto& up   = partition(u);      auto& zp = partition(z);

    for( size_type k=0; k<jvp.numVectors(); ++k ) {
      con_->applyAdjointJacobian_2_time( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
    }
  } // applyAdjointJacobian_2


  virtual void applyAdjointJacobian_2( V& ajv, const V& v, const V& u, const V& z, 
                                       const V& dualv, Real& tol) override {

    auto& ajvp = partition(ajv);    auto& vp = partition(v);
    auto& up   = partition(u);      auto& zp = partition(z);

    for( size_type k=0; k<jvp.numVectors(); ++k ) {
      con_->applyAdjointJacobian_2_time( *(jvp.get(k)), *(vp.get(k)), *(up.get(k)), *(up.get(k-1)), *(zp.get(k)), tol );
    }
  }


  virtual void applyInverseAdjointJacobian_1( V& iajv, const V& v, const V& u, 
                                              const V& z, Real& tol) override {

    auto& iajvp = partition(iajv);   auto& vp = partition(v);
    auto& up    = partition(u);      auto& zp = partition(z);

    size_type N = iajvp.numVectors()-1;

    // Last Step
    con_->applyInverseAdjointJacobian_1_new( *(iajvp.get(N)), *(vp.get(N)), *(up.get(N-1)), *(up.get(N)), *(zp.get(N)), tol );
    
    auto tmp = workspace_.clone( vp.get(0) );

    for( size_type k=N-1; k>0; --k ) {
      con_->applyAdjointJacobian_1_old( *tmp, *(iajvp.get(k+1)), *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );

      tmp->scale( -1.0 );
      tmp->plus( *(vp.get(k) );

      con_->applyAdjointJacobian_1_new( *(iajvp.get(k)), *tmp, *(up.get(k-1)), *(up.get(k)), *(zp.get(k)), tol );
    }

    // "First step"
    con_->applyAdjointJacobian_1_old( *tmp, *(iajvp.get(1)), *u0_, *(up.get(0)), *(zp.get(0)), tol );

    tmp->scale( -1.0 );
    tmp->plus( *(vp.get(0) );

    con_->applyAdjointJacobian_1_new( *(iajvp.get(0)), *tmp, *u0_, *(up.get(0)), *(zp.get(0)), tol );
  }


  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v, 
                                       const V& u, const V& z, Real& tol) override { 
     // TODO 
  }


  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v, 
                                       const V& u, const V& z, Real& tol) override {
    ahwv.zero();
  }

  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v, 
                                       const V& u, const V& z, Real& tol) override { 
    ahwv.zero();
  }

  virtual void applyAdjointHessian_11( V& ahwv, const V& w, const V& v, 
                                       const V& u, const V& z, Real& tol) override {
    ahwv.zero();
  }


 

   

  
  


}; // class Constraint_SerialSimOpt



} // namespace ROL


#endif // ROL_CONSTRAINT_SERIALSIMOPT_HPP