test_13.cpp

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_OptimizationSolver.hpp"

#include "ROL_BinaryConstraint.hpp"
#include "ROL_DiagonalOperator.hpp"
#include "ROL_QuadraticObjective.hpp"
#include "ROL_RandomVector.hpp"
#include "ROL_StdVector.hpp"

#include "Teuchos_oblackholestream.hpp"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"



// Creates f(x) = <x,Dx>/2 - <x,b> where D is a diagonal matrix
// If D_{ii}>0 for all i, then the minimizer is the solution to 
// the linear system x_i=b_i/d_i
template<class Real> 
Teuchos::RCP<ROL::Objective<Real>> 
createDiagonalQuadraticObjective( const ROL::Vector<Real> &a, 
                                  const ROL::Vector<Real> &b ) {
  using Teuchos::rcp;
  auto op = rcp( new ROL::DiagonalOperator<Real>(a) );
  auto vec = b.clone();
  vec->set(b);
  auto obj = rcp( new ROL::QuadraticObjective<Real>(op,vec) );
  return obj;
}

typedef double RealT;

int main(int argc, char *argv[]) {
  
//  typedef ROL::Vector<RealT>    V;
  typedef ROL::StdVector<RealT> SV;

  using Teuchos::RCP; using Teuchos::rcp;

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

  // This little trick lets us print to std::cout only if a
  // (dummy) command-line argument is provided.
  int iprint     = argc - 1;
  Teuchos::RCP<std::ostream> outStream;
  Teuchos::oblackholestream bhs; // outputs nothing
  if (iprint > 0)
    outStream = Teuchos::rcp(&std::cout, false);
  else
    outStream = Teuchos::rcp(&bhs, false);

  int errorFlag = 0;

  try { 
 
    auto parlist = Teuchos::getParametersFromXmlFile("binary_constraint.xml");

    // Penalty parameter
    RealT gamma = 1.0;

    RealT INF  = ROL::ROL_INF<RealT>();
    RealT NINF = ROL::ROL_NINF<RealT>();
  
    /*-----  Optimization Vector -----*/

    // Initial guess
    auto x0_rcp = rcp( new std::vector<RealT>(4) );  
    auto x0 = rcp( new SV(x0_rcp) );
    ROL::RandomizeVector(*x0);

    auto x = x0->clone(); x->set(*x0);

    /*----- Objective Function -----*/

    // Diagonal quadratic objective scaling vector
    auto d_rcp = rcp( new std::vector<RealT>(4) );  
    auto d = rcp( new SV(d_rcp) );

    // Quadratic objective offset vector
    auto b_rcp = rcp( new std::vector<RealT>(4) );  
    auto b = rcp( new SV(b_rcp) );

    // Set values for objective
    (*b_rcp)[0] = 1.0;  (*b_rcp)[1] = 1.0;  
    (*b_rcp)[2] = 1.0;  (*b_rcp)[3] = 1.0;   

    (*d_rcp)[0] = 1.0;  (*d_rcp)[1] = 2.0;  
    (*d_rcp)[2] = 4.0;  (*d_rcp)[3] = 8.0;

    auto obj = createDiagonalQuadraticObjective( *d, *b ); 

    // Unconstrained minimizer: x = [ 1.0, 0.5, 0.25, 0.125 ]


    /*----- Bound Constraint -----*/

    // Lower bound vector
    auto xl_rcp = rcp( new std::vector<RealT>(4) );  
    auto xl = rcp( new SV(xl_rcp) );

    // Upper bound vector
    auto xu_rcp = rcp( new std::vector<RealT>(4) );  
    auto xu = rcp( new SV(xu_rcp) );
    
    // Set bounds
    (*xl_rcp)[0] = 0.0;   (*xl_rcp)[1] = 0.0;  
    (*xl_rcp)[2] = NINF;  (*xl_rcp)[3] = NINF;

    (*xu_rcp)[0] = 1.0;   (*xu_rcp)[1] = INF;  
    (*xu_rcp)[2] = 1.0;   (*xu_rcp)[3] = INF;

//    ROL::BoundConstraint<RealT> bnd(xl,xu);
    auto bnd = rcp( new ROL::BoundConstraint<RealT>(xl,xu) );

    /*---- Equality Constraint and Lagrange Multiplier -----*/

    auto con = rcp( new ROL::BinaryConstraint<RealT>( bnd, gamma ) );

    // Lagrange multiplier
    auto l = x->dual().clone();
    ROL::Elementwise::Fill<RealT> FillWithOnes(1.0);
    l->applyUnary( ROL::Elementwise::Fill<RealT>(1.0) );

    // Constrained minimizer set X = { [ 0, 0, 1, 0.125 ], [ 1, 0, 1, 0.125 ] }
   
    // Create Optimization problems
    ROL::OptimizationProblem<RealT> problem_E( obj, x, Teuchos::null, con, l ); 
    ROL::OptimizationProblem<RealT> problem_EB( obj, x, bnd, con, l ); 
   
    // Perform linear algebra and finite difference checks
    problem_E.check( *outStream );
    

    // Solve using Composite Step where the bound is not enforced and
    // equality constraints are satisfied asymptotically
    parlist->sublist("Step").set("Type","Composite Step");
    ROL::OptimizationSolver<RealT> solver_cs( problem_E, *parlist );
    solver_cs.solve( *outStream );    
    *outStream << "\n\nFinal optimization vector:";
    x->print(*outStream);

    // Reset optimization vector and Lagrange multiplier to initial values 
    x->set(*x0); l->applyUnary(FillWithOnes);
  

    // Solve using Augmented Lagrangian where the bound is enforced explicitly 
    // and equality constraints are enforced through penalization
    parlist->sublist("Step").set("Type","Augmented Lagrangian");
    ROL::OptimizationSolver<RealT> solver_al( problem_EB, *parlist );
    solver_al.solve( *outStream );    
    *outStream << "\n\nFinal optimization vector:";
    x->print(*outStream);
    

    // Reset optimization vector and Lagrange multiplier to initial values 
    x->set(*x0); l->applyUnary(FillWithOnes);


    // Solve using Moreau-Yosida where the bound is enforced through penalization 
    // and equality constraints are satisfied asymptotically
    parlist->sublist("Step").set("Type","Moreau-Yosida Penalty");
    ROL::OptimizationSolver<RealT> solver_my( problem_EB, *parlist );
    solver_my.solve( *outStream );    
    *outStream << "\n\nFinal optimization vector:";
    x->print(*outStream);
    

  }
  
  catch (std::logic_error err) {
    *outStream << err.what() << "\n";
    errorFlag = -1000;
  }; // end try

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

  return 0;
}