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

#include "ROL_StdVector.hpp"
#include "Sacado.hpp"

using namespace ROL;

template<class ScalarT>
class FunctionZakharov {

  

  public:

    ScalarT eval(const std::vector<ScalarT> &x);

};

template<class ScalarT>
ScalarT FunctionZakharov<ScalarT>::eval(const std::vector<ScalarT> & x) {

    typedef typename std::vector<ScalarT>::size_type uint;

    ScalarT xdotx = 0;
    ScalarT kdotx = 0;
    ScalarT J = 0;
   
    // Compute dot products 
    for(uint i=0; i<x.size(); ++i) {
        xdotx += pow(x[i],2);       // (k,x)
        kdotx += ScalarT(i+1)*x[i];  // (x,x)
    }

    // Sum terms in objective function
    J = xdotx + pow(kdotx,2)/4.0 + pow(kdotx,4)/16.0;
    
    return J;
}


template<class Real>
class Zakharov_Sacado_Objective : public Objective<Real> {

  typedef std::vector<Real>  vector;
  typedef Vector<Real>       V;
  typedef StdVector<Real>    SV;

  typedef Sacado::Fad::DFad<Real>          GradType;
  typedef Sacado::Fad::SFad<Real,1>        DirDerivType;
  typedef Sacado::Fad::DFad<DirDerivType>  HessVecType;

  typedef typename vector::size_type uint;

  // In C++11, we could use template typedefs:
  // template <typename T>       using  GradTypeT = Sacado::Fad::DFad<T>;
  // typedef Sacado::Fad::SFad<Real,1>  DirDerivType;
  // typedef GradTypeT<Real>            GradType;
  // typedef GradTypeT<DirDerivType>    HessVecType;

  private:

    FunctionZakharov<Real>        zfunc_; 
    FunctionZakharov<GradType>    zfuncGrad_; 
    FunctionZakharov<HessVecType> zfuncHessVec_; 

    ROL::Ptr<const vector> getVector( const V& x ) {
      
      return dynamic_cast<const SV&>(x).getVector();
    }

    ROL::Ptr<vector> getVector( V& x ) {
      
      return dynamic_cast<SV&>(x).getVector();
    }

  public:

    Zakharov_Sacado_Objective() {}

    /* Evaluate the objective function at x */
    Real value( const Vector<Real> &x, Real &tol ) {
      
      ROL::Ptr<const vector> xp = getVector(x);      
      return zfunc_.eval(*xp);
    }

    /* Evaluate the gradient at x */
    void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
      
      ROL::Ptr<const vector> xp = getVector(x);
      ROL::Ptr<vector> gp = getVector(g);

      uint n = xp->size();

      std::vector<GradType> x_grad(n);

      for(uint i=0; i<n; ++i) {
        x_grad[i] = (*xp)[i]; // Set values x(i).
        x_grad[i].diff(i,n);  // Choose canonical directions.
      }

      GradType J_grad = zfuncGrad_.eval(x_grad);

      for(uint i=0; i<n; ++i) {
        (*gp)[i] = J_grad.dx(i);
      }

    }

    /* Compute the action of the Hessian evaluated at x on a vector v */
    void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {

      
      ROL::Ptr<vector> hvp = getVector(hv);
      ROL::Ptr<const vector> vp = getVector(v);
      ROL::Ptr<const vector> xp = getVector(x);

      uint n = xp->size();

      std::vector<HessVecType>  x_hessvec(n);

      for(uint i=0; i<n; ++i) {
        DirDerivType tmp(1,(*xp)[i]);  // Set values x(i).
        tmp.fastAccessDx(0)= (*vp)[i]; // Set direction values v(i).
        x_hessvec[i] = tmp;            // Use tmp to define hessvec-able x.
        x_hessvec[i].diff(i,n);        // Choose directions.
      }

      // Compute Hessian-vector product (and other currently irrelevant things).
      HessVecType J_hessvec = zfuncHessVec_.eval(x_hessvec);

      for(uint i=0; i<n; ++i) {
        (*hvp)[i]   =  (J_hessvec.dx(i)).fastAccessDx(0);
        // hessvec  =  get gradient      get dir deriv
      }
    }
};