ROL_HS24.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef ROL_HS24_HPP
#define ROL_HS24_HPP

#include "ROL_StdVector.hpp"
#include "ROL_TestProblem.hpp"
#include "ROL_Bounds.hpp"

namespace ROL {
namespace ZOO {

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

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

private:
  const Real rt3_;

public:

  Objective_HS24() : rt3_(std::sqrt(3)) {}

  Real value( const Vector<Real> &x, Real &tol ) {

    Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 

    return rt3_*(*xp)[0]*std::pow((*xp)[1],3)*((*xp)[0]-6)/81.0;
  }

  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {

    Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 

    Ptr<vector> gp = dynamic_cast<SV&>(g).getVector(); 

    (*gp)[0] = 2*rt3_*std::pow((*xp)[1],3)*((*xp)[0]-3)/81.0;
    (*gp)[1] = rt3_*(*xp)[0]*std::pow((*xp)[1],2)*((*xp)[0]-6)/27.0;

  }  

  
  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {

    Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector(); 
    Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector(); 
    Ptr<vector> hvp = dynamic_cast<SV&>(hv).getVector(); 

    Real a00 = pow((*xp)[1],3)/81.0;
    Real a01 = pow((*xp)[1],2)*((*xp)[0]-3)/27.0;
    Real a11 = (*xp)[1]*(std::pow((*xp)[0]-3,2)-9)/27.0;

    (*hvp)[0] = a00*(*vp)[0] + a01*(*vp)[1];
    (*hvp)[1] = a01*(*vp)[0] + a11*(*vp)[1]; 
    hv.scale(2*rt3_);

  }
}; // class Objective_HS24


template<class Real> 
class Constraint_HS24 : public Constraint<Real>  {

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

private:
  const Real rt3_;

public:
  Constraint_HS24() : rt3_(std::sqrt(3)) {}

  void value( Vector<Real> &c, const Vector<Real> &x, Real &tol ) {

    Ptr<const vector> xp = dynamic_cast<const SV&>(x).getVector();
    Ptr<vector> cp = dynamic_cast<SV&>(c).getVector();

    (*cp)[0] =  (*xp)[0]/rt3_ -      (*xp)[1];
    (*cp)[1] =  (*xp)[0]      + rt3_*(*xp)[1];
    (*cp)[2] = -(*xp)[0]      - rt3_*(*xp)[1] + 6;

  }

  void applyJacobian( Vector<Real> &jv, const Vector<Real> &v,
                      const Vector<Real> &x, Real &tol ) {

    Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
    Ptr<vector> jvp = dynamic_cast<SV&>(jv).getVector();

    (*jvp)[0] =  (*vp)[0]/rt3_ -      (*vp)[1];
    (*jvp)[1] =  (*vp)[0]      + rt3_*(*vp)[1];
    (*jvp)[2] = -(*vp)[0]      - rt3_*(*vp)[1];


  }

  void applyAdjointJacobian( Vector<Real> &ajv, const Vector<Real> &v, 
                             const Vector<Real> &x, Real &tol ) {

    Ptr<const vector> vp = dynamic_cast<const SV&>(v).getVector();
    Ptr<vector> ajvp = dynamic_cast<SV&>(ajv).getVector();

    (*ajvp)[0] = rt3_*(*vp)[0]/3 + (*vp)[1] - (*vp)[2];
    (*ajvp)[1] = -(*vp)[0] + rt3_*(*vp)[1] - rt3_*(*vp)[2];

  }   

  
  void applyAdjointHessian( Vector<Real> &ahuv, const Vector<Real> &u,
                            const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
    ahuv.zero();
  }

}; // class Constraint_HS24



template<class Real>
class getHS24 : public TestProblem<Real> {
public:
  getHS24(void) {}

  Ptr<Objective<Real> > getObjective( void ) const {
    return makePtr<Objective_HS24<Real>>();
  }

  Ptr<Constraint<Real> > getInequalityConstraint( void ) const {
    return makePtr<Constraint_HS24<Real>>();
  }

  Ptr<BoundConstraint<Real>> getBoundConstraint( void ) const {
    // Lower bound is zero  
    Ptr<std::vector<Real>> lp = makePtr<std::vector<Real>>(2,0.0);
    
    // No upper bound
    Ptr<std::vector<Real>> up = makePtr<std::vector<Real>>(2,ROL_INF<Real>());
   
    Ptr<Vector<Real>> l = makePtr<StdVector<Real>>(lp);
    Ptr<Vector<Real>> u = makePtr<StdVector<Real>>(up);
  
    return makePtr<Bounds<Real>>(l,u);
  }

  Ptr<Vector<Real>> getInitialGuess( void ) const {
    Ptr<std::vector<Real> > x0p = makePtr<std::vector<Real>>(2);
    (*x0p)[0] = 1.0;
    (*x0p)[1] = 0.5;
  
    return makePtr<StdVector<Real>>(x0p);
  }

  Ptr<Vector<Real>> getSolution( const int i = 0 ) const {
    Ptr<std::vector<Real> > xp = makePtr<std::vector<Real>>(2);
    (*xp)[0] = 3.0;
    (*xp)[1] = std::sqrt(3.0);
  
    return makePtr<StdVector<Real>>(xp);
  }

  Ptr<Vector<Real>> getInequalityMultiplier( void ) const {
    Ptr<std::vector<Real> > lp = makePtr<std::vector<Real>>(3,0.0);
    return makePtr<StdVector<Real>>(lp);
  }

  Ptr<BoundConstraint<Real>> getSlackBoundConstraint(void) const {
    // Lower bound is zero  
    Ptr<std::vector<Real> > lp = makePtr<std::vector<Real>>(3,0.0);
    
    // No upper bound
    Ptr<std::vector<Real> > up = makePtr<std::vector<Real>>(3,ROL_INF<Real>());
   
    Ptr<Vector<Real> > l = makePtr<StdVector<Real>>(lp);
    Ptr<Vector<Real> > u = makePtr<StdVector<Real>>(up);
  
    return makePtr<Bounds<Real>>(l,u);
  }
};


} // namespace ZOO
} // namespace ROL

#endif // ROL_HS24_HPP