step/interiorpoint/test_02.cpp

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

#define OPTIMIZATION_PROBLEM_REFACTOR

#include "Teuchos_GlobalMPISession.hpp"

#include "ROL_RandomVector.hpp"
#include "ROL_StdVector.hpp"
#include "ROL_NonlinearProgram.hpp"
#include "ROL_OptimizationProblem.hpp"
#include "ROL_InteriorPointPenalty.hpp"
#include "ROL_PrimalDualInteriorPointResidual.hpp"

#include "HS_Problem_041.hpp"

#include <iomanip>

template<class Real> 
void printVector( const ROL::Vector<Real> &x, std::ostream &outStream ) {

  try {
    ROL::Ptr<const std::vector<Real> > xp = 
      dynamic_cast<const ROL::StdVector<Real>&>(x).getVector();

    outStream << "Standard Vector" << std::endl;
    for( size_t i=0; i<xp->size(); ++i ) {
      outStream << (*xp)[i] << std::endl;
    }
  }
  catch( const std::bad_cast& e ) {
    outStream << "Partitioned Vector" << std::endl;
    
    typedef ROL::PartitionedVector<Real>    PV;
    typedef typename PV::size_type          size_type;

    const PV &xpv = dynamic_cast<const PV&>(x);

    for( size_type i=0; i<xpv.numVectors(); ++i ) {
      outStream << "--------------------" << std::endl;
      printVector( *(xpv.get(i)), outStream );
    }
    outStream << "--------------------" << std::endl;
  }

}


// Exact residual for H&S Problem 41
template<class Real>
void value( ROL::Vector<Real> &c, const ROL::Vector<Real> &sol, const Real &mu ) {

  typedef std::vector<Real>            vector;
  typedef ROL::StdVector<Real>         SV;
  typedef ROL::PartitionedVector<Real> PV;

  typedef typename PV::size_type size_type;

  
  
  using ROL::dynamicPtrCast;

  const size_type OPT   = 0;
  const size_type EQUAL = 1;
  const size_type LOWER = 2;
  const size_type UPPER = 3;

  const PV &sol_pv = dynamic_cast<const PV&>(sol);
  const vector &x  = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(OPT))->getVector());
  const vector &l  = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(EQUAL))->getVector());
  const vector &zl = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(LOWER))->getVector());
  const vector &zu = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(UPPER))->getVector());

  PV &c_pv = dynamic_cast<PV&>(c);
  vector &cx  = *(ROL::dynamicPtrCast<SV>(c_pv.get(OPT))->getVector());
  vector &cl  = *(ROL::dynamicPtrCast<SV>(c_pv.get(EQUAL))->getVector());
  vector &czl = *(ROL::dynamicPtrCast<SV>(c_pv.get(LOWER))->getVector());
  vector &czu = *(ROL::dynamicPtrCast<SV>(c_pv.get(UPPER))->getVector());
 
  cx[0] = -x[1]*x[2] +   l[0] - zl[0] + zu[0];
  cx[1] = -x[0]*x[1] + 2*l[0] - zl[1] + zu[1];
  cx[2] = -x[0]*x[1] + 2*l[0] - zl[2] + zu[2];
  cx[3] =            -   l[0] - zl[3] + zu[3];

  cl[0] = x[0] + 2*x[1] + 2*x[2] - x[3];

  czl[0] = x[0]*zl[0] - mu;
  czl[1] = x[1]*zl[1] - mu;
  czl[2] = x[2]*zl[2] - mu;
  czl[3] = x[3]*zl[3] - mu;
 
  czu[0] = (1.0-x[0])*zu[0] - mu;
  czu[1] = (1.0-x[1])*zu[1] - mu;
  czu[2] = (1.0-x[2])*zl[2] - mu;
  czu[3] = (2.0-x[3])*zl[3] - mu;
}
 


// Exact residual for H&S Problem 41
template<class Real>
void applyJacobian( ROL::Vector<Real> &jv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &sol ) { 

  typedef std::vector<Real>            vector;
  typedef ROL::StdVector<Real>         SV;
  typedef ROL::PartitionedVector<Real> PV;

  typedef typename PV::size_type size_type;

  
  
  using ROL::dynamicPtrCast;

  const size_type OPT   = 0;
  const size_type EQUAL = 1;
  const size_type LOWER = 2;
  const size_type UPPER = 3;

  const PV &sol_pv = dynamic_cast<const PV&>(sol);
  const vector &x  = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(OPT))->getVector());
//const vector &l  = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(EQUAL))->getVector());
  const vector &zl = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(LOWER))->getVector());
  const vector &zu = *(ROL::dynamicPtrCast<const SV>(sol_pv.get(UPPER))->getVector());

  const PV &v_pv = dynamic_cast<const PV&>(v);
  const vector &vx  = *(ROL::dynamicPtrCast<const SV>(v_pv.get(OPT))->getVector());
  const vector &vl  = *(ROL::dynamicPtrCast<const SV>(v_pv.get(EQUAL))->getVector());
  const vector &vzl = *(ROL::dynamicPtrCast<const SV>(v_pv.get(LOWER))->getVector());
  const vector &vzu = *(ROL::dynamicPtrCast<const SV>(v_pv.get(UPPER))->getVector());

  PV &jv_pv = dynamic_cast<PV&>(jv);
  vector &jvx  = *(ROL::dynamicPtrCast<SV>(jv_pv.get(OPT))->getVector());
  vector &jvl  = *(ROL::dynamicPtrCast<SV>(jv_pv.get(EQUAL))->getVector());
  vector &jvzl = *(ROL::dynamicPtrCast<SV>(jv_pv.get(LOWER))->getVector());
  vector &jvzu = *(ROL::dynamicPtrCast<SV>(jv_pv.get(UPPER))->getVector());
 
  jvx[0] = -x[1]*vx[2] - x[2]*vx[1] +   vl[0] - vzl[0] + vzu[0];
  jvx[1] = -x[0]*vx[2] - x[2]*vx[0] + 2*vl[0] - vzl[1] + vzu[1];
  jvx[2] = -x[0]*vx[1] - x[1]*vx[0] + 2*vl[0] - vzl[2] + vzu[2];
  jvx[3] =                          -   vl[0] - vzl[3] + vzu[3];

  jvl[0] = vx[0] + 2*vx[1] + 2*vx[2] - vx[3];

  jvzl[0] = zl[0]*vx[0] + vzl[0]*x[0];
  jvzl[1] = zl[1]*vx[1] + vzl[1]*x[1];
  jvzl[2] = zl[2]*vx[2] + vzl[2]*x[2];
  jvzl[3] = zl[3]*vx[3] + vzl[3]*x[3];

  jvzu[0] = -zu[0]*vx[0] + vzu[0]*(1.0-x[0]);
  jvzu[1] = -zu[1]*vx[1] + vzu[1]*(1.0-x[1]);
  jvzu[2] = -zu[2]*vx[2] + vzu[2]*(1.0-x[2]);
  jvzu[3] = -zu[3]*vx[3] + vzu[3]*(2.0-x[3]);

}
 

typedef double RealT;

int main(int argc, char *argv[]) {
 
//  typedef std::vector<RealT>                          vector;

  typedef ROL::ParameterList                      PL;

  typedef ROL::Vector<RealT>                          V;
  typedef ROL::PartitionedVector<RealT>               PV;
  typedef ROL::Objective<RealT>                       OBJ;
  typedef ROL::Constraint<RealT>                      CON;
  typedef ROL::BoundConstraint<RealT>                 BND;
  typedef ROL::OptimizationProblem<RealT>             OPT;
  typedef ROL::NonlinearProgram<RealT>                NLP;

  typedef ROL::InteriorPointPenalty<RealT>            PENALTY;
  typedef ROL::PrimalDualInteriorPointResidual<RealT> RESIDUAL;

   

  Teuchos::GlobalMPISession mpiSession(&argc, &argv);

  int iprint = argc - 1;
  ROL::Ptr<std::ostream> outStream;
  ROL::nullstream bhs;
  if( iprint > 0 ) 
    outStream = ROL::makePtrFromRef(std::cout);
  else
    outStream = ROL::makePtrFromRef(bhs);

  int errorFlag = 0;
   
  try {

    RealT mu = 0.1;

    RealT tol = std::sqrt(ROL::ROL_EPSILON<RealT>());

    PL parlist;

    PL &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
    PL &lblist = iplist.sublist("Barrier Objective");

    iplist.set("Symmetrize Primal Dual System",false);

    lblist.set("Use Linear Damping", true);
    lblist.set("Linear Damping Coefficient",1.e-4);
    lblist.set("Initial Barrier Parameter", mu);

    // Need an example problem that satisfies the following criteria:
    // 1) Has an equality constraint 
    // 2) Has a bound constraint where all variables have finite upper and lower bounds

    ROL::Ptr<NLP> nlp = ROL::makePtr<HS::Problem_041<RealT>>(); 
    ROL::Ptr<OPT> opt = nlp->getOptimizationProblem();
 
    ROL::Ptr<V>   x   = opt->getSolutionVector();
    ROL::Ptr<V>   l   = opt->getMultiplierVector();
    ROL::Ptr<V>   zl  = x->clone();
    ROL::Ptr<V>   zu  = x->clone();

    ROL::Ptr<V>   scratch = x->clone();

    ROL::Ptr<PV>  x_pv = ROL::dynamicPtrCast<PV>(x);

    ROL::Ptr<V>   sol = CreatePartitionedVector(x,l,zl,zu);

    ROL::Ptr<V>   c = sol->clone();
    ROL::Ptr<V>   v = sol->clone();
    ROL::Ptr<V>  jv = sol->clone();

    ROL::Ptr<V>   c_exact = c->clone();
    ROL::Ptr<V>  jv_exact = jv->clone();

    ROL::RandomizeVector(*l, -1.0, 1.0);
    ROL::RandomizeVector(*v,  0.0, 1.0);


    ROL::Ptr<OBJ> obj = opt->getObjective();
    ROL::Ptr<CON> con = opt->getConstraint();
    ROL::Ptr<BND> bnd = opt->getBoundConstraint();

    PENALTY penalty(obj,bnd,parlist);
 
    ROL::Ptr<const V> maskL = penalty.getLowerMask();
    ROL::Ptr<const V> maskU = penalty.getUpperMask();

    zl->set(*maskL);
    zu->set(*maskU);

    ROL::Ptr<CON> res = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL,maskU,scratch,mu,false);


    *outStream << "\n[x|lambda|zl|zu]" << std::endl;
    printVector(*sol,*outStream);

    res->value(*c,*sol,tol);
    value(*c_exact,*sol,mu);

    *outStream << "\n[cx|clambda|czl|czu]" << std::endl;

    printVector(*c,*outStream);

    c->axpy(-1.0,*c_exact);

    RealT cerror = c->norm();

    if( cerror>tol ) {
      ++errorFlag;
    }

    *outStream << "\n\n||c-c_exact|| = " << cerror << std::endl;

    res->applyJacobian(*jv,*v,*sol,tol);
    applyJacobian(*jv_exact,*v,*sol);

    *outStream << "\n[vx|vlambda|vzl|vzu]" << std::endl;

    printVector(*v,*outStream);

    *outStream << "\n[jvx|jvlambda|jvzl|jvzu]" << std::endl;

    printVector(*jv,*outStream);

    jv->axpy(-1.0,*jv_exact);

    RealT jverror = jv->norm();
    
    if( jverror > tol ) {
      ++errorFlag;  
    }

    *outStream << "\n\n||jv-jv_exact|| = " << jverror << std::endl;

    *outStream << "Residual Jacobian Finite Difference Check" << std::endl;
    res->checkApplyJacobian(*sol,*v,*v);



  }
  catch (std::logic_error& err) {
    *outStream << err.what() << std::endl;
    errorFlag = -1000;
  }

  if (errorFlag != 0)
    std::cout << "End Result: TEST FAILED" << std::endl;
  else
    std::cout << "End Result: TEST PASSED" << std::endl;

  return 0;
}