diff --git a/packages/rol/CMakeLists.txt b/packages/rol/CMakeLists.txt index ee8430a96c21..98c276b8c796 100644 --- a/packages/rol/CMakeLists.txt +++ b/packages/rol/CMakeLists.txt @@ -24,6 +24,12 @@ TRIBITS_ADD_OPTION_AND_DEFINE(${PACKAGE_NAME}_ENABLE_TIMERS OFF ) +TRIBITS_ADD_OPTION_AND_DEFINE(${PACKAGE_NAME}_ENABLE_PYROL + ENABLE_PYBIND11_PYROL + "Build ROL with PyROL interface." + OFF + ) + # Build Options SET( CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR}/cmake) INCLUDE(BuildOptions) @@ -49,6 +55,10 @@ IF (ROL_ENABLE_ArrayFireCPU) ADD_SUBDIRECTORY(adapters/arrayfire) ENDIF() +IF( ROL_ENABLE_PYROL ) + ADD_SUBDIRECTORY(pyrol) +ENDIF() + # Adapters which require use of Teuchos #GET_PROPERTY( PTR_STRING GLOBAL PROPERTY PTR_IMPL ) diff --git a/packages/rol/cmake/PythonTest.cmake b/packages/rol/cmake/PythonTest.cmake new file mode 100644 index 000000000000..2568fa4b6377 --- /dev/null +++ b/packages/rol/cmake/PythonTest.cmake @@ -0,0 +1,22 @@ +MACRO(PYTHON_ADD_TEST TEST_NAME) + + TRIBITS_ADD_TEST( + ${PYTHON_EXECUTABLE} + NOEXEPREFIX + NOEXESUFFIX + NAME PYTHON_${TEST_NAME} + ARGS "-m unittest ${TEST_NAME}.py" + NUM_MPI_PROCS 1 + PASS_REGULAR_EXPRESSION "OK" + ${ARGN} + ) + + IF(TPL_ENABLE_MPI) + tribits_set_tests_properties(ROL_PYTHON_${TEST_NAME}_MPI_1 + PROPERTIES ENVIRONMENT "${PyROL_PYTHONPATH}") + ELSE() + tribits_set_tests_properties(ROL_PYTHON_${TEST_NAME} + PROPERTIES ENVIRONMENT "${PyROL_PYTHONPATH}") + ENDIF() + +ENDMACRO(PYTHON_ADD_TEST TEST_NAME) diff --git a/packages/rol/cmake/ROL_config.h.in b/packages/rol/cmake/ROL_config.h.in index f267fb4fac7e..80df89cfa5fd 100644 --- a/packages/rol/cmake/ROL_config.h.in +++ b/packages/rol/cmake/ROL_config.h.in @@ -51,3 +51,6 @@ /* Define for performance timer support. */ #cmakedefine ROL_TIMERS + +/* Define for python interface support. */ +#cmakedefine ENABLE_PYBIND11_PYROL diff --git a/packages/rol/example/PDE-OPT/dynamic/adv_diff/CMakeLists.txt b/packages/rol/example/PDE-OPT/dynamic/adv_diff/CMakeLists.txt index 26a273f1a5f1..f34b376b327d 100644 --- a/packages/rol/example/PDE-OPT/dynamic/adv_diff/CMakeLists.txt +++ b/packages/rol/example/PDE-OPT/dynamic/adv_diff/CMakeLists.txt @@ -19,6 +19,11 @@ IF(${PROJECT_NAME}_ENABLE_Intrepid AND SOURCES example_02.cpp ADD_DIR_TO_NAME ) + TRIBITS_ADD_EXECUTABLE( + example_03 + SOURCES example_03.cpp + ADD_DIR_TO_NAME + ) TRIBITS_COPY_FILES_TO_BINARY_DIR( ParabolicDataCopy SOURCE_FILES diff --git a/packages/rol/example/PDE-OPT/dynamic/adv_diff/example_03.cpp b/packages/rol/example/PDE-OPT/dynamic/adv_diff/example_03.cpp new file mode 100644 index 000000000000..2b0850eb5a50 --- /dev/null +++ b/packages/rol/example/PDE-OPT/dynamic/adv_diff/example_03.cpp @@ -0,0 +1,283 @@ +// @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 + +/*! \file example_01.cpp + \brief Solves a source inversion problem governed by the + advection-diffusion equation. +*/ + +#include "Teuchos_Comm.hpp" +#include "Teuchos_GlobalMPISession.hpp" +#include "Tpetra_Core.hpp" +#include "Tpetra_Version.hpp" + +#include "ROL_Stream.hpp" +#include "ROL_ParameterList.hpp" +#include "ROL_Solver.hpp" +#include "ROL_ReducedDynamicObjective.hpp" +#include "ROL_Bounds.hpp" +#include "ROL_DynamicConstraintCheck.hpp" +#include "ROL_DynamicObjectiveCheck.hpp" +#include "ROL_TypeBIndicatorObjective.hpp" +#include +//#include + +#include "../../TOOLS/meshmanager.hpp" +#include "../../TOOLS/lindynconstraint.hpp" +#include "../../TOOLS/ltiobjective.hpp" +#include "../../TOOLS/pdevector.hpp" +#include "../../TOOLS/pdeobjective.hpp" +#include "dynpde_adv_diff.hpp" +#include "obj_adv_diff.hpp" +#include "mesh_adv_diff.hpp" +#include "l1penaltydynamic.hpp" + +#include "ROL_TypeP_TrustRegionAlgorithm.hpp" + + +int main(int argc, char *argv[]) { +// feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); + + using RealT = double; + + /*** Initialize communicator. ***/ + Teuchos::GlobalMPISession mpiSession(&argc, &argv); + ROL::Ptr> comm + = Tpetra::getDefaultComm(); + + // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided. + const int myRank = comm->getRank(); + ROL::Ptr outStream = ROL::makeStreamPtr( std::cout, (argc > 1) && (myRank==0) ); + + int errorFlag = 0; + + // *** Example body. + try { + + /*** Read in XML input ***/ + ROL::Ptr parlist = ROL::getParametersFromXmlFile("input.xml"); + int nt = parlist->sublist("Time Discretization").get("Number of Time Steps", 100); + RealT T = parlist->sublist("Time Discretization").get("End Time", 1.0); + RealT dt = T/static_cast(nt); + int controlDim = 9; + + /*************************************************************************/ + /***************** BUILD GOVERNING PDE ***********************************/ + /*************************************************************************/ + /*** Initialize main data structure. ***/ + ROL::Ptr> meshMgr + = ROL::makePtr>(*parlist); + // Initialize PDE describing advection-diffusion equation + ROL::Ptr> pde + = ROL::makePtr>(*parlist); + + /*************************************************************************/ + /***************** BUILD CONSTRAINT **************************************/ + /*************************************************************************/ + bool isLTI = !parlist->sublist("Problem").get("Time Varying Coefficients",false); + ROL::Ptr> dyn_con + = ROL::makePtr>(pde,meshMgr,comm,*parlist,isLTI,*outStream); + dyn_con->getAssembler()->printMeshData(*outStream); + + /*************************************************************************/ + /***************** BUILD STATE VECTORS ***********************************/ + /*************************************************************************/ + ROL::Ptr> u0_ptr, uo_ptr, un_ptr, ck_ptr; + u0_ptr = dyn_con->getAssembler()->createStateVector(); + uo_ptr = dyn_con->getAssembler()->createStateVector(); + un_ptr = dyn_con->getAssembler()->createStateVector(); + ck_ptr = dyn_con->getAssembler()->createResidualVector(); + ROL::Ptr> u0, uo, un, ck, zk; + u0 = ROL::makePtr>(u0_ptr,pde,*dyn_con->getAssembler()); + uo = ROL::makePtr>(uo_ptr,pde,*dyn_con->getAssembler()); + un = ROL::makePtr>(un_ptr,pde,*dyn_con->getAssembler()); + ck = ROL::makePtr>(ck_ptr,pde,*dyn_con->getAssembler()); + zk = ROL::makePtr>(ROL::makePtr>(controlDim)); + ROL::Ptr> z + = ROL::PartitionedVector::create(*zk, nt); + + /*************************************************************************/ + /***************** BUILD COST FUNCTIONAL *********************************/ + /*************************************************************************/ + std::vector>> qoi_vec(2,ROL::nullPtr); + qoi_vec[0] = ROL::makePtr>(pde->getFE()); + qoi_vec[1] = ROL::makePtr>(pde->getFE()); + RealT stateCost = parlist->sublist("Problem").get("State Cost",1e5); + RealT myCost = 1e0; + std::vector wts = {stateCost, myCost}; + ROL::Ptr> obj_k + = ROL::makePtr>(qoi_vec,wts,dyn_con->getAssembler()); + ROL::Ptr> dyn_obj + = ROL::makePtr>(*parlist,obj_k,false); + + /*************************************************************************/ + /***************** BUILD REDUCED COST FUNCTIONAL *************************/ + /*************************************************************************/ + std::vector> timeStamp(nt); + for( int k=0; ksublist("Reduced Dynamic Objective"); + ROL::Ptr> obj + = ROL::makePtr>(dyn_obj, dyn_con, u0, zk, ck, timeStamp, rpl, outStream); + + /*************************************************************************/ + /***************** BUILD BOUND CONSTRAINT AND L1 PENALTY *****************/ + /*************************************************************************/ + ROL::Ptr> zlo = ROL::PartitionedVector::create(*zk, nt); + ROL::Ptr> zhi = ROL::PartitionedVector::create(*zk, nt); + zlo->setScalar(-static_cast(1)); zhi->setScalar(static_cast(1)); + ROL::Ptr> nobj + = ROL::makePtr>(*parlist, timeStamp, zlo, zhi); + + /*************************************************************************/ + /***************** RUN VECTOR AND DERIVATIVE CHECKS **********************/ + /*************************************************************************/ + bool checkDeriv = parlist->sublist("Problem").get("Check Derivatives",false); + if ( checkDeriv ) { + ROL::Ptr> dz = ROL::PartitionedVector::create(*zk, nt); + ROL::Ptr> hz = ROL::PartitionedVector::create(*zk, nt); + z->randomize(); dz->randomize(); hz->randomize(); + zk->randomize(); uo->randomize(); un->randomize(); + ROL::ValidateFunction validate(1,13,20,11,true,*outStream); + ROL::DynamicObjectiveCheck::check(*dyn_obj,validate,*uo,*un,*zk,timeStamp[0]); + ROL::DynamicConstraintCheck::check(*dyn_con,validate,*uo,*un,*zk,timeStamp[0]); + obj->checkGradient(*z,*dz,true,*outStream); + obj->checkHessVec(*z,*dz,true,*outStream); + obj->checkHessSym(*z,*dz,*hz,true,*outStream); + } + + /*************************************************************************/ + /***************** SOLVE OPTIMIZATION PROBLEM ****************************/ + /*************************************************************************/ + RealT gtol(1e-2), zerr(0), ztol0(1e-6), ztol(1); + z->zero(); + ROL::Ptr> zprev = z->clone(), zdiff = z->clone(); + ROL::Ptr> algo; + parlist->sublist("Status Test").set("Use Relative Tolerances",false); + for (unsigned i = 0; i < 20; ++i) { + *outStream << std::endl << "Moreau-Yosida Parameter: " + << myCost << std::endl + << "Gradient Tolerance: " + << gtol << std::endl << std::endl; + + ztol = ztol0 / std::max(static_cast(1),z->norm()); + zprev->set(*z); + parlist->sublist("Status Test").set("Gradient Tolerance",gtol); + parlist->sublist("Status Test").set("Step Tolerance",static_cast(1e-4)*gtol); + algo = ROL::makePtr>(*parlist); + std::clock_t timer = std::clock(); + algo->run(*z, *obj, *nobj, *outStream); + *outStream << "Optimization time: " + << static_cast(std::clock()-timer)/static_cast(CLOCKS_PER_SEC) + << " seconds." << std::endl << std::endl; + parlist->sublist("Reduced Dynamic Objective").set("State Rank", static_cast(obj->getStateRank())); + parlist->sublist("Reduced Dynamic Objective").set("Adjoint Rank", static_cast(obj->getAdjointRank())); + parlist->sublist("Reduced Dynamic Objective").set("State Sensitivity Rank", static_cast(obj->getStateSensitivityRank())); + + zdiff->set(*z); zdiff->axpy(static_cast(-1),*zprev); + zerr = zdiff->norm(); + *outStream << std::endl << "Control Difference: " + << zerr << std::endl << std::endl; + if (zerr < ztol && algo->getState()->gnorm < static_cast(1e-8)) break; + + myCost *= static_cast(2e0); + gtol *= static_cast(1e-1); gtol = std::max(gtol,static_cast(1e-10)); + wts = {stateCost, myCost}; + obj_k = ROL::makePtr>(qoi_vec,wts,dyn_con->getAssembler()); + dyn_obj = ROL::makePtr>(*parlist,obj_k,false); + obj = ROL::makePtr>(dyn_obj, dyn_con, u0, zk, ck, timeStamp, rpl, outStream); + } + + /*************************************************************************/ + /***************** OUTPUT RESULTS ****************************************/ + /*************************************************************************/ + std::clock_t timer_print = std::clock(); + // Output control to file + if ( myRank == 0 ) { + ROL::Ptr> zn; + std::ofstream zfile; + zfile.open("control.txt"); + zfile << std::scientific << std::setprecision(15); + for (int n = 0; n < nt; ++n) { + zn = ROL::dynamicPtrCast>(z->get(n))->getParameter()->getVector(); + for (int i = 0; i < controlDim; ++i) { + zfile << std::right << std::setw(25) << (*zn)[i]; + } + zfile << std::endl; + } + zfile.close(); + } + // Output state to file + uo->set(*u0); un->zero(); + for (int k = 1; k < nt; ++k) { + // Print previous state to file + std::stringstream ufile; + ufile << "state." << k-1 << ".txt"; + dyn_con->outputTpetraVector(uo_ptr, ufile.str()); + // Advance time stepper + dyn_con->solve(*ck, *uo, *un, *z->get(k), timeStamp[k]); + uo->set(*un); + } + // Print previous state to file + std::stringstream ufile; + ufile << "state." << nt-1 << ".txt"; + dyn_con->outputTpetraVector(uo_ptr, ufile.str()); + *outStream << "Output time: " + << static_cast(std::clock()-timer_print)/static_cast(CLOCKS_PER_SEC) + << " seconds." << std::endl << std::endl; + } + 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; +} diff --git a/packages/rol/example/PDE-OPT/dynamic/adv_diff/l1penaltydynamic.hpp b/packages/rol/example/PDE-OPT/dynamic/adv_diff/l1penaltydynamic.hpp index c77096f0d966..6e395e27cded 100644 --- a/packages/rol/example/PDE-OPT/dynamic/adv_diff/l1penaltydynamic.hpp +++ b/packages/rol/example/PDE-OPT/dynamic/adv_diff/l1penaltydynamic.hpp @@ -13,12 +13,11 @@ class L1_Dyn_Objective : public ROL::Objective { using size_type = typename std::vector::size_type; private: Real theta_, beta_;//, T_; - const size_type Nt_; + const size_type Nt_; const std::vector> ts_; - const ROL::Ptr> zl_; - const ROL::Ptr> zu_; - - ROL::PartitionedVector &partition ( ROL::Vector& x ) const { + const ROL::Ptr> zl_, zu_; + + ROL::PartitionedVector &partition ( ROL::Vector& x ) const { return static_cast&>(x); } @@ -26,7 +25,7 @@ class L1_Dyn_Objective : public ROL::Objective { return static_cast&>(x); } - ROL::Ptr > getConstParameter(const ROL::Vector &x) const { + ROL::Ptr > getConstParameter(const ROL::Vector &x) const { ROL::Ptr > xp; try { xp = dynamic_cast&>(x).getVector(); @@ -72,104 +71,93 @@ class L1_Dyn_Objective : public ROL::Objective { return xp; } - - public: - L1_Dyn_Objective(ROL::ParameterList &parlist, - const std::vector> &timeStamp, - const ROL::Ptr> &zl, - const ROL::Ptr> &zu - ) : Nt_(timeStamp.size()), ts_(timeStamp), zl_(zl), zu_(zu) - { +public: + L1_Dyn_Objective(ROL::ParameterList &parlist, + const std::vector> &timeStamp, + const ROL::Ptr> &zl, + const ROL::Ptr> &zu) + : Nt_(timeStamp.size()), ts_(timeStamp), zl_(zl), zu_(zu) { theta_ = parlist.sublist("Reduced Dynamic Objective").sublist("Time Discretization").get("Theta", 1.0); - beta_ = parlist.sublist("Problem").get("Control Cost", 1.e0); - - } + beta_ = parlist.sublist("Problem").get("Control Cost", 1.e0); + } - // value - Real value(const ROL::Vector &z, - Real &tol - ){ + // value + Real value(const ROL::Vector &z,Real &tol) { + const ROL::PartitionedVector &zp = partition(z); + const Real one(1); + Real dt(0), val(0); + bool isinf(false); + + for (size_type k = 0; k < Nt_; ++k){ // dynamic obj + ROL::Ptr> zlk = getConstParameter(*zl_->get(k)); + ROL::Ptr> zuk = getConstParameter(*zu_->get(k)); + ROL::Ptr> zpk = getConstParameter(*zp.get(k)); // ROL vector of partition/kth timestep + for (size_type i = 0; i < zpk->size(); ++i){ + if ((*zpk)[i] < (*zlk)[i] || (*zpk)[i]>(*zuk)[i]){ + val = ROL::ROL_INF(); + isinf = true; + break; + } + } + } - const ROL::PartitionedVector &zp = partition(z); - const Real one(1); - Real dt(0), val(0); - bool isinf(false); - - for (size_type k = 0; k < Nt_; ++k){//dynamic obj - ROL::Ptr> zlk = getConstParameter(*zl_->get(k)); - ROL::Ptr> zuk = getConstParameter(*zu_->get(k)); - ROL::Ptr> zpk = getConstParameter(*zp.get(k)); // ROL vector of partition/kth timestep - - for (size_type i = 0; i < zpk->size(); ++i){ - if ((*zpk)[i] < (*zlk)[i] || (*zpk)[i]>(*zuk)[i]){ - val = ROL::ROL_INF(); - isinf=true; - break; - } - } - } - - if (!isinf) { - for (size_type k = 1; k < Nt_; ++k){//dynamic obj - - ROL::Ptr> zpn = getConstParameter(*zp.get(k)); // ROL vector of partition/kth timestep - ROL::Ptr> zpo = getConstParameter(*zp.get(k-1)); // ROL vector of partition/kth timestep - dt = ts_[k].t[1] - ts_[k].t[0]; - - for (size_type i = 0; i < zpn->size(); ++i){ - val += dt*((one - theta_)*std::abs((*zpo)[i]) + theta_*std::abs((*zpn)[i])); - } // end i for - }// end k for - } - return beta_*val; - } //end value - - // prox - void prox(ROL::Vector &Pz, - const ROL::Vector &z, - Real t, - Real &tol ) - { - - //partitioned vectors - const ROL::PartitionedVector &zp = partition(z); + if (!isinf) { + for (size_type k = 1; k < Nt_; ++k){ // dynamic obj + ROL::Ptr> zpn = getConstParameter(*zp.get(k)); // ROL vector of partition/kth timestep + ROL::Ptr> zpo = getConstParameter(*zp.get(k-1)); // ROL vector of partition/kth timestep + dt = ts_[k].t[1] - ts_[k].t[0]; + for (size_type i = 0; i < zpn->size(); ++i){ + val += dt*((one - theta_)*std::abs((*zpo)[i]) + theta_*std::abs((*zpn)[i])); + } // end i for + } // end k for + } + return beta_*val; + } // end value + + // prox + void prox(ROL::Vector &Pz, + const ROL::Vector &z, + Real t, + Real &tol) { + //partitioned vectors + const ROL::PartitionedVector &zp = partition(z); ROL::PartitionedVector &Pzp = partition(Pz); - //constants - const Real one(1), zero(0); - Real hk(0), hkplus(0); - Real l1param(0); + //constants + const Real one(1), zero(0); + Real hk(0), hkplus(0); + Real l1param(0); - for (size_type k = 0; kNt inclusive + for (size_type k = 0; kNt inclusive ROL::Ptr> zlk = getConstParameter(*zl_->get(k)); - ROL::Ptr> zuk = getConstParameter(*zu_->get(k)); - ROL::Ptr> zk = getConstParameter(zp[k]); //kth timestep - ROL::Ptr> Pzk = getParameter(Pzp[k]); + ROL::Ptr> zuk = getConstParameter(*zu_->get(k)); + ROL::Ptr> zk = getConstParameter(zp[k]); //kth timestep + ROL::Ptr> Pzk = getParameter(Pzp[k]); - // update prox parameter - if (k == 0){ + // update prox parameter + if (k == 0){ hkplus = ts_[k+1].t[1] - ts_[k+1].t[0]; - l1param = t*hkplus*beta_*(one - theta_); - } else if (k == (Nt_ - 1)) { - hk = ts_[k].t[1] - ts_[k].t[0]; + l1param = t*hkplus*beta_*(one - theta_); + } + else if (k == (Nt_ - 1)) { + hk = ts_[k].t[1] - ts_[k].t[0]; l1param = t*hk*beta_*theta_; - } else { - hk = ts_[k].t[1] - ts_[k].t[0]; - hkplus = ts_[k+1].t[1] - ts_[k+1].t[0]; + } + else { + hk = ts_[k].t[1] - ts_[k].t[0]; + hkplus = ts_[k+1].t[1] - ts_[k+1].t[0]; l1param = t*beta_*(hkplus*(one - theta_) + theta_*hk); - } + } - for (size_type i = 0; i < zk->size(); ++i){ - if ((*zk)[i] > l1param) { (*Pzk)[i] = (*zk)[i] - l1param; } - else if ((*zk)[i] < -l1param) {(*Pzk)[i] = (*zk)[i] + l1param;} - else { (*Pzk)[i] = zero;} - - (*Pzk)[i] = std::min((*zuk)[i], std::max((*zlk)[i], (*Pzk)[i])); - } // end i loop - - }//end k for - - }//end prox + for (size_type i = 0; i < zk->size(); ++i){ + if ((*zk)[i] > l1param) { (*Pzk)[i] = (*zk)[i] - l1param; } + else if ((*zk)[i] < -l1param) {(*Pzk)[i] = (*zk)[i] + l1param;} + else { (*Pzk)[i] = zero;} + + (*Pzk)[i] = std::min((*zuk)[i], std::max((*zlk)[i], (*Pzk)[i])); + } // end i loop + } // end k for + } // end prox }; #endif diff --git a/packages/rol/example/PDE-OPT/dynamic/adv_diff/obj_adv_diff.hpp b/packages/rol/example/PDE-OPT/dynamic/adv_diff/obj_adv_diff.hpp index 547b8a78d941..0a5cd8e4bc07 100644 --- a/packages/rol/example/PDE-OPT/dynamic/adv_diff/obj_adv_diff.hpp +++ b/packages/rol/example/PDE-OPT/dynamic/adv_diff/obj_adv_diff.hpp @@ -201,6 +201,205 @@ class QoI_State_Cost_adv_diff : public QoI { }; // QoI_State_Cost +template +class QoI_State_MY_adv_diff : public QoI { +private: + ROL::Ptr> fe_; + ROL::Ptr> ubnd_; + + Real computeUpperBound(const std::vector &x) const { + return static_cast(0.04); + } + +public: + QoI_State_MY_adv_diff(const ROL::Ptr> &fe) : fe_(fe) { + int c = fe_->cubPts()->dimension(0); + int p = fe_->cubPts()->dimension(1); + int d = fe_->cubPts()->dimension(2); + std::vector pt(d); + ubnd_ = ROL::makePtr>(c,p); + for (int i = 0; i < c; ++i) { + for (int j = 0; j < p; ++j) { + for (int k = 0; k < d; ++k) { + pt[k] = (*fe_->cubPts())(i,j,k); + } + (*ubnd_)(i,j) = computeUpperBound(pt); + } + } + } + + Real value(ROL::Ptr> & val, + const ROL::Ptr> & u_coeff, + const ROL::Ptr> & z_coeff = ROL::nullPtr, + const ROL::Ptr> & z_param = ROL::nullPtr) { + // Get relevant dimensions + const int c = fe_->gradN()->dimension(0); + const int p = fe_->gradN()->dimension(2); + // Initialize output val + val = ROL::makePtr>(c); + // Evaluate state on FE basis + ROL::Ptr> valU_eval = + ROL::makePtr>(c, p); + fe_->evaluateValue(valU_eval, u_coeff); + // Compute integrand + ROL::Ptr> minU_eval = + ROL::makePtr>(c, p); + for (int i = 0; i < c; ++i) { + for (int j = 0; j < p; ++j) { + (*minU_eval)(i,j) = std::max(static_cast(0),(*valU_eval)(i,j)-(*ubnd_)(i,j)); + } + } + // Compute squared L2-norm of diff + fe_->computeIntegral(val,minU_eval,minU_eval); + // Scale by one half + Intrepid::RealSpaceTools::scale(*val,static_cast(0.5)); + return static_cast(0); + } + + void gradient_1(ROL::Ptr > & grad, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + // Get relevant dimensions + const int c = fe_->gradN()->dimension(0); + const int f = fe_->gradN()->dimension(1); + const int p = fe_->gradN()->dimension(2); + // Initialize output grad + grad = ROL::makePtr>(c, f); + // Evaluate state on FE basis + ROL::Ptr> valU_eval = + ROL::makePtr>(c, p); + fe_->evaluateValue(valU_eval, u_coeff); + // Compute integrand + ROL::Ptr> minU_eval = + ROL::makePtr>(c, p); + for (int i = 0; i < c; ++i) { + for (int j = 0; j < p; ++j) { + (*minU_eval)(i,j) = std::max(static_cast(0),(*valU_eval)(i,j)-(*ubnd_)(i,j)); + } + } + // Compute gradient of squared L2-norm of diff + Intrepid::FunctionSpaceTools::integrate(*grad, + *minU_eval, + *(fe_->NdetJ()), + Intrepid::COMP_CPP, false); + } + + void gradient_2(ROL::Ptr > & grad, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::gradient_2 is zero."); + } + + std::vector gradient_3(std::vector > > & grad, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::gradient_3 is zero."); + } + + void HessVec_11(ROL::Ptr > & hess, + const ROL::Ptr > & v_coeff, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + const int c = fe_->gradN()->dimension(0); + const int f = fe_->gradN()->dimension(1); + const int p = fe_->gradN()->dimension(2); + ROL::Ptr> valU_eval = + ROL::makePtr>(c, p); + ROL::Ptr> valV_eval = + ROL::makePtr>(c, p); + hess = ROL::makePtr>(c, f); + fe_->evaluateValue(valU_eval, u_coeff); + fe_->evaluateValue(valV_eval, v_coeff); + // Compute integrand + ROL::Ptr> minV_eval = + ROL::makePtr>(c, p); + for (int i = 0; i < c; ++i) { + for (int j = 0; j < p; ++j) { + if ( (*valU_eval)(i,j) > (*ubnd_)(i,j) ) { + (*minV_eval)(i,j) = (*valV_eval)(i,j); + } + else { + (*minV_eval)(i,j) = static_cast(0); + } + } + } + Intrepid::FunctionSpaceTools::integrate(*hess, + *minV_eval, + *(fe_->NdetJ()), + Intrepid::COMP_CPP, false); + } + + void HessVec_12(ROL::Ptr > & hess, + const ROL::Ptr > & v_coeff, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_12 is zero."); + } + + void HessVec_13(ROL::Ptr > & hess, + const ROL::Ptr > & v_param, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_13 is zero."); + } + + void HessVec_21(ROL::Ptr > & hess, + const ROL::Ptr > & v_coeff, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_21 is zero."); + } + + void HessVec_22(ROL::Ptr > & hess, + const ROL::Ptr > & v_coeff, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_22 is zero."); + } + + void HessVec_23(ROL::Ptr > & hess, + const ROL::Ptr > & v_param, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_23 is zero."); + } + + std::vector HessVec_31(std::vector > > & hess, + const ROL::Ptr > & v_coeff, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_31 is zero."); + } + + std::vector HessVec_32(std::vector > > & hess, + const ROL::Ptr > & v_coeff, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_32 is zero."); + } + + std::vector HessVec_33(std::vector > > & hess, + const ROL::Ptr > & v_param, + const ROL::Ptr > & u_coeff, + const ROL::Ptr > & z_coeff = ROL::nullPtr, + const ROL::Ptr > & z_param = ROL::nullPtr) { + throw Exception::Zero(">>> QoI_State_MY_adv_diff::HessVec_33 is zero."); + } + +}; // QoI_State_Cost + + template class QoI_Control_Cost_adv_diff : public QoI { public: diff --git a/packages/rol/pyrol/CMakeLists.txt b/packages/rol/pyrol/CMakeLists.txt new file mode 100644 index 000000000000..be85d4bf11bf --- /dev/null +++ b/packages/rol/pyrol/CMakeLists.txt @@ -0,0 +1,457 @@ +# -*- cmake -*- + +# @HEADER +# *********************************************************************** +# +# PyROL: Automatic Python Interfaces to ROL +# Copyright (2022) Sandia Corporation +# +# Under the terms of Contract DE-AC04-94AL85000 with Sandia +# Corporation, the U.S. Government retains certain rights in this +# software. +# +# 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 Kim Liegeois (knliege@sandia.gov) +# +# *********************************************************************** +# @HEADER + + +IF(NOT BUILD_SHARED_LIBS) + MESSAGE(FATAL_ERROR "PyROL can only be built with shared libraries. Building of shared libraries is currently set to OFF. To enable shared libraries please set the cache variable \"BUILD_SHARED_LIBS\" to ON") +ENDIF() + +# Set the package version number +SET(PyROL_VERSION ${Trilinos_VERSION}) + +TRIBITS_ADD_OPTION_AND_DEFINE(PYROL_PIP_INSTALL + PYROL_SCIKIT + "Logic specific to the scikit-build-core backend." + OFF ) + +TRIBITS_ADD_OPTION_AND_DEFINE(PYROL_ENABLE_BINDER + PYROL_GENERATE_SRC + "Enable regeneration of source files with Binder." + OFF ) + +TRIBITS_ADD_OPTION_AND_DEFINE(PYROL_BINDER_SUPPRESS_ERRORS + PYROL_SUPPRESS_ERRORS + "Enable the suppress errors option of Binder." + OFF ) + +TRIBITS_ADD_OPTION_AND_DEFINE(PYROL_BINDER_CMAKE_ERROR + PYROL_CMAKE_ERROR + "Stop the configuration if binder fails." + ON ) + +TRIBITS_ADD_OPTION_AND_DEFINE(PYROL_BINDER_VERBOSE + PYROL_B_VERBOSE + "Increase the verbosity of binder" + OFF ) + +TRIBITS_ADD_OPTION_AND_DEFINE(PYROL_ENABLE_BINDER_UPDATE + PYROL_UPDATE_GENERATED_SRC + "Enable the update of the generated source files with Binder." + OFF ) + +MESSAGE("-- PYTHON_EXECUTABLE:") +IF(NOT DEFINED ${PYTHON_EXECUTABLE}) + find_program(PYTHON_EXECUTABLE + NAMES python3 python + ) + MESSAGE(" -- CMake has set: PYTHON_EXECUTABLE = ${PYTHON_EXECUTABLE}") +ELSE() + MESSAGE(" -- User has set: PYTHON_EXECUTABLE = ${PYTHON_EXECUTABLE}") +ENDIF() + + +function(get_all_include_dirs LIBRARY_NAME all_include_dirs all_visited_libs) + if (TARGET ${LIBRARY_NAME}) + get_property(depend_libs TARGET ${LIBRARY_NAME} PROPERTY INTERFACE_LINK_LIBRARIES) + foreach(depend_lib IN LISTS depend_libs) + if (TARGET ${depend_lib} AND (NOT ${depend_lib} IN_LIST all_visited_libs)) + list(APPEND all_visited_libs "${depend_lib}") # Update list in the current scope only + get_property(current_includes TARGET ${depend_lib} PROPERTY INTERFACE_INCLUDE_DIRECTORIES) + foreach(include IN LISTS current_includes) + STRING(REPLACE "$" "" new_include ${new_tmp_include}) + list(APPEND all_include_dirs "${new_include}") # Update list in the current scope only + endforeach() + get_all_include_dirs(${depend_lib} "${all_include_dirs}" "${all_visited_libs}") + endif() + endforeach() + set(all_include_dirs ${all_include_dirs} PARENT_SCOPE) + set(all_visited_libs ${all_visited_libs} PARENT_SCOPE) + endif() +endfunction() + +# Python files to install +FILE(GLOB PyROLPyFiles ${CMAKE_CURRENT_SOURCE_DIR}/python/*.py) +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/pyrol) + +GET_PROPERTY( PTR_IMPL GLOBAL PROPERTY PTR_IMPL ) +GET_PROPERTY( PTR_DIR GLOBAL PROPERTY PTR_DIR ) +GET_PROPERTY( STACKTRACE_DIR GLOBAL PROPERTY STACKTRACE_DIR ) +GET_PROPERTY( PARAMETERLIST_DIR GLOBAL PROPERTY PARAMETERLIST_DIR ) +GET_PROPERTY( LA_DIR GLOBAL PROPERTY LA_DIR ) +GET_PROPERTY( LAPACK_DIR GLOBAL PROPERTY LAPACK_DIR ) +GET_PROPERTY( BLAS_DIR GLOBAL PROPERTY BLAS_DIR ) + +set(ROL_all_include_dirs "") +list(APPEND ROL_all_include_dirs "${CMAKE_CURRENT_BINARY_DIR}/../src") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeB") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeE") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeG") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeU") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeB/pqn") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeG/augmentedlagrangian/") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeG/fletcher") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeG/interiorpoint") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeG/moreauyosida") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeG/stabilizedlcl") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeU/bundle") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeU/linesearch") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeU/linesearch/descent") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/algorithm/TypeU/trustregion") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/elementwise") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/dynamic") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/boundconstraint") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/constraint") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/nlls") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/polyproj") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/objective") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/operator") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/simopt") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/sketching") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/function/std") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/distribution") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/expectationquad") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/risk") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/risk/fdivergence") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/risk/spectral") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/error") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/regret") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/deviation") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/randvarfunctional/probability") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/function/progressivehedging") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/algorithm") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/algorithm/PrimalDual") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/vector") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/sampler") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/sampler/SROM") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/sol/status") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/status") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/augmentedlagrangian") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/bundle") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/fletcher") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/interiorpoint") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/krylov") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/linesearch") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/moreauyosidapenalty") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/nonlinearcg") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/secant") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/step/trustregion") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/utils") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/utils/function_bindings") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/vector") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/zoo") +list(APPEND ROL_all_include_dirs "${PTR_DIR}") +list(APPEND ROL_all_include_dirs "${STACKTRACE_DIR}") +list(APPEND ROL_all_include_dirs "${PARAMETERLIST_DIR}") +list(APPEND ROL_all_include_dirs "${LA_DIR}") +list(APPEND ROL_all_include_dirs "${LAPACK_DIR}") +list(APPEND ROL_all_include_dirs "${BLAS_DIR}") + +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed/constraint") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed/factors") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed/objective") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed/objective/het") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed/objective/hom") +list(APPEND ROL_all_include_dirs "${${PACKAGE_NAME}_SOURCE_DIR}/src/oed/utilities") + +IF (PYROL_GENERATE_SRC) + MESSAGE("-- PyROL_BINDER_EXECUTABLE:") + IF(NOT DEFINED PyROL_BINDER_EXECUTABLE) + find_program(PyROL_BINDER_EXECUTABLE + NAMES binder + ) + MESSAGE(" -- CMake has set: PyROL_BINDER_EXECUTABLE = ${PyROL_BINDER_EXECUTABLE}") + ELSE() + MESSAGE(" -- User has set: PyROL_BINDER_EXECUTABLE = ${PyROL_BINDER_EXECUTABLE}") + ENDIF() + file(REMOVE_RECURSE ${CMAKE_CURRENT_BINARY_DIR}/include_tmp) + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include_tmp) + file(REMOVE_RECURSE ${CMAKE_CURRENT_BINARY_DIR}/binder) + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/binder) + + MESSAGE("PTR_IMPL = ${PTR_IMPL}") + IF( PTR_IMPL STREQUAL "std::shared_ptr") + SET( BINDER_CFG "${CMAKE_CURRENT_SOURCE_DIR}/scripts/PyROL_shared_ptr.cfg" ) + ELSE() + SET( BINDER_CFG "${CMAKE_CURRENT_SOURCE_DIR}/scripts/PyROL_RCP.cfg" ) + ENDIF() + + set(binder_include_name "${CMAKE_CURRENT_BINARY_DIR}/trilinos_includes.hpp") + set(all_header_with_dir_list "${CMAKE_CURRENT_BINARY_DIR}/list_with_dir.txt") + set(all_header_without_dir_list "${CMAKE_CURRENT_BINARY_DIR}/list_without_dir.txt") + set(all_ETI_classes_list "${CMAKE_CURRENT_BINARY_DIR}/all_ETI_classes_list.txt") + set(all_include_list "${CMAKE_CURRENT_BINARY_DIR}/list_include.txt") + + set(all_include_dirs "") + set(all_visited_libs "") + foreach(depPkg IN LISTS ROL_LIB_ENABLED_DEPENDENCIES) + get_all_include_dirs(${depPkg}::all_libs "${all_include_dirs}" "${all_visited_libs}") + endforeach() + + list(REMOVE_DUPLICATES all_include_dirs) + list(REMOVE_ITEM all_include_dirs "") + + FOREACH(include_dir IN LISTS ROL_all_include_dirs) + list(APPEND all_include_dirs "${include_dir}") + ENDFOREACH(include_dir) + #MESSAGE("all_include_dirs = ${all_include_dirs}") + + set(PyROL_all_include_files_with_dir "") + set(PyROL_all_include_files_without_dir "") + foreach(include_dir IN LISTS all_include_dirs) + file(GLOB include_files + "${include_dir}/*.hpp" + "${include_dir}/*.h" + "${include_dir}/*/*.hpp" + "${include_dir}/*/*.h" + "${include_dir}/*/*/*.hpp" + "${include_dir}/*/*/*.h" + "${include_dir}/*/*/*/*.hpp" + "${include_dir}/*/*/*/*.h" + "${include_dir}/*/*/*/*.inc" + "${include_dir}/*/*/*/*.inc_predicate" + ) + foreach(include_file IN LISTS include_files) + list(APPEND PyROL_all_include_files_with_dir "${include_file}") + string(REPLACE "${include_dir}/" "" include_file_without_dir "${include_file}") + list(APPEND PyROL_all_include_files_without_dir "${include_file_without_dir}") + endforeach() + endforeach() + + + #list(REMOVE_DUPLICATES PyROL_all_include_files_without_dir) + #list(REMOVE_ITEM PyROL_all_include_files_without_dir "") + + #list(REMOVE_DUPLICATES PyROL_all_include_files_with_dir) + #list(REMOVE_ITEM PyROL_all_include_files_with_dir "") + + #MESSAGE("PyROL_all_include_files_with_dir = ${PyROL_all_include_files_with_dir}") + + SET(CONTENTS "") + FOREACH(line IN LISTS all_include_dirs) + SET(CONTENTS "${CONTENTS}${line}\n") + ENDFOREACH(line) + file(WRITE ${all_include_list} ${CONTENTS}) + + SET(CONTENTS "") + FOREACH(line IN LISTS PyROL_all_include_files_with_dir) + SET(CONTENTS "${CONTENTS}${line}\n") + ENDFOREACH(line) + file(WRITE ${all_header_with_dir_list} ${CONTENTS}) + + SET(CONTENTS "") + FOREACH(line IN LISTS PyROL_all_include_files_without_dir) + SET(CONTENTS "${CONTENTS}${line}\n") + ENDFOREACH(line) + file(WRITE ${all_header_without_dir_list} ${CONTENTS}) + + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/python) + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/src) + + file (GLOB PyROLPyFiles2 "${CMAKE_CURRENT_BINARY_DIR}/python/*.py") + list (APPEND PyROLPyFiles ${PyROLPyFiles2}) + + EXECUTE_PROCESS(COMMAND + ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/scripts/gather_includes.py ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${all_header_with_dir_list} ${all_header_without_dir_list} ${binder_include_name} + ) + + set(BINDER_OPTIONS "") + list(APPEND BINDER_OPTIONS ${binder_include_name}) + list(APPEND BINDER_OPTIONS --root-module pyrol) + list(APPEND BINDER_OPTIONS --prefix ${CMAKE_CURRENT_BINARY_DIR}/binder) + list(APPEND BINDER_OPTIONS -max-file-size=1000000) + list(APPEND BINDER_OPTIONS --bind Teuchos) + list(APPEND BINDER_OPTIONS --bind ROL) + list(APPEND BINDER_OPTIONS --bind pyrol) + IF(PYROL_B_VERBOSE) + list(APPEND BINDER_OPTIONS -v) + ENDIF() + list(APPEND BINDER_OPTIONS --config ${BINDER_CFG}) + IF(PYROL_SUPPRESS_ERRORS) + list(APPEND BINDER_OPTIONS --suppress-errors) + ENDIF() + list(APPEND BINDER_OPTIONS --) + IF(TPL_ENABLE_CUDA) + list(APPEND BINDER_OPTIONS -x cuda --cuda-host-only) + ENDIF() + list(APPEND BINDER_OPTIONS ${PyROL_BINDER_FLAGS}) + list(APPEND BINDER_OPTIONS -std=c++${CMAKE_CXX_STANDARD}) + list(APPEND BINDER_OPTIONS -I${CMAKE_CURRENT_BINARY_DIR}/include_tmp) + list(APPEND BINDER_OPTIONS -I${CMAKE_CURRENT_BINARY_DIR}/src) + list(APPEND BINDER_OPTIONS -I${CMAKE_CURRENT_SOURCE_DIR}/src) + + list(APPEND BINDER_OPTIONS -I${PyROL_BINDER_clang_include_dirs}) + list(APPEND BINDER_OPTIONS -iwithsysroot${PyROL_BINDER_LibClang_include_dir}) + IF(PyROL_BINDER_GCC_TOOLCHAIN) + list(APPEND BINDER_OPTIONS --gcc-toolchain=${PyROL_BINDER_GCC_TOOLCHAIN}) + ENDIF() + list(APPEND BINDER_OPTIONS -DNDEBUG) + + message("BINDER_OPTIONS='${BINDER_OPTIONS}'") + + EXECUTE_PROCESS(COMMAND + ${PyROL_BINDER_EXECUTABLE} ${BINDER_OPTIONS} + RESULT_VARIABLE STATUS + OUTPUT_VARIABLE OUTPUT_BINDER + ) + + if(STATUS AND NOT STATUS EQUAL 0) + message("${OUTPUT_BINDER}") + if(PYROL_CMAKE_ERROR) + message(FATAL_ERROR "BINDER FAILED: ${STATUS}") + else() + message("BINDER FAILED: ${STATUS}") + endif() + else() + message(STATUS "BINDER SUCCESS:") + message("${OUTPUT_BINDER}") + endif() + + IF(PYROL_UPDATE_GENERATED_SRC) + FILE(REMOVE_RECURSE ${CMAKE_CURRENT_SOURCE_DIR}/binder) + FILE(COPY ${CMAKE_CURRENT_BINARY_DIR}/binder DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}) + ENDIF() +ELSE() + # Copy binder files to ${CMAKE_CURRENT_BINARY_DIR} + FILE(COPY ${CMAKE_CURRENT_SOURCE_DIR}/binder DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +ENDIF() + +# Get the python version +EXECUTE_PROCESS(COMMAND ${PYTHON_EXECUTABLE} -c + "import sys; print(sys.version_info.major)" + OUTPUT_VARIABLE PYTHON_MAJOR_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE +) +EXECUTE_PROCESS(COMMAND ${PYTHON_EXECUTABLE} -c + "import sys; print(sys.version_info.minor)" + OUTPUT_VARIABLE PYTHON_MINOR_VERSION + OUTPUT_STRIP_TRAILING_WHITESPACE +) + +SET(PYTHON_VERSION ${PYTHON_MAJOR_VERSION}.${PYTHON_MINOR_VERSION}) + +SET(PYBIND11_PYTHON_VERSION ${PYTHON_VERSION}) + +# Determine the install directory +IF(PYROL_SCIKIT) + SET(PyROL_INSTALL_DIR + "." + CACHE STRING "The path where PyROL will be installed." + ) +ELSE() + SET(PyROL_INSTALL_DIR + lib/python${PYTHON_VERSION}/site-packages/pyrol + CACHE STRING "The path where PyROL will be installed" + ) +ENDIF() + +cmake_path(IS_RELATIVE PyROL_INSTALL_DIR is_pyrol_install_relative) + +IF(NOT is_pyrol_install_relative) + MESSAGE(FATAL_ERROR "PyROL install dir is not relative.") +ENDIF() + +cmake_path(GET PyROL_INSTALL_DIR RELATIVE_PART PyROL_INSTALL_DIR_RELATIVE) + +MESSAGE(STATUS "PyROL installation path: ${PyROL_INSTALL_DIR_RELATIVE}") + +INSTALL(FILES + ${PyROLPyFiles} + DESTINATION ${PyROL_INSTALL_DIR}) + +MACRO(SUBDIRLIST result curdir) + FILE(GLOB children RELATIVE ${curdir} ${curdir}/*) + SET(dirlist "") + FOREACH(child ${children}) + IF(IS_DIRECTORY ${curdir}/${child}) + LIST(APPEND dirlist ${curdir}/${child}) + ENDIF() + ENDFOREACH() + SET(${result} ${dirlist}) +ENDMACRO() + +SUBDIRLIST(PyROL_SUBDIRS ${CMAKE_CURRENT_SOURCE_DIR}/python) + +FOREACH(PyROL_subdir ${PyROL_SUBDIRS}) + INSTALL(DIRECTORY + ${PyROL_subdir} + DESTINATION ${PyROL_INSTALL_DIR}) +ENDFOREACH() + +# Find the pybind11 CMake module +EXECUTE_PROCESS(COMMAND + ${PYTHON_EXECUTABLE} -c "import pybind11; print(pybind11.get_cmake_dir())" + OUTPUT_VARIABLE pybind11_DIR + ERROR_VARIABLE pybind11_CMAKE_ERROR + OUTPUT_STRIP_TRAILING_WHITESPACE + ) +MESSAGE(STATUS "pybind11 CMake path: ${pybind11_DIR}") + +find_package(pybind11 REQUIRED) + +EXECUTE_PROCESS(COMMAND + ${PYTHON_EXECUTABLE} -c "import mpi4py; print(mpi4py.get_include())" + OUTPUT_VARIABLE Mpi4Py_INCLUDE_DIR + ERROR_VARIABLE Mpi4Py_INCLUDE_ERROR + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + +ADD_SUBDIRECTORY( src ) + +file(COPY ${PyROLPyFiles} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/pyrol/.) + +FOREACH(PyROL_subdir ${PyROL_SUBDIRS}) + FILE(COPY ${PyROL_subdir} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/pyrol) +ENDFOREACH() + +SET(PyROL_PYTHONPATH "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}") + +TRIBITS_ADD_TEST_DIRECTORIES(test) diff --git a/packages/rol/pyrol/example/.gitignore b/packages/rol/pyrol/example/.gitignore new file mode 100644 index 000000000000..541cb64f9b85 --- /dev/null +++ b/packages/rol/pyrol/example/.gitignore @@ -0,0 +1 @@ +test.txt \ No newline at end of file diff --git a/packages/rol/pyrol/example/dynamic_interface/forward_euler_1d.py b/packages/rol/pyrol/example/dynamic_interface/forward_euler_1d.py new file mode 100644 index 000000000000..d4205df724e2 --- /dev/null +++ b/packages/rol/pyrol/example/dynamic_interface/forward_euler_1d.py @@ -0,0 +1,221 @@ +# TO-DO: put in separate file(s) ############################################## + +from abc import ABC, abstractmethod + + +class VectorField(ABC): + + @abstractmethod + def value(self, u, z, t): + pass + + @abstractmethod + def applyJacobian_u(self, jv, v, u, z, t): + pass + + @abstractmethod + def applyJacobian_z(self, jv, v, u, z, t): + pass + + +class DecayingExponential(VectorField): + + def value(self, u, z, t): + return z[:] - u[:] + + def applyJacobian_u(self, jv, v, u, z, t): + jv[:] = -v[:] + + def applyJacobian_z(self, jv, v, u, z, t): + jv[:] = +v[:] + + +# TO-DO: put in separate file ################################################# + +from pyrol.unsupported import DynamicConstraint + + +class ForwardEuler(DynamicConstraint): + + def __init__(self, f): + super().__init__() + self.f = f + + def value(self, c, uo, un, z, ts): + dt = ts.t[1] - ts.t[0] + c[:] = un[:] - uo[:] - dt*self.f.value(uo, z, ts.t[0]) + + def solve(self, c, uo, un, z, ts): + dt = ts.t[1] - ts.t[0] + un[:] = uo[:] + dt*self.f.value(uo, z, ts.t[0]) + self.value(c, uo, un, z, ts) + + def applyJacobian_uo(self, jv, vo, uo, un, z, ts): + self.f.applyJacobian_u(jv, vo, uo, z, ts.t[0]) + dt = ts.t[1] - ts.t[0] + jv[:] = - vo[:] - dt*jv[:] + + def applyJacobian_un(self, jv, vn, uo, un, z, ts): + jv[:] = vn[:] + + def applyJacobian_z(self, jv, vz, uo, un, z, ts): + self.f.applyJacobian_z(jv, vz, uo, z, ts.t[0]) + dt = ts.t[1] - ts.t[0] + jv[:] = - dt*jv[:] + + def applyAdjointJacobian_uo(self, ajv, vo, uo, un, z, ts): + self.applyJacobian_uo(ajv, vo, uo, un, z, ts) + + def applyAdjointJacobian_un(self, ajv, vn, uo, un, z, ts): + self.applyJacobian_un(ajv, vn, uo, un, z, ts) + + def applyAdjointJacobian_z(self, ajv, vz, uo, un, z, ts): + self.applyJacobian_z(ajv, vz, uo, un, z, ts) + + def applyInverseJacobian_un(self, ijv, vn, uo, un, z, ts): + self.applyJacobian_un(ijv, vn, uo, un, z, ts) + + def applyInverseAdjointJacobian_un(self, iajv, vn, uo, un, z, ts): + self.applyJacobian_un(iajv, vn, uo, un, z, ts) + + # Gauss-Newton for now + + +# TO-DO: put in separate file ################################################# + +from pyrol.unsupported import DynamicObjective + + +class SquaredLoss(DynamicObjective): + + def __init__(self, T, y): + super().__init__() + self.T = T + self.y = y + + def value(self, uo, un, z, ts): + v = 0 + if ts.t[1] == self.T: + v = 0.5*np.sum((un[:] - self.y)**2) + return v + + def gradient_uo(self, g, uo, un, z, ts): + g[:] = 0 + + def gradient_un(self, g, uo, un, z, ts): + g[:] = 0 + if ts.t[1] == self.T: + g[:] = un[:] - self.y + + def gradient_z(self, g, uo, un, z, ts): + g[:] = 0 + + +# TO-DO: put in separate file ################################################# + +import numpy as np +import matplotlib.pyplot as plt + +from pyrol.vectors import NumPyVector +from pyrol import getCout, Problem, Solver + +from pyrol.unsupported import TimeStamp +from pyrol.unsupported import ReducedDynamicObjective +from pyrol.unsupported import PartitionedVector +from pyrol.unsupported import SerialConstraint + +from pyrol.pyrol.Teuchos import ParameterList +import pyrol + + +def as_numpy(v): + # v : ROL:PartitionedVector + n = v.dimension() + w = np.full(n, np.nan) + for k in range(n): + w[k] = v[k][:][0] + return w + + +def get_trajectory(constraint, u0, z, timestamps): + serial_constraint = SerialConstraint(constraint, u0, timestamps) + n = z.dimension() + c = PartitionedVector.create(NumPyVector(np.array([np.nan])), n) + u = PartitionedVector.create(NumPyVector(np.array([np.nan])), n) + serial_constraint.solve(c, u, z, np.nan) + t = np.full(n, np.nan) + for k in range(n): + t[k] = timestamps[k].t[0] + u = np.hstack([u0[:][0], as_numpy(u)[:-1]]) + return t, u + + +def plot_trajectory(constraint, u0, z, timestamps): + t, u = get_trajectory(constraint, u0, z, timestamps) + plt.plot(t, np.exp(-t), label='continuous') + plt.plot(t, u, '.-', label='numerical solution') + plt.grid() + plt.title(f'# of timesteps = {z.dimension()}') + plt.ylabel('state') + plt.xlabel('t') + plt.legend() + plt.show() + + +def optimize(n): + + T = 1 # end time + + # configure timestamps + dt = T/n + timestamps = pyrol.pyrol.std.vector_ROL_TimeStamp_double_t() + for k in range(n): + ts = TimeStamp() + ts.t.push_back(k*dt) + ts.t.push_back((k + 1)*dt) + timestamps.push_back(ts) + + # configure DynamicConstraint + f = DecayingExponential() + y = np.exp(-1) # the output we want to match (an everywhere zero control) + constraint = ForwardEuler(f) + + # configure DynamicObjective + objective = SquaredLoss(T, y) + + parameters = ParameterList() + cout = getCout() + + # configure ReducedDynamicObjective + u0 = NumPyVector(np.ones(1)) + zk = NumPyVector(np.zeros(1)) + ck = NumPyVector(np.zeros(1)) + reduced_objective = ReducedDynamicObjective(objective, constraint, u0, zk, ck, timestamps, parameters, cout) + + # Problem + z = PartitionedVector.create(zk, n) + problem = Problem(reduced_objective, z) + # problem.check(True, cout) + + parameters['General'] = ParameterList() + parameters['General']['Output Level'] = 1 + parameters['Step'] = ParameterList() + parameters['Step']['Trust Region'] = ParameterList() + parameters['Step']['Trust Region']['Subproblem Solver'] = 'Dogleg' + + # Solver + solver = Solver(problem, parameters) + solver.solve(cout) + + return constraint, u0, z, timestamps + + +def main(): + n = 100 # number of time steps + constraint, u0, z, timestamps = optimize(n) + print(f'\n||z||_inf = {np.linalg.norm(as_numpy(z), np.inf)}\n') + plot_trajectory(constraint, u0, z, timestamps) + + +if __name__ == '__main__': + main() diff --git a/packages/rol/pyrol/example/old/regression/definitions.py b/packages/rol/pyrol/example/old/regression/definitions.py new file mode 100644 index 000000000000..783a4a7c07c8 --- /dev/null +++ b/packages/rol/pyrol/example/old/regression/definitions.py @@ -0,0 +1,74 @@ +import numpy as np + +from pyrol import * +from pyrol.vectors import npVector + + +class myVector(npVector): + + def norm(self): + return np.sqrt(self.dot(self)) + + def __sub__(self, other): + return myVector(self.values - other.values) + + +class ForwardProblem(Constraint_SimOpt): + + def __init__(self, X): + self.X = np.reshape(X, (X.size, 1)) + self.X = np.hstack([self.X, np.ones_like(self.X)]) + super().__init__() + + def value(self, c, u, z, tol): + c.values[:] = u.values - self.X@z.values + + def applyJacobian_1(self, jv, v, u, z, tol): + jv.values[:] = v.values + + def applyJacobian_2(self, jv, v, u, z, tol): + jv.values[:] = - self.X@v.values + + def applyInverseJacobian_1(self, ijv, v, u, z, tol): + ijv.values[:] = v.values + + def applyAdjointJacobian_1(self, ajv, v, u, z, tol): + ajv.values[:] = v.values + + def applyAdjointJacobian_2(self, ajv, v, u, z, tol): + ajv.values[:] = - self.X.T@v.values + + def applyInverseAdjointJacobian_1(self, iajv, v, u, z, tol): + iajv.values[:] = v.values + + def solve(self, c, u, z, tol): + u.values[:] = self.X@z.values + self.value(c, u, z, tol) + + +class Loss(Objective_SimOpt): + + def __init__(self, Y): + self.Y = Y + super().__init__() + + def value(self, u, z, tol): + return (self.Y - u).norm()**2/2 + + def gradient_1(self, g, u, z, tol): + g.values[:] = u.values - self.Y.values + + def gradient_2(self, g, u, z, tol): + g.zero() + + def hessVec_11(self, hv, v, u, z, tol): + hv.values[:] = v.values + + def hessVec_12(self, hv, v, u, z, tol): + hv.zero() + + def hessVec_21(self, hv, v, u, z, tol): + hv.zero() + + def hessVec_22(self, hv, v, u, z, tol): + hv.zero() diff --git a/packages/rol/pyrol/example/old/regression/input.xml b/packages/rol/pyrol/example/old/regression/input.xml new file mode 100644 index 000000000000..8823aab2dffc --- /dev/null +++ b/packages/rol/pyrol/example/old/regression/input.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + \ No newline at end of file diff --git a/packages/rol/pyrol/example/old/regression/testPolymorphism.py b/packages/rol/pyrol/example/old/regression/testPolymorphism.py new file mode 100644 index 000000000000..e0739b0f035f --- /dev/null +++ b/packages/rol/pyrol/example/old/regression/testPolymorphism.py @@ -0,0 +1,50 @@ +import numpy as np +from numpy.random import default_rng + +from pyrol import * +from definitions import myVector, ForwardProblem + +class Loss(Objective_SimOpt): + + def __init__(self, Y): + self.Y = Y + super().__init__() + + def value_3(self, *args): + u, z, tol = args + return (self.Y - u).norm()**2/2 + +def main(): + + a, b = 1, 0 # the "true" parameter values from which we generate data + n = 100 # the size of our data set + + # Create predictors. ####################### + rng = default_rng(23040901) + X = rng.uniform(-1, 1, n) + + # Create the forward problem. ############## + constraint = ForwardProblem(X) + + # Set up vectors. ########################## + u = myVector(np.empty(n)) + z = myVector(np.empty(2)) + + # Create responses. ######################## + Y = u.clone() + constraint.solve(u.clone(), Y, myVector(np.array([a, b])), np.nan) + noise = 0.1*rng.normal(0, 1, n) + Y.values += noise + + # Set up the objective. #################### + objective = Loss(Y) + + # Evaluate the objective. ################## + x = Vector_SimOpt(u, z) + tol = 0.0 + print(objective.value(u, z, tol)) + print(objective.value(x, tol)) + + +if __name__ == "__main__": + main() diff --git a/packages/rol/pyrol/example/old/regression/testRegression.py b/packages/rol/pyrol/example/old/regression/testRegression.py new file mode 100644 index 000000000000..f880e0aa90ae --- /dev/null +++ b/packages/rol/pyrol/example/old/regression/testRegression.py @@ -0,0 +1,55 @@ +import numpy as np +from numpy.random import default_rng + +from pyrol import * +from definitions import myVector, ForwardProblem, Loss + + +def main(): + + a, b = 1, 0 # the "true" parameter values from which we generate data + n = 100 # the size of our data set + + # Configure parameter list. ################ + params = getParametersFromXmlFile("input.xml") + + # Set the output stream. ################### + stream = getCout() + + # Create predictors. ####################### + rng = default_rng(23040901) + X = rng.uniform(-1, 1, n) + + # Create the forward problem. ############## + constraint = ForwardProblem(X) + + # Set up vectors. ########################## + u = myVector(np.empty(n)) + z = myVector(np.empty(2)) + g = z.dual() # gradient with respect to z + l = u.dual() # adjoint variables + + # Create responses. ######################## + Y = u.clone() + constraint.solve(u.clone(), Y, myVector(np.array([a, b])), np.nan) + noise = 0.1*rng.normal(0, 1, n) + Y.values += noise + + # Set up the problem. ###################### + objective = Loss(Y) + reducedObjective = Reduced_Objective_SimOpt(objective, constraint, u, z, l) + problem = Problem(reducedObjective, z, g) + problem.check(True, stream) + + # Solve. ################################### + solver = Solver(problem, params) + solver.solve(stream) + + # Check the solution. ###################### + target, _, _, _ = np.linalg.lstsq(constraint.X, objective.Y.values, rcond=None) + np.testing.assert_allclose(z.values, target) + print('Test passed: Correctly optimized!') + + +if __name__ == "__main__": + main() diff --git a/packages/rol/pyrol/example/old/test.py b/packages/rol/pyrol/example/old/test.py new file mode 100644 index 000000000000..fe96196fa11d --- /dev/null +++ b/packages/rol/pyrol/example/old/test.py @@ -0,0 +1,94 @@ +from PyROL.vectors.npVector import npVector as vector_type +from PyROL import * +from PyROL.PyROL.ROL import getParametersFromXmlFile +import numpy as np + + +class norm2Obj(Objective): + def __init__(self, H, g, c=0): + self.H = H + self.g = g + self.c = c + super().__init__() + + def value(self, x, tol): + tmp = x.clone() + x.zero() + self.H.apply(tmp, x, tol) + tmp.scale(0.5) + tmp.plus(self.g) + return x.apply(tmp) + self.c + + def gradient(self, g, x, tol): + self.H.apply(g, x, tol) + g.plus(self.g) + + def hessVec(self, hv, v, x, tol): + self.H.apply(hv, v, tol) + + def invHessVec(self, hv, v, x, tol): + self.H.applyInverse(hv, v, tol) + + +# Matrix from rol/example/quadratic/example_01.cpp +class matrix(LinearOperator): + def __init__(self, dim): + self.dim = dim + super().__init__() + + def apply(self, Hv, v, tol): + for i in range(0, self.dim): + Hv[i] = 2.*v[i] + if i > 0: + Hv[i] -= v[i-1] + if i < self.dim - 1: + Hv[i] -= v[i+1] + + +op = matrix(10) +g = vector_type.full(10, 0.) +x = vector_type.full(10, 1.) +obj = norm2Obj(op, g) +c = Constraint() +a = vector_type.full(10, 1.) +b = vector_type.full(10, 1.) +a.scale(2.) +np.testing.assert_equal(a.norm(), 2*np.sqrt(10)) +a.zero() +np.testing.assert_equal(a.norm(), 0.) + +a.axpy(1., b) +np.testing.assert_equal(a.norm(), np.sqrt(10)) +np.testing.assert_equal(a.dot(b), 10.) + +b.setScalar(2.) +a = b +np.testing.assert_equal(a.apply(b), 40.) +np.testing.assert_equal(b.norm(), 2*np.sqrt(10)) +np.testing.assert_equal(b[0], 2.) +b[0:2] = [1., 3.] +np.testing.assert_equal(b[0:2], [1., 3.]) + +print(obj.value(a, 1e-8)) + +op.apply(b,a,1e-8) +print(b[0:3]) +g = vector_type.full(10, 0.) + +params = getParametersFromXmlFile("input.xml") + +problem = Problem(obj, x) +#solver = Solver(problem, params) + +#print(params) + +obj_q = QuadraticObjective(op, g) +step = TrustRegionStep(params) +status = StatusTest(params) +algo = Algorithm(step,status,False) + +print("before optimization: "+str(obj_q.value(x, 1e-8))) +algo.run_void(x=x, obj=obj_q) +print("after optimization: "+str(obj_q.value(x, 1e-8))) +print(g[:]) +print(x[:]) diff --git a/packages/rol/pyrol/example/old/test_bounds.py b/packages/rol/pyrol/example/old/test_bounds.py new file mode 100644 index 000000000000..9ee3734cbd25 --- /dev/null +++ b/packages/rol/pyrol/example/old/test_bounds.py @@ -0,0 +1,46 @@ +from PyROL.vectors.npVector import npVector as vector_type +from PyROL import * +from PyROL.PyROL.ROL import getParametersFromXmlFile +from PyROL.PyROL.ROL import Bounds_double_t, Problem_double_t, Solver_double_t +from PyROL.PyROL.ROL import getCout +import sys +import numpy as np + +# Matrix from rol/example/quadratic/example_01.cpp +class matrix(LinearOperator): + def __init__(self, dim): + self.dim = dim + super().__init__() + + def apply(self, Hv, v, tol): + for i in range(0, self.dim): + Hv[i] = 2.*v[i] + if i > 0: + Hv[i] -= v[i-1] + if i < self.dim - 1: + Hv[i] -= v[i+1] + +op = matrix(10) +g = vector_type.full(10, 0.) +x = vector_type.full(10, 1.) + +# passes: +# bnd = Bounds_double_t(g.clone(), x.clone()) +# fails: +bnd = Bounds_double_t(vector_type.full(10, 0.), x.clone()) +obj = QuadraticObjective(op, g) + +params = getParametersFromXmlFile("input.xml") +status = StatusTest(params) + +problem = Problem_double_t(obj, x) +problem.addBoundConstraint(bnd) +solver = Solver_double_t(problem, params) +cout = openOfstream('test.txt') +solver.solve(cout, status) +closeOfstream(cout) + +state = solver.getAlgorithmState() +print(state.iter) +print(state.value) +print(state.nfval) diff --git a/packages/rol/pyrol/example/old/test_reduce.py b/packages/rol/pyrol/example/old/test_reduce.py new file mode 100644 index 000000000000..b5760b25fc41 --- /dev/null +++ b/packages/rol/pyrol/example/old/test_reduce.py @@ -0,0 +1,9 @@ +from PyROL.PyROL import ROL + +op = ROL.Elementwise.ReductionMin_double_t() + +v = 10. +w = 20. +minVal = min(v, w) +op.reduce(v, w) +assert abs(w-minVal) < 1e-7, (v, w, minVal) diff --git a/packages/rol/pyrol/example/old/typeB_test_02.py b/packages/rol/pyrol/example/old/typeB_test_02.py new file mode 100644 index 000000000000..0a6b041dc2a8 --- /dev/null +++ b/packages/rol/pyrol/example/old/typeB_test_02.py @@ -0,0 +1,9 @@ +from PyROL.PyROL import ROL + +paramlist = ROL.ParameterList +paramlist.sublist("Status Test").set("Gradient Tolerance",1e-8) +paramlist.sublist("Status Test").set("Constraint Tolerance",1e-8) +paramlist.sublist("Status Test").set("Step Tolerance",1e-12) +paramlist.sublist("Status Test").set("Iteration Limit", 250) +paramlist.sublist("Step").set("Type","Line Search") +paramlist.sublist("General").set("Output Level",iprint) \ No newline at end of file diff --git a/packages/rol/pyrol/example/pytorch/Models.py b/packages/rol/pyrol/example/pytorch/Models.py new file mode 100644 index 000000000000..e69680764058 --- /dev/null +++ b/packages/rol/pyrol/example/pytorch/Models.py @@ -0,0 +1,23 @@ +import torch + + +class LinearLayer(torch.nn.Module): + '''A linear transformation.''' + + def __init__(self, input_size): + super().__init__() + self.linear = torch.nn.Linear(input_size, 1) + self.input_size = input_size + + def forward(self, x): + return self.linear(x) + +class NN(torch.nn.Module): + + def __init__(self, input_size): + super().__init__() + self.linear = torch.nn.Linear(input_size, 1) + self.input_size = input_size + + def forward(self, x): + return torch.tanh(self.linear(x)) diff --git a/packages/rol/pyrol/example/pytorch/TorchObjectives.py b/packages/rol/pyrol/example/pytorch/TorchObjectives.py new file mode 100644 index 000000000000..8d3c86dd6526 --- /dev/null +++ b/packages/rol/pyrol/example/pytorch/TorchObjectives.py @@ -0,0 +1,54 @@ +from pyrol import Objective + +import torch + + +class TorchObjective(Objective): + # https://pytorch.org/docs/stable/func.html + + # @staticmethod + # def _copy(source, target): + # target.zero() + # target.plus(source) + + def __init__(self): + super().__init__() + self.torch_gradient = torch.func.grad(self.torch_value) + + def value(self, x, tol): + return self.torch_value(x.torch_object).item() + + def gradient(self, g, x, tol): + ans = self.torch_gradient(x.torch_object) + g.torch_object = ans + + def _forward_over_reverse(self, input, x, v): + # https://github.com/google/jax/blob/main/docs/notebooks/autodiff_cookbook.ipynb + return torch.func.jvp(input, (x,), (v,)) + + def hessVec(self, hv, v, x, tol): + input = torch.func.grad(self.torch_value) + _, ans = self._forward_over_reverse(input, x.torch_object, v.torch_object) + hv.torch_object = ans + + def torch_value(self, x): + # Returns a scalar torch Tensor + raise NotImplementedError + + +class SquaredNorm(TorchObjective): + + def torch_value(self, x): + return 0.5*x.squeeze()**2 + + +class LeastSquaresObjective(TorchObjective): + + def __init__(self, data, model): + super().__init__() + self.x, self.y = data + self.model = model + self.loss = torch.nn.MSELoss(reduction='sum') + + def torch_value(self, x): + return 0.5*self.loss(torch.func.functional_call(self.model, x, self.x), self.y) diff --git a/packages/rol/pyrol/example/pytorch/TorchVectors.py b/packages/rol/pyrol/example/pytorch/TorchVectors.py new file mode 100644 index 000000000000..6c58dd8ea6fe --- /dev/null +++ b/packages/rol/pyrol/example/pytorch/TorchVectors.py @@ -0,0 +1,258 @@ +from pyrol.pyrol import ROL +from pyrol.getTypeName import * + +import torch + +import copy + + +class PythonVector(getTypeName('Vector')): + + def __iadd__(self, other): + self.plus(other) + + # TO-DO: Add other operations. + + +class TensorVector(PythonVector): + + @torch.no_grad() + def __init__(self, tensor): + super().__init__() + assert isinstance(tensor, torch.Tensor) + self.torch_object = tensor + + @property + def tensor(self): + return self.torch_object + + @torch.no_grad() + def axpy(self, alpha, other): + self.tensor.add_(other.tensor, alpha=alpha) + + @torch.no_grad() + def scale(self, alpha): + self.tensor.mul_(alpha) + + @torch.no_grad() + def zero(self): + self.tensor.zero_() + + @torch.no_grad() + def dot(self, other): + ans = torch.sum(torch.mul(self.tensor, other.tensor)) + return ans.item() + + @torch.no_grad() + def clone(self): + tensor = copy.deepcopy(self.tensor) + ans = TensorVector(tensor) + ans.zero() + return ans + + @torch.no_grad() + def dimension(self): + return self.tensor.numel() + + @torch.no_grad() + def setScalar(self, alpha): + self.fill_(alpha) + + @torch.no_grad() + def __getitem__(self, index): + flat = self.tensor.view(-1) + return flat[index].item() + + @torch.no_grad() + def __setitem__(self, index, value): + flat = self.tensor.view(-1) + flat[index] = value + + # Derived methods ######################################################### + + @torch.no_grad() + def reduce(self, op): + reduction_type = op.reductionType() + match reduction_type: + case ROL.Elementwise.REDUCE_MIN: + ans = self.tensor.min() + ans = ans.item() + case ROL.Elementwise.REDUCE_MAX: + ans = self.tensor.max() + ans = ans.item() + case ROL.Elementwise.REDUCE_SUM: + ans = torch.sum(self.tensor) + ans = ans.item() + case ROL.Elementwise.REDUCE_AND: + ans = self.tensor.all() + ans = ans.item() + case ROL.Elementwise.REDUCE_BOR: + ans = 0 + for i in range(self.dimension()): + ans = ans | int(self[i].item()) + case _: + raise NotImplementedError(reduction_type) + return ans + + @torch.no_grad() + def plus(self, other): + self.axpy(1, other) + + @torch.no_grad() + def norm(self): + return self.dot(self)**0.5 + + @torch.no_grad() + def basis(self, i): + b = self.clone() + b.zero() + b[i] = 1 + return b + + #### + + @torch.no_grad() + def applyUnary(self, op): + for i in range(self.dimension()): + self[i] = op.apply(self[i]) + + @torch.no_grad() + def applyBinary(self, other, op): + for i in range(self.dimension()): + self[i] = op.apply(self[i], other[i]) + + +class TensorDictVector(PythonVector): + + @torch.no_grad() + def __init__(self, tensor_dict): + super().__init__() + assert isinstance(tensor_dict, dict) + self.torch_object = tensor_dict + + @property + def tensor_dict(self): + return self.torch_object + + @torch.no_grad() + def axpy(self, alpha, other): + for k, v in self.tensor_dict.items(): + v.add_(other.tensor_dict[k], alpha=alpha) + + @torch.no_grad() + def scale(self, alpha): + for _, v in self.tensor_dict.items(): + v.mul_(alpha) + + @torch.no_grad() + def zero(self): + for _, v in self.tensor_dict.items(): + v.zero_() + + @torch.no_grad() + def dot(self, other): + ans = 0 + for k, v in self.tensor_dict.items(): + ans += torch.sum(torch.mul(v, other.tensor_dict[k])) + return ans.item() + + @torch.no_grad() + def clone(self): + tensor_dict = copy.deepcopy(self.tensor_dict) + ans = TensorDictVector(tensor_dict) + ans.zero() + return ans + + @torch.no_grad() + def dimension(self): + # TO-DO: Cache value + ans = 0 + for _, v in self.tensor_dict.items(): + ans += v.numel() + return ans + + @torch.no_grad() + def setScalar(self, alpha): + for _, v in self.tensor_dict.items(): + v.fill_(alpha) + + @torch.no_grad() + def __getitem__(self, index): + total = 0 + for _, v in self.tensor_dict.items(): + numel = v.numel() + if index < numel + total: + flat = v.view(-1) + return flat[index - total].item() + total += numel + + @torch.no_grad() + def __setitem__(self, index, value): + total = 0 + for _, v in self.tensor_dict.items(): + numel = v.numel() + if index < numel + total: + flat = v.view(-1) + flat[index - total] = value + return + total += numel + + # Derived methods ######################################################### + + @torch.no_grad() + def reduce(self, op): + reduction_type = op.reductionType() + match reduction_type: + case ROL.Elementwise.REDUCE_MIN: + ans = float('+inf') + for _, v in self.torch_dict.items(): + ans = min(ans, v.min().item()) + case ROL.Elementwise.REDUCE_MAX: + ans = float('-inf') + for _, v in self.torch_dict.items(): + ans = max(ans, v.max().item()) + case ROL.Elementwise.REDUCE_SUM: + ans = 0 + for _, v in self.torch_dict.items(): + ans += torch.sum(v) + ans = ans.item() + case ROL.Elementwise.REDUCE_AND: + ans = True + for _, v in self.torch_dict.items(): + ans = ans and v.all().item() + if ans == False: + break + case ROL.Elementwise.REDUCE_BOR: + ans = 0 + for i in range(self.dimension()): + ans = ans | int(self[i].item()) + case _: + raise NotImplementedError(reduction_type) + return ans + + @torch.no_grad() + def plus(self, other): + self.axpy(1, other) + + @torch.no_grad() + def norm(self): + return self.dot(self)**0.5 + + @torch.no_grad() + def basis(self, i): + b = self.clone() + b.zero() + b[i] = 1 + return b + + #### + + @torch.no_grad() + def applyUnary(self, op): + for i in range(self.dimension()): + self[i] = op.apply(self[i]) + + @torch.no_grad() + def applyBinary(self, other, op): + for i in range(self.dimension()): + self[i] = op.apply(self[i], other[i]) diff --git a/packages/rol/pyrol/example/pytorch/nn_example.py b/packages/rol/pyrol/example/pytorch/nn_example.py new file mode 100644 index 000000000000..ccc9c6e987f7 --- /dev/null +++ b/packages/rol/pyrol/example/pytorch/nn_example.py @@ -0,0 +1,74 @@ +from TorchVectors import TensorDictVector +from TorchObjectives import LeastSquaresObjective +from Models import NN + +from pyrol import getCout, Problem, Solver +from pyrol.pyrol.Teuchos import ParameterList + +import numpy as np + +import torch + + +def generate_data(model): + n = 10 # number of data points + torch.manual_seed(24010101) + x = torch.normal(0, 1, size=(n, model.input_size)) + x.requires_grad = False + # TO-DO: Consider adding noise to y. + model.eval() + with torch.no_grad(): + y = model.forward(x) + return x, y + + +def build_parameter_list(): + params = ParameterList() + params['General'] = ParameterList() + params['General']['Output Level'] = 1 + params['Step'] = ParameterList() + params['Step']['Trust Region'] = ParameterList() + params['Step']['Trust Region']['Subproblem Solver'] = 'Truncated CG' + params['Step']['Trust Region']['Initial Radius'] = 1e2 + return params + + +def check_result(x, alpha): + target = [] + actual = [] + for i in range(x.dimension()): + target.append(alpha) + actual.append(x[i]) + np.testing.assert_allclose(actual, target) + print('Check passed! Correctly optimized.') + + +def main(): + torch.set_default_dtype(torch.float64) + + input_size = 2 + model = NN(input_size) + m = TensorDictVector(model.state_dict()) # model parameters + alpha = 1.1 + m.setScalar(alpha) + data = generate_data(model) + + x = m.clone() # our optimization variable + + objective = LeastSquaresObjective(data, model) + g = x.clone() + + stream = getCout() + + problem = Problem(objective, x, g) + problem.checkDerivatives(True, stream) + + params = build_parameter_list() + solver = Solver(problem, params) + solver.solve(stream) + + check_result(x, alpha) + + +if __name__ == '__main__': + main() diff --git a/packages/rol/pyrol/example/pytorch/tensor_example.py b/packages/rol/pyrol/example/pytorch/tensor_example.py new file mode 100644 index 000000000000..537e68b686de --- /dev/null +++ b/packages/rol/pyrol/example/pytorch/tensor_example.py @@ -0,0 +1,43 @@ +from TorchVectors import TensorVector +from TorchObjectives import SquaredNorm + +from pyrol import getCout, Problem, Solver +from pyrol.pyrol.Teuchos import ParameterList + +import torch + + +def build_parameter_list(): + params = ParameterList() + params['General'] = ParameterList() + params['General']['Output Level'] = 1 + params['Step'] = ParameterList() + params['Step']['Trust Region'] = ParameterList() + params['Step']['Trust Region']['Subproblem Solver'] = 'Truncated CG' + params['Step']['Trust Region']['Initial Radius'] = 1e2 + return params + + +def main(): + torch.set_default_dtype(torch.float64) + + x = torch.tensor([[3.]],requires_grad=False) + x = TensorVector(x) + + objective = SquaredNorm() + g = x.clone() + + stream = getCout() + + problem = Problem(objective, x, g) + problem.checkDerivatives(True, stream) + + params = build_parameter_list() + solver = Solver(problem, params) + solver.solve(stream) + + print(f'x = {x.torch_object}') + + +if __name__ == '__main__': + main() diff --git a/packages/rol/pyrol/example/rosenbrock.py b/packages/rol/pyrol/example/rosenbrock.py new file mode 100644 index 000000000000..3b1bd50ea1ae --- /dev/null +++ b/packages/rol/pyrol/example/rosenbrock.py @@ -0,0 +1,70 @@ +from pyrol import getCout, Objective, Problem, Solver +from pyrol.vectors import NumPyVector + +from pyrol.pyrol.Teuchos import ParameterList + +import numpy as np + + +class RosenbrockObjective(Objective): + + def __init__(self): + self.alpha = 100 + super().__init__() + + def value(self, x, tol): + v = self.alpha*(x[0]**2 - x[1])**2 + (x[0] - 1)**2 + return v + + def gradient(self, g, x, tol): + g[0] = 2*self.alpha*(x[0]**2 - x[1])*2*x[0] + g[0] += 2*(x.array[0] - 1) + g[1] = -2*self.alpha*(x[0]**2 - x[1]) + + def hessVec(self, hv, v, x, tol): + h11 = 12*self.alpha*x[0]**2 - 4*self.alpha*x[1] + 2 + h12 = -4*self.alpha*x[0] + h22 = 2*self.alpha + hv[0] = h11*v[0] + h12*v[1] + hv[1] = h12*v[0] + h22*v[1] + + +def build_parameter_list(): + params = ParameterList() + params['General'] = ParameterList() + params['General']['Output Level'] = 1 + params['Step'] = ParameterList() + params['Step']['Trust Region'] = ParameterList() + params['Step']['Trust Region']['Subproblem Solver'] = 'Truncated CG' + return params + + +def main(): + + # Configure parameter list. ################ + params = build_parameter_list() + + # Set the output stream. ################### + stream = getCout() + + # Set up vectors. ########################## + x = NumPyVector(np.array([-3., -4.])) + g = x.dual() + + # Set up the problem. ###################### + objective = RosenbrockObjective() + problem = Problem(objective, x, g) + problem.check(True, stream) + + # Solve. ################################### + solver = Solver(problem, params) + solver.solve(stream) + + # Check the solution. ###################### + target = np.array([1., 1.]) + np.testing.assert_allclose(x.array, target) + print('Test passed: Correctly optimized!') + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/packages/rol/pyrol/pyproject.toml b/packages/rol/pyrol/pyproject.toml new file mode 100644 index 000000000000..eadbd789546b --- /dev/null +++ b/packages/rol/pyrol/pyproject.toml @@ -0,0 +1,49 @@ +[build-system] +requires = ["scikit-build-core>=0.3.3","pybind11 @ git+https://github.com/pybind/pybind11.git@smart_holder"] +build-backend = "scikit_build_core.build" + + +[project] +name = "pyrol" +version = "0.0.1" +description="A Python interface to the Rapid Optimization Library (ROL)" +readme = "README.md" +authors = [ + { name = "Christian Glusa" }, + { name = "Kim Liegeois" }, + { name = "Aurya Javeed" }, +] +requires-python = ">=3.7" +dependencies = ["numpy"] +classifiers = [ + "Development Status :: 4 - Beta", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: 3 :: Only", + "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", +] + +[tool.scikit-build] +wheel.expand-macos-universal-tags = true +cmake.verbose = true +logging.level = "DEBUG" +cmake.minimum-version = "3.23.0" +wheel.install-dir = "pyrol" + +[tool.scikit-build.cmake.define] +CMAKE_BUILD_TYPE = "RELEASE" +Trilinos_ASSERT_DEFINED_DEPENDENCIES = "OFF" +BUILD_SHARED_LIBS = "ON" +Trilinos_ENABLE_TESTS = "OFF" +Trilinos_ENABLE_EXAMPLES = "OFF" +Trilinos_ENABLE_Fortran = "OFF" +Trilinos_ENABLE_ALL_OPTIONAL_PACKAGES = "OFF" +Trilinos_ENABLE_ROL = "ON" +ROL_Ptr = "Teuchos::RCP" +ROL_ENABLE_PYROL = "ON" +PYROL_ENABLE_BINDER = "OFF" +PYROL_PIP_INSTALL = "ON" +CMAKE_INSTALL_RPATH="$ORIGIN/lib64;$ORIGIN/lib;$ORIGIN;@loader_path/lib64;@loader_path/lib;@loader_path" diff --git a/packages/rol/pyrol/python/__init__.py b/packages/rol/pyrol/python/__init__.py new file mode 100644 index 000000000000..0aa9682944c4 --- /dev/null +++ b/packages/rol/pyrol/python/__init__.py @@ -0,0 +1,45 @@ +import importlib + +from . getTypeName import getTypeName, getDefaultScalarType, ROL_classes, ROL_members + +__version__ = '0.1.0' +def version(): + return 'PyROL version: ' + __version__ + + +def getWrapper(classname): + def wrapper(scalarType=getDefaultScalarType()): + actualClass = getTypeName(classname, scalarType) + return actualClass + return wrapper + + +defaultScalarType = getDefaultScalarType() + +supported_objects = {"Bounds", "Constraint", "LinearOperator", + "LinearConstraint", "Objective", "Problem", "Solver", + "Vector", "getCout", "getParametersFromXmlFile"} + +unsupported = importlib.import_module('.unsupported', 'pyrol') + +for classnameLong in ROL_members: + class_obj, _ = ROL_members[classnameLong] + pos = classnameLong.find('_'+defaultScalarType+'_t') + if pos <= 0: + if classnameLong in supported_objects: + locals().update({classnameLong: class_obj}) + else: + setattr(unsupported, classnameLong, class_obj) + continue + + classname = classnameLong[:pos] + + if classname in supported_objects: + locals().update({classname: getTypeName(classname, defaultScalarType)}) + locals().update({classname+'_forScalar': getWrapper(classname)}) + else: + setattr(unsupported, classname, getTypeName(classname, defaultScalarType)) + setattr(unsupported, classname + '_forScalar', getWrapper(classname)) + +del importlib, getTypeName, getDefaultScalarType, ROL_classes, ROL_members +del class_obj, classname, classnameLong, getWrapper, pos diff --git a/packages/rol/pyrol/python/getTypeName.py b/packages/rol/pyrol/python/getTypeName.py new file mode 100644 index 000000000000..7838425e1f0a --- /dev/null +++ b/packages/rol/pyrol/python/getTypeName.py @@ -0,0 +1,59 @@ +from pyrol.pyrol import ROL +import inspect +import sys + +trackedTypes = [] + +def track(self, *args, **kwargs): + for arg in args: + if isinstance(arg, tuple(trackedTypes)): + self._tracked_constructor_args.append(arg) + for key in kwargs: + val = kwargs[key] + if isinstance(val, tuple(trackedTypes)): + self._tracked_constructor_args.append(val) + +def tracking_method(cls_method): + def new_method(self, *args, **kwargs): + track(self, *args, **kwargs) + return cls_method(self, *args, **kwargs) + return new_method + +def withTrackingConstructor(cls): + class newCls(cls): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + self._tracked_constructor_args = [] + + track(self, *args, **kwargs) + method_names = [m for m in dir(cls) if callable(getattr(cls, m))] + + for name in method_names: + if name.startswith('add'): + setattr(newCls, name, tracking_method(getattr(cls, name))) + + return newCls + + +ROL_members = {} +for cls_name, cls_obj in inspect.getmembers(sys.modules['pyrol.pyrol.ROL']): + if inspect.isclass(cls_obj): + cls_obj = withTrackingConstructor(cls_obj) + trackedTypes.append(cls_obj) + setattr(sys.modules['pyrol.pyrol.ROL'], cls_name, cls_obj) + ROL_members[cls_name] = (cls_obj, inspect.isclass(cls_obj)) + + +ROL_classes = [(cls_name , cls_obj) for cls_name, cls_obj in inspect.getmembers(sys.modules['pyrol.pyrol.ROL']) if inspect.isclass(cls_obj)] + +def getDefaultScalarType(): + return "double" + +def getTypeName(class_name, scalar_type=getDefaultScalarType()): + class_name_scalar = class_name.lower()+"_"+scalar_type.lower() + for i in range(0, len(ROL_classes)): + if ROL_classes[i][0].lower().find(class_name_scalar) == 0: + return ROL_classes[i][1] + print("Warning: Unknown type \"{}\", the function returns None.".format(class_name)) + return None diff --git a/packages/rol/pyrol/python/unsupported.py b/packages/rol/pyrol/python/unsupported.py new file mode 100644 index 000000000000..fd1415d630cc --- /dev/null +++ b/packages/rol/pyrol/python/unsupported.py @@ -0,0 +1,2 @@ +# A sub-module that holds the unsupported classes. +# Content gets populated automatically. diff --git a/packages/rol/pyrol/python/vectors/NumPyVector.py b/packages/rol/pyrol/python/vectors/NumPyVector.py new file mode 100644 index 000000000000..64b1273067ba --- /dev/null +++ b/packages/rol/pyrol/python/vectors/NumPyVector.py @@ -0,0 +1,86 @@ +import numpy as np +from numpy import linalg as LA +from pyrol.pyrol import ROL +from pyrol.getTypeName import * + + +class NumPyVector(getTypeName('Vector')): + def __init__(self, array=None): + assert isinstance(array, np.ndarray) + assert array.ndim == 1 + self.array = array + super().__init__() + + @staticmethod + def full(dimension=1, default_value=0.): + array = np.full((dimension,), fill_value=default_value) + return NumPyVector(array) + + def plus(self, b): + self.array += b.array + + def scale(self, scale_factor): + self.array *= scale_factor + + def dot(self, b): + return np.vdot(self.array, b.array) + + def norm(self): + return np.sqrt(self.dot(self)) + + def zero(self): + self.setScalar(0.) + + def clone(self): + tmp = type(self)(np.full(self.array.shape, fill_value=0.)) + return tmp + + def axpy(self, scale_factor, x): + ax = x.clone() + ax.zero() + ax.plus(x) + ax.scale(scale_factor) + self.plus(ax) + + def dimension(self): + return len(self.array) + + def setScalar(self, new_value): + self.array[:] = new_value + + def reduce(self, op): + reductionType = op.reductionType() + if reductionType == ROL.Elementwise.REDUCE_MIN: + return self.array.min() + elif reductionType == ROL.Elementwise.REDUCE_MAX: + return self.array.max() + elif reductionType == ROL.Elementwise.REDUCE_SUM: + return self.array.sum() + elif reductionType == ROL.Elementwise.REDUCE_AND: + return np.logical_and.reduce(self.array) + elif reductionType == ROL.Elementwise.REDUCE_BOR: + return np.bitwise_or.reduce(self.array) + else: + raise NotImplementedError(reductionType) + + def applyUnary(self, op): + for i in range(self.dimension()): + self.array[i] = op.apply(self.array[i]) + + def applyBinary(self, op, other): + assert self.dimension() == other.dimension() + for i in range(self.dimension()): + self.array[i] = op.apply(self.array[i], other.array[i]) + + def __getitem__(self, index): + return self.array[index] + + def __setitem__(self, index, val): + self.array[index] = val + + def basis(self, i): + b = self.clone() + b.array[:] = 0. + b.array[i] = 1. + self._basis = b + return b diff --git a/packages/rol/pyrol/python/vectors/__init__.py b/packages/rol/pyrol/python/vectors/__init__.py new file mode 100644 index 000000000000..76f8813a315c --- /dev/null +++ b/packages/rol/pyrol/python/vectors/__init__.py @@ -0,0 +1,4 @@ +from . NumPyVector import NumPyVector +# from . tVector import tVector + +__all__ = [NumPyVector] diff --git a/packages/rol/pyrol/python/vectors/tVector.py b/packages/rol/pyrol/python/vectors/tVector.py new file mode 100644 index 000000000000..2e330063e127 --- /dev/null +++ b/packages/rol/pyrol/python/vectors/tVector.py @@ -0,0 +1,83 @@ +from mpi4py import MPI +from PyROL.PyROL import ROL +from PyTrilinos2.PyTrilinos2 import Teuchos +from PyTrilinos2.getTpetraTypeName import * + +class tVector(ROL.Vector_double_t): + def __init__(self, dimension=1, default_value=0., map=None, comm=None, scalar_type=getDefaultScalarType(), local_ordinal_type=getDefaultLocalOrdinalType(), global_ordinal_type=getDefaultGlobalOrdinalType(), node_type=getDefaultNodeType()): + self.scalar_type = scalar_type + self.local_ordinal_type = local_ordinal_type + self.global_ordinal_type = global_ordinal_type + self.node_type = node_type + self.mapType = getTypeName('Map', local_ordinal_type=self.local_ordinal_type, global_ordinal_type=self.global_ordinal_type, node_type=self.node_type) + self.vectorType = getTypeName('Vector', scalar_type=self.scalar_type, local_ordinal_type=self.local_ordinal_type, global_ordinal_type=self.global_ordinal_type, node_type=self.node_type) + if map is None: + if comm is None: + comm = Teuchos.getTeuchosComm(MPI.COMM_WORLD) + map = self.mapType(dimension, 0, comm) + self.tvector = self.vectorType(map, False) + self.tvector.putScalar(default_value) + super().__init__() + + def plus(self, b): + self.tvector.update(1., b.tvector, 1.) + + def scale(self, scale_factor): + self.tvector.scale(scale_factor) + + def dot(self, b): + return self.tvector.dot(b.tvector) + + def norm(self): + return self.tvector.norm2() + + def clone(self): + tmp = tVector(map=self.tvector.getMap(), scalar_type=self.scalar_type, local_ordinal_type=self.local_ordinal_type, global_ordinal_type=self.global_ordinal_type, node_type=self.node_type) + return tmp + + def axpy(self, scale_factor, x): + self.tvector.update(scale_factor, x.tvector, 1.) + + def dimension(self): + return self.tvector.getMap().getGlobalNumElements() + + def setScalar(self, new_value): + self.tvector.putScalar(new_value) + + def __getitem__(self, index): + if isinstance( index, int ): + map = self.tvector.getMap() + if map.isNodeGlobalElement(index): + local_index = map.getLocalElement(index) + view = self.tvector.getLocalViewHost() + return view[local_index] + if isinstance( index, slice ): + map = self.tvector.getMap() + view = self.tvector.getLocalViewHost() + global_indices = range(*index.indices(self.dimension())) + local_indices = np.empty(np.size(global_indices), dtype=int) + for i in range(0, len(global_indices)): + if map.isNodeGlobalElement(global_indices[i]): + local_indices[i] = map.getLocalElement(global_indices[i]) + else: + local_indices[i] = 0 + return view[local_indices] + + def __setitem__(self, index, val): + if isinstance( index, int ): + map = self.tvector.getMap() + if map.isNodeGlobalElement(index): + self.tvector.replaceGlobalValue(index, val) + if isinstance( index, slice ): + map = self.tvector.getMap() + global_indices = range(*index.indices(self.dimension())) + for i in range(0, len(global_indices)): + if map.isNodeGlobalElement(global_indices[i]): + if len(val) > 1: + self.tvector.replaceGlobalValue(global_indices[i], val[i]) + else: + self.tvector.replaceGlobalValue(global_indices[i], val) + + def reduce(self, op): + reductionType = op.reductionType() + raise NotImplementedError(reductionType) diff --git a/packages/rol/pyrol/scripts/PyROL_RCP.cfg b/packages/rol/pyrol/scripts/PyROL_RCP.cfg new file mode 100644 index 000000000000..f587597d046c --- /dev/null +++ b/packages/rol/pyrol/scripts/PyROL_RCP.cfg @@ -0,0 +1,161 @@ +#+include +#+include + + +################################################## +# Teuchos # +################################################## + +# Tell binder about Teuchos RCP ++include ++custom_shared Teuchos::RCP + +-class Teuchos::CommandLineProcessor +-class Teuchos::Ptr >> +-class Teuchos::ArrayView >> +-class Teuchos::RCPNodeTmpl +-class Teuchos::DeallocDelete +-class Teuchos::MpiComm +-class Teuchos::is_comparable +-class Teuchos::ArrayRCP +-class Teuchos::LAPACK +-class Teuchos::BLAS +-class Teuchos::is_printable +-class Teuchos::compare +-class Teuchos::print +-namespace Teuchos::Details +-function Teuchos::TimeMonitor::computeGlobalTimerStatistics +-function Teuchos::mpiErrorCodeToString +-class Teuchos::CommandLineProcessor::enum_opt_data_t +-class Teuchos::CommandLineProcessor::TimeMonitorSurrogate +-class Teuchos::RawWorkspace +-class Teuchos::Describable +-class Teuchos::MpiCommRequestBase +-function Teuchos::getRawMpiComm +-class Teuchos::OpaqueWrapper +-class Teuchos::OpaqueWrapper +# -function Teuchos::Details::setMpiReductionOp +# -function Teuchos::Details::getMpiOpForEReductionType +-class Teuchos::OpaqueWrapper +-function Teuchos::rcp_dynamic_cast +-class Teuchos::VerboseObjectBase +-function Teuchos::ParameterList::sublist +-function Teuchos::ParameterList::set +-function Teuchos::ParameterList::get ++include_for_namespace Teuchos ++add_on_binder Teuchos::ParameterList def_ParameterList_member_functions + +-class Teuchos::FilteredIterator + +-class Teuchos::TypeNameTraits +-function Teuchos::typeName +-class Teuchos::GetBaseObjVoidPtrImpl +-class Teuchos::RCPNodeHandle +-class Teuchos::ArrayView +-function Teuchos::getArrayFromStringParameter + +-function Teuchos::typeName +-function Teuchos::rcp +-function Teuchos::is_null +-class Teuchos::EmbeddedObjDealloc + +################################################## +# ROL # +################################################## + +-function ROL::Secant::get_state +-function ROL::lBFGS::get_state +-function ROL::lDFP::get_state +-function ROL::lSR1::get_state +-function ROL::BarzilaiBorwein::get_state +-function ROL::NonlinearCG::get_state +-function ROL::makePtr +-function ROL::SecantFactory +-function ROL::DescentDirectionUFactory +-function ROL::TypeB::AlgorithmFactory +-function ROL::TypeE::AlgorithmFactory +-function ROL::TypeG::AlgorithmFactory +-function ROL::TypeU::AlgorithmFactory +-class ROL::ConstraintAssembler + ++smart_holder ROL::Constraint_SimOpt ++smart_holder ROL::Objective_SimOpt ++smart_holder ROL::Reduced_Objective_SimOpt ++smart_holder ROL::DynamicFunction ++smart_holder ROL::DynamicConstraint ++smart_holder ROL::DynamicObjective ++smart_holder ROL::ReducedDynamicObjective ++smart_holder ROL::SerialFunction ++smart_holder ROL::SerialConstraint ++smart_holder ROL::Constraint ++smart_holder ROL::Objective ++smart_holder ROL::Vector_SimOpt ++smart_holder ROL::PartitionedVector ++smart_holder ROL::Vector + +# not required +-function ROL::Objective_SimOpt::value_2 + +# special handling of ROL::Vector::clone ++include_for_class ROL::Vector ++trampoline_member_function_binder ROL::Vector::clone customClone + ++include_for_namespace ROL::PyROL +-namespace ROL::details + +#################################################r +# std library # +################################################## + ++module_local_namespace @all_namespaces + +-class std::complex +#-class std::ostream +-class std::basic_ios +-class std::map +-class std::set +-class std::basic_streambuf + +-class std::filebuf +-class std::stringbuf +-class std::exception +-class std::invalid_argument +-class std::logic_error +-class std::invalid_arugment +-class std::runtime_error +-class std::type_info +-class std::bad_cast + +-class std::basic_ofstream +-class std::basic_istream +-class std::basic_istringstream +-class std::basic_ostringstream +#-class std::basic_ostream +-class std::fpos +-function std::ostream::_M_write +-function std::ostream::seekp +-function std::basic_ostream::_M_write +-function std::basic_ostream::seekp + +-function std::basic_ostream::operator<< + +-function std::streambuf::stossc +-function std::streambuf::__safe_gbump +-function std::streambuf::__safe_pbump +-function std::stringbuf::stossc +-function std::stringbuf::__safe_gbump +-function std::stringbuf::__safe_pbump +-function std::filebuf::stossc +-function std::filebuf::__safe_gbump +-function std::filebuf::__safe_pbump + +-class std::iterator +-class std::reverse_iterator + +-class std::vector +-class std::vector ++class std::vector ++class std::vector ++class std::vector> + +-namespace __gnu_cxx diff --git a/packages/rol/pyrol/scripts/PyROL_shared_ptr.cfg b/packages/rol/pyrol/scripts/PyROL_shared_ptr.cfg new file mode 100644 index 000000000000..6ab13c297cba --- /dev/null +++ b/packages/rol/pyrol/scripts/PyROL_shared_ptr.cfg @@ -0,0 +1,77 @@ ++include ++include + + +################################################## +# Teuchos # +################################################## + +-class Teuchos::CommandLineProcessor +-class Teuchos::Ptr >> +-class Teuchos::ArrayView >> +-class Teuchos::RCPNodeTmpl +-class Teuchos::DeallocDelete +-class Teuchos::MpiComm +-class Teuchos::is_comparable +-class Teuchos::ArrayRCP +-class Teuchos::LAPACK +-class Teuchos::BLAS +-class Teuchos::is_printable +-class Teuchos::compare +-class Teuchos::print +-namespace Teuchos::Details +-function Teuchos::TimeMonitor::computeGlobalTimerStatistics +-function Teuchos::mpiErrorCodeToString +-class Teuchos::CommandLineProcessor::enum_opt_data_t +-class Teuchos::CommandLineProcessor::TimeMonitorSurrogate +-class Teuchos::RawWorkspace +-class Teuchos::Describable +-class Teuchos::MpiCommRequestBase +-function Teuchos::getRawMpiComm +-class Teuchos::OpaqueWrapper +-class Teuchos::OpaqueWrapper +# -function Teuchos::Details::setMpiReductionOp +# -function Teuchos::Details::getMpiOpForEReductionType +-class Teuchos::OpaqueWrapper +-function Teuchos::rcp_dynamic_cast +-class Teuchos::VerboseObjectBase +-function Teuchos::ParameterList::sublist +-function Teuchos::ParameterList::set +-function Teuchos::ParameterList::get ++include_for_namespace Teuchos ++add_on_binder Teuchos::ParameterList def_ParameterList_member_functions + +################################################## +# ROL # +################################################## + +-function ROL::Secant::get_state +-function ROL::lBFGS::get_state +-function ROL::lDFP::get_state +-function ROL::lSR1::get_state +-function ROL::BarzilaiBorwein::get_state +-function ROL::NonlinearCG::get_state +-function ROL::makePtr +-function ROL::SecantFactory +-class ROL::ConstraintAssembler + +# special handling of ROL::Vector::clone ++include_for_class ROL::Vector ++trampoline_member_function_binder ROL::Vector::clone customClone + ++include_for_namespace ROL::PyROL +-namespace ROL::details + + +################################################## +# std library # +################################################## + ++module_local_namespace @all_namespaces + +-class std::complex +-class std::ostream +-class std::basic_ios +-class std::vector +-class std::map +-class std::set diff --git a/packages/rol/pyrol/scripts/create_sdist b/packages/rol/pyrol/scripts/create_sdist new file mode 100755 index 000000000000..a19ba13ba4fc --- /dev/null +++ b/packages/rol/pyrol/scripts/create_sdist @@ -0,0 +1,76 @@ +##!/bin/sh + +# This script builds an SDist from the ROL repository. It's meant to be copied +# to and then run from the directory containing the ROL repository. + + +# - Users of this script should check that the variables below are defined +# properly. +############################################################################## +TRILINOS_VERSION="14.5" +REPO_NAME="ROL-Trilinos" + +LLVM_PREFIX=$(spack location -i llvm) +LLVM_VERSION=$(echo ${LLVM_PREFIX} | awk -F[\-\-] '{print $5}') +GCC_PREFIX=$(spack location -i gcc) +############################################################################## + +## Other prerequisites: + +# * Binder (after the changes on its smart_holder branch) +if [ ! command -v binder &> /dev/null ] +then + echo "Binder not available!" +fi + +# * pybind11 built from its smart_holder branch. +# * pipx + +## + +cp ./${REPO_NAME}/packages/rol/pyrol/pyproject.toml ./${REPO_NAME}/pyproject.toml + +## Step 1: Use CMake to create a build directory. +[ -d build ] && rm -rf build +cmake -G Ninja \ +-D CMAKE_BUILD_TYPE:STRING=RELEASE \ +-D Trilinos_ENABLE_CPACK_PACKAGING=ON \ +-D Trilinos_CPACK_SOURCE_GENERATOR="TGZ" \ +-D BUILD_SHARED_LIBS:BOOL=ON \ +-D Trilinos_ENABLE_TESTS:BOOL=OFF \ +-D Trilinos_ENABLE_EXAMPLES:BOOL=OFF \ +-D Trilinos_ENABLE_Fortran:BOOL=OFF \ +-D Trilinos_ENABLE_ALL_OPTIONAL_PACKAGES:BOOL=OFF \ +-D Trilinos_ENABLE_ROL:BOOL=ON \ +-D ROL_ENABLE_PYROL:BOOL=ON \ +-D ROL_ENABLE_TESTS:BOOL=OFF \ +-D ROL_ENABLE_EXAMPLES:BOOL=OFF \ +-D PYROL_ENABLE_BINDER:BOOL=ON \ +-D PYROL_ENABLE_BINDER_UPDATE:BOOL=ON \ +-D PyROL_BINDER_clang_include_dirs:PATH="${LLVM_PREFIX}/include" \ +-D PyROL_BINDER_LibClang_include_dir:PATH="${LLVM_PREFIX}/lib/clang/${LLVM_VERSION}/include" \ +-D PyROL_BINDER_GCC_TOOLCHAIN:PATH="${GCC_PREFIX}" \ +-D CPACK_SOURCE_IGNORE_FILES="/packages/rol/example;/packages/rol/test;/packages/rol/tutorial;/packages/teuchos/core/test;/packages/teuchos/core/example;/packages/teuchos/numerics/test;/packages/teuchos/numerics/example;/packages/teuchos/parameterlist/test;/packages/teuchos/parameterlist/example;/packages/teuchos/parser/test;/packages/teuchos/parser/example;/packages/teuchos/remainder/test;/packages/teuchos/remainder/example;/packages/teuchos/comm/test;/packages/teuchos/comm/example" \ +-D CMAKE_INSTALL_PREFIX:PATH=install \ +./${REPO_NAME} -B./build + +## Step 2: Run Binder. +make -C build + +## Step 3: Create the reduced tarball. +make package_source -C build +TARBALL_NAME="trilinos-${TRILINOS_VERSION}-Source" + +## Step 4: Unpack the reduced tarball. +[ -d ${TARBALL_NAME} ] && rm -rf ${TARBALL_NAME} +tar -zxf "build/${TARBALL_NAME}.tar.gz" +mv ${TARBALL_NAME} pyrol + +## Step 5: Create an SDist from the tarball. +python -m pipx run build --sdist pyrol +cp -r pyrol/dist/* . + +## Step 6: Clean up. +rm -rf build +rm -rf ${REPO_NAME}/packages/rol/pyrol/binder +rm -rf pyrol diff --git a/packages/rol/pyrol/scripts/gather_includes.py b/packages/rol/pyrol/scripts/gather_includes.py new file mode 100644 index 000000000000..c1e60c71742c --- /dev/null +++ b/packages/rol/pyrol/scripts/gather_includes.py @@ -0,0 +1,114 @@ +import glob +import os +import sys + +def get_without_subfolder(line): + last_index=-1 + for i in range(len(line)): + if line[i] == '/': + last_index = i + if line[i] == '<': + first_index = i + if last_index == -1: + return line + return line[:first_index+1]+line[last_index+1:] + +def get_angular_include(line, remove_subfolder=False): + first=True + i0 = 0 + i1 = 0 + newline = line + for i in range(len(newline)): + if newline[i] == '"': + if first: + newline = newline[:i] + '<' + newline[i+1:] + first = False + i0 = i+1 + else: + newline = newline[:i] + '>' + newline[i+1:] + i1 = i + if newline[i0:i1] in ['storage_class.h', 'cuda_cc7_asm_atomic_op.inc_predicate', 'cuda_cc7_asm_atomic_fetch_op.inc_predicate']: + return line + if remove_subfolder: + return get_without_subfolder(newline) + return newline + +def make_all_includes(all_include_filename, folders): + all_includes = [] + for folder in folders: + for filename in (glob.glob(f'{folder}/**/*.hpp', recursive=True) + + glob.glob(f'{folder}/**/*.cpp', recursive=True) + + glob.glob(f'{folder}/**/*.h', recursive=True) + + glob.glob(f'{folder}/**/*.cc', recursive=True) + + glob.glob(f'{folder}/**/*.c', recursive=True)): + with open(filename, 'r') as fh: + for line in fh: + if line.startswith('#include'): + all_includes.append(get_angular_include(line).strip()) + all_includes = list(set(all_includes)) + # This is to ensure that the list is always the same and doesn't + # depend on the filesystem state. Not technically necessary, but + # will cause inconsistent errors without it. + all_includes.sort() + with open(all_include_filename, 'w') as fh: + for include in all_includes: + fh.write(f'{include}\n') + return all_include_filename + +def make_all_includes_from_filenames(all_include_filename, filenames): + all_includes = [] + for filename in filenames: + with open(filename, 'r') as fh: + for line in fh: + if line.startswith('#include'): + all_includes.append(get_angular_include(line).strip()) + all_includes = list(set(all_includes)) + # This is to ensure that the list is always the same and doesn't + # depend on the filesystem state. Not technically necessary, but + # will cause inconsistent errors without it. + all_includes.sort() + with open(all_include_filename, 'w') as fh: + for include in all_includes: + fh.write(f'{include}\n') + return all_include_filename + +# https://github.com/RosettaCommons/binder/issues/212 + +def copy_and_angular_includes(filenames, filenames_witout_dir, to_dir): + # loops over the files, replace include " " by include < > and write them in the to_dir: + for i in range(len(filenames)): + filename = filenames[i] + filename_witout_dir = filenames_witout_dir[i] + path=os.path.dirname(filename_witout_dir) + if not os.path.exists(to_dir+'/'+path): + os.makedirs(to_dir+'/'+path) + try: + with open(filename, 'r') as from_f: + lines = from_f.readlines() + except UnicodeDecodeError: + with open(filename, 'r', encoding='iso-8859-1') as from_f: + lines = from_f.readlines() + except PermissionError: + continue + + with open(to_dir+'/'+filename_witout_dir, 'w') as to_f: + for line in lines: + if line.startswith('#include'): + line = get_angular_include(line, False) + to_f.write(f'{line}') + +if __name__ == '__main__': + CMAKE_CURRENT_SOURCE_DIR = sys.argv[1] + CMAKE_CURRENT_BINARY_DIR = sys.argv[2] + all_header_list_with_dir = sys.argv[3] + all_header_list_without_dir = sys.argv[4] + binder_include_name = sys.argv[5] + + with open(all_header_list_with_dir, 'r') as fh: + all_include_filenames_with_dir = fh.read().splitlines() + + with open(all_header_list_without_dir, 'r') as fh: + all_include_filenames_without_dir = fh.read().splitlines() + + copy_and_angular_includes(all_include_filenames_with_dir, all_include_filenames_without_dir, CMAKE_CURRENT_BINARY_DIR+'/include_tmp') + make_all_includes_from_filenames(binder_include_name, [CMAKE_CURRENT_SOURCE_DIR+'/src/Binder_Input.hpp']) diff --git a/packages/rol/pyrol/src/Binder_Input.hpp b/packages/rol/pyrol/src/Binder_Input.hpp new file mode 100644 index 000000000000..17c20f2dbd12 --- /dev/null +++ b/packages/rol/pyrol/src/Binder_Input.hpp @@ -0,0 +1,3 @@ +#include +#include +#include \ No newline at end of file diff --git a/packages/rol/pyrol/src/CMakeLists.txt b/packages/rol/pyrol/src/CMakeLists.txt new file mode 100644 index 000000000000..1ebc1e8ab97f --- /dev/null +++ b/packages/rol/pyrol/src/CMakeLists.txt @@ -0,0 +1,29 @@ +FILE(GLOB PYROL_SRC "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp") +MESSAGE("PYROL_SRC = ${PYROL_SRC}") +FILE(COPY ${PYROL_SRC} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + +MESSAGE("CMAKE_CURRENT_BINARY_DIR = ${CMAKE_CURRENT_BINARY_DIR}") + +file(STRINGS ${CMAKE_CURRENT_BINARY_DIR}/../binder/pyrol.sources BINDER_SRCS) +list(TRANSFORM BINDER_SRCS PREPEND "${CMAKE_CURRENT_BINARY_DIR}/../binder/") +list(APPEND PYROL_SRC ${BINDER_SRCS}) + +MESSAGE("PYROL_SRC with binder = ${PYROL_SRC}") + +pybind11_add_module(pyrol ${PYROL_SRC}) +target_include_directories(pyrol PUBLIC ${Mpi4Py_INCLUDE_DIR} ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}/../binder "$") +target_compile_features(pyrol PUBLIC cxx_std_14) + +foreach(depPkg IN LISTS ROL_LIB_ENABLED_DEPENDENCIES) + target_link_libraries(pyrol PUBLIC ${depPkg}::all_libs) +endforeach() +target_link_libraries(pyrol PUBLIC ${Trilinos_EXTRA_LINK_FLAGS}) +set_target_properties(pyrol PROPERTIES SUFFIX ".so") + +INSTALL(TARGETS pyrol + DESTINATION ${PyROL_INSTALL_DIR}) + +add_custom_command(TARGET pyrol POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_BINARY_DIR}/pyrol.so ${CMAKE_CURRENT_BINARY_DIR}/../pyrol/. + COMMENT "Copy ${PROJECT_BINARY_DIR}/src/PyROL.so" +) diff --git a/packages/rol/pyrol/src/PyROL_ETI.hpp b/packages/rol/pyrol/src/PyROL_ETI.hpp new file mode 100644 index 000000000000..5cb0efd345b9 --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_ETI.hpp @@ -0,0 +1,68 @@ +#ifndef PYROL_ETI +#define PYROL_ETI + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#define BINDER_ETI_ABSTRACT(CLASS_NAME) \ + template class CLASS_NAME; + +#define BINDER_ETI_WITH_FOO(CLASS_NAME) \ + template class CLASS_NAME; \ + template <> inline void PyROL::foo(CLASS_NAME a){} + +#define BINDER_ROL_VECTOR(SCALAR) \ + BINDER_ETI_ABSTRACT(Vector) \ + BINDER_ETI_ABSTRACT(Vector_SimOpt) + +#define BINDER_ROL_OBJECTIVE(SCALAR) \ + BINDER_ETI_ABSTRACT(Objective) \ + BINDER_ETI_ABSTRACT(Objective_SimOpt) \ + BINDER_ETI_ABSTRACT(Reduced_Objective_SimOpt) \ + BINDER_ETI_ABSTRACT(ReducedDynamicObjective) + +#define BINDER_ROL_CONSTRAINT(SCALAR) \ + BINDER_ETI_ABSTRACT(Constraint) \ + BINDER_ETI_ABSTRACT(SimConstraint) \ + BINDER_ETI_ABSTRACT(BoundConstraint_SimOpt) \ + BINDER_ETI_ABSTRACT(SerialConstraint) + +#define BINDER_ROL_SOLVER(SCALAR) \ + BINDER_ETI_ABSTRACT(Solver) + +#define BINDER_ROL_PROBLEM(SCALAR) \ + BINDER_ETI_ABSTRACT(Problem) + +#define BINDER_ROL_OED(SCALAR) \ + BINDER_ETI_ABSTRACT(Factory) + +namespace ROL { + + BINDER_ROL_VECTOR(double) + BINDER_ROL_OBJECTIVE(double) + BINDER_ROL_CONSTRAINT(double) + BINDER_ROL_SOLVER(double) + BINDER_ROL_PROBLEM(double) + +namespace OED { + BINDER_ROL_OED(double) +} + +} + +#endif // PYROL_ETI diff --git a/packages/rol/pyrol/src/PyROL_ETI_helper.hpp b/packages/rol/pyrol/src/PyROL_ETI_helper.hpp new file mode 100644 index 000000000000..596c83b96e25 --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_ETI_helper.hpp @@ -0,0 +1,13 @@ +#ifndef PYROL_ETI_HELPER +#define PYROL_ETI_HELPER + +namespace ROL { +namespace PyROL { + + template + inline void foo(T a){} + +} +} + +#endif // PYROL_ETI_HELPER diff --git a/packages/rol/pyrol/src/PyROL_Teuchos_Custom.cpp b/packages/rol/pyrol/src/PyROL_Teuchos_Custom.cpp new file mode 100644 index 000000000000..2252cd96586e --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_Teuchos_Custom.cpp @@ -0,0 +1,174 @@ +#include + +#if PY_VERSION_HEX >= 0x03000000 + +#define PyClass_Check(obj) PyObject_IsInstance(obj, (PyObject *)&PyType_Type) +#define PyInt_Check(x) PyLong_Check(x) +#define PyInt_AsLong(x) PyLong_AsLong(x) +#define PyInt_FromLong(x) PyLong_FromLong(x) +#define PyInt_FromSize_t(x) PyLong_FromSize_t(x) +#define PyString_Check(name) PyBytes_Check(name) +#define PyString_FromString(x) PyUnicode_FromString(x) +#define PyString_FromStringAndSize(x,s) PyUnicode_FromStringAndSize(x,s) +#define PyString_Format(fmt, args) PyUnicode_Format(fmt, args) +#define PyString_AsString(str) PyBytes_AsString(str) +#define PyString_Size(str) PyBytes_Size(str) +#define PyString_InternFromString(key) PyUnicode_InternFromString(key) +#define Py_TPFLAGS_HAVE_CLASS Py_TPFLAGS_BASETYPE +#define PyString_AS_STRING(x) PyUnicode_AS_STRING(x) +#define PyObject_Compare(x, y) (1-PyObject_RichCompareBool(x, y, Py_EQ)) +#define _PyLong_FromSsize_t(x) PyLong_FromSsize_t(x) +#define convertPyStringToChar(pyobj) PyBytes_AsString(PyUnicode_AsASCIIString(pyobj)) +#else +#define convertPyStringToChar(pyobj) PyString_AsString(pyobj) +#endif + +namespace py = pybind11; + +// Implementation based on: +// https://github.com/trilinos/Trilinos/tree/master/packages/PyTrilinos/src/PyTrilinos_Teuchos_Util.cpp +bool setPythonParameter(Teuchos::RCP plist, + const std::string & name, + py::object value) +{ + py::handle h = value; + + // Boolean values + if (PyBool_Check(value.ptr ())) + { + plist->set(name, h.cast()); + } + + // Integer values + else if (PyInt_Check(value.ptr ())) + { + plist->set(name, h.cast()); + } + + // Floating point values + else if (PyFloat_Check(value.ptr ())) + { + plist->set(name, h.cast()); + } + + // Unicode values + else if (PyUnicode_Check(value.ptr ())) + { + PyObject * pyBytes = PyUnicode_AsASCIIString(value.ptr ()); + if (!pyBytes) return false; + plist->set(name, std::string(PyBytes_AsString(pyBytes))); + Py_DECREF(pyBytes); + } + + // String values + else if (PyString_Check(value.ptr ())) + { + plist->set(name, h.cast()); + } + + // None object not allowed: this is a python type not usable by + // Trilinos solver packages, so we reserve it for the + // getPythonParameter() function to indicate that the requested + // parameter does not exist in the given Teuchos::ParameterList. + // For logic reasons, this check must come before the check for + // Teuchos::ParameterList + else if (value.ptr () == Py_None) + { + return false; + } + + // All other value types are unsupported + else + { + return false; + } + + // Successful type conversion + return true; +} // setPythonParameter + + +// Implementation based on: +// https://github.com/trilinos/Trilinos/tree/master/packages/PyTrilinos/src/PyTrilinos_Teuchos_Util.cpp +py::object getPythonParameter(Teuchos::RCP plist, + const std::string & name) +{ + // Get the parameter entry. I now deal with the Teuchos::ParameterEntry + // objects so that I can query the Teuchos::ParameterList without setting + // the "used" flag to true. + const Teuchos::ParameterEntry * entry = plist->getEntryPtr(name); + // Boolean parameter values + if (entry->isType< bool >()) + { + bool value = Teuchos::any_cast< bool >(entry->getAny(false)); + return py::cast(value); + } + // Integer parameter values + else if (entry->isType< int >()) + { + int value = Teuchos::any_cast< int >(entry->getAny(false)); + return py::cast(value); + } + // Double parameter values + else if (entry->isType< double >()) + { + double value = Teuchos::any_cast< double >(entry->getAny(false)); + return py::cast(value); + } + // String parameter values + else if (entry->isType< std::string >()) + { + std::string value = Teuchos::any_cast< std::string >(entry->getAny(false)); + return py::cast(value.c_str()); + } + // Char * parameter values + else if (entry->isType< char * >()) + { + char * value = Teuchos::any_cast< char * >(entry->getAny(false)); + return py::cast(value); + } + + else if (entry->isArray()) + { + try + { + Teuchos::Array< int > tArray = + Teuchos::any_cast< Teuchos::Array< int > >(entry->getAny(false)); + return copyTeuchosArrayToNumPy(tArray); + } + catch(Teuchos::bad_any_cast &e) + { + try + { + Teuchos::Array< long > tArray = + Teuchos::any_cast< Teuchos::Array< long > >(entry->getAny(false)); + return copyTeuchosArrayToNumPy(tArray); + } + catch(Teuchos::bad_any_cast &e) + { + try + { + Teuchos::Array< float > tArray = + Teuchos::any_cast< Teuchos::Array< float > >(entry->getAny(false)); + return copyTeuchosArrayToNumPy(tArray); + } + catch(Teuchos::bad_any_cast &e) + { + try + { + Teuchos::Array< double > tArray = + Teuchos::any_cast< Teuchos::Array< double > >(entry->getAny(false)); + return copyTeuchosArrayToNumPy(tArray); + } + catch(Teuchos::bad_any_cast &e) + { + return py::none(); + } + } + } + } + } + + // All other types are unsupported + return py::none(); +} // getPythonParameter diff --git a/packages/rol/pyrol/src/PyROL_Teuchos_Custom.hpp b/packages/rol/pyrol/src/PyROL_Teuchos_Custom.hpp new file mode 100644 index 000000000000..f4e96ba42404 --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_Teuchos_Custom.hpp @@ -0,0 +1,99 @@ +#ifndef PYROL_TEUCHOS_CUSTOM +#define PYROL_TEUCHOS_CUSTOM + +#include +#include +#include +//#include +#include +#include + +#if PY_VERSION_HEX >= 0x03000000 + +#define PyClass_Check(obj) PyObject_IsInstance(obj, (PyObject *)&PyType_Type) +#define PyInt_Check(x) PyLong_Check(x) +#define PyInt_AsLong(x) PyLong_AsLong(x) +#define PyInt_FromLong(x) PyLong_FromLong(x) +#define PyInt_FromSize_t(x) PyLong_FromSize_t(x) +#define PyString_Check(name) PyBytes_Check(name) +#define PyString_FromString(x) PyUnicode_FromString(x) +#define PyString_FromStringAndSize(x,s) PyUnicode_FromStringAndSize(x,s) +#define PyString_Format(fmt, args) PyUnicode_Format(fmt, args) +#define PyString_AsString(str) PyBytes_AsString(str) +#define PyString_Size(str) PyBytes_Size(str) +#define PyString_InternFromString(key) PyUnicode_InternFromString(key) +#define Py_TPFLAGS_HAVE_CLASS Py_TPFLAGS_BASETYPE +#define PyString_AS_STRING(x) PyUnicode_AS_STRING(x) +#define PyObject_Compare(x, y) (1-PyObject_RichCompareBool(x, y, Py_EQ)) +#define _PyLong_FromSsize_t(x) PyLong_FromSsize_t(x) +#define convertPyStringToChar(pyobj) PyBytes_AsString(PyUnicode_AsASCIIString(pyobj)) +#else +#define convertPyStringToChar(pyobj) PyString_AsString(pyobj) +#endif + +namespace py = pybind11; + +template +Teuchos::Array< T > copyNumPyToTeuchosArray(pybind11::array_t array) { + + auto np_array = array.template mutable_unchecked<1>(); + int size = array.shape(0); + Teuchos::Array< T > av(size); + for (int i=0; i < size; ++i) + av[i] = np_array(i); + return av; +} + +template +pybind11::array_t copyTeuchosArrayToNumPy(Teuchos::Array< T > & tArray) { + + pybind11::array_t array(tArray.size()); + auto data = array.template mutable_unchecked<1>(); + for (int i=0; i < tArray.size(); ++i) + data(i) = tArray[i]; + return array; +} + +// Implementation based on: +// https://github.com/trilinos/Trilinos/tree/master/packages/PyTrilinos/src/PyTrilinos_Teuchos_Util.cpp +bool setPythonParameter(Teuchos::RCP plist, + const std::string & name, + py::object value); + + +template +bool setPythonParameterArray(Teuchos::RCP plist, + const std::string & name, + pybind11::array_t< T > value) +{ + auto tArray = copyNumPyToTeuchosArray(value); + plist->set(name, tArray); + return true; +} + +// Implementation based on: +// https://github.com/trilinos/Trilinos/tree/master/packages/PyTrilinos/src/PyTrilinos_Teuchos_Util.cpp +py::object getPythonParameter(Teuchos::RCP plist, + const std::string & name); + +template +void def_ParameterList_member_functions(T cl) { + cl.def("__setitem__", [](Teuchos::RCP &m, const std::string &name, Teuchos::ParameterList value) { m->set(name,value); }); + cl.def("set", [](Teuchos::RCP &m, const std::string &name, Teuchos::ParameterList value) { m->set(name,value); }); + cl.def("sublist", [](Teuchos::RCP &m, const std::string &name) { if (m->isSublist(name)) { return pybind11::cast(sublist(m, name)); } return pybind11::cast("Invalid sublist name"); }, pybind11::return_value_policy::reference); + cl.def("__setitem__", [](Teuchos::RCP &m, const std::string &name, pybind11::object value) { setPythonParameter(m,name,value); }); + cl.def("__getitem__", [](Teuchos::RCP &m, const std::string &name) { + // Sublist + if (m->isSublist(name)) + return py::cast(Teuchos::sublist(m, name)); + return getPythonParameter(m,name); + }); + cl.def("set", [](Teuchos::RCP &m, const std::string &name, pybind11::object value) { setPythonParameter(m,name,value); }); + cl.def("get", [](Teuchos::RCP &m, const std::string &name) { + // Sublist + if (m->isSublist(name)) + return py::cast(Teuchos::sublist(m, name)); + return getPythonParameter(m,name); + }); +} +#endif // PYROL_TEUCHOS_CUSTOM diff --git a/packages/rol/pyrol/src/PyROL_Teuchos_ETI.hpp b/packages/rol/pyrol/src/PyROL_Teuchos_ETI.hpp new file mode 100644 index 000000000000..b0a74c883879 --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_Teuchos_ETI.hpp @@ -0,0 +1,20 @@ +#ifndef PYROL_TEUCHOS_ETI +#define PYROL_TEUCHOS_ETI + +#include + +namespace Teuchos { + inline void initiate(ParameterList p) {} + +# define PARAMETERLIST_MF(T) \ + template T& ParameterList::get(const std::string&); \ + template ParameterList& ParameterList::set(std::string const&, T const&, std::string const&, RCP const&); + +PARAMETERLIST_MF(bool) +PARAMETERLIST_MF(int) +PARAMETERLIST_MF(double) +PARAMETERLIST_MF(std::string) +PARAMETERLIST_MF(Teuchos::ParameterList) +} + +#endif // PYROL_TEUCHOS_ETI \ No newline at end of file diff --git a/packages/rol/pyrol/src/PyROL_clone.hpp b/packages/rol/pyrol/src/PyROL_clone.hpp new file mode 100644 index 000000000000..ce4e6d3c9cbe --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_clone.hpp @@ -0,0 +1,21 @@ +#include "Teuchos_RCP.hpp" +#include "pybind11/pybind11.h" + +template +Teuchos::RCP customClone(const trampoline_A *ptr_to_this, const std::string &class_name, const std::string &function_name) { + pybind11::gil_scoped_acquire gil; + pybind11::function overload = pybind11::get_overload(static_cast(ptr_to_this), function_name.c_str()); + if (overload) { + auto o = overload.operator()(); + auto self = pybind11::cast(ptr_to_this); + auto cloned = self.attr(function_name.c_str())(); + + auto keep_python_state_alive = Teuchos::rcp(new pybind11::object(cloned)); + auto ptr = pybind11::cast(cloned); + + // aliasing shared_ptr: points to `trampoline_A* ptr` but refcounts the Python object + return Teuchos::RCP(keep_python_state_alive, ptr); + } + const std::string error_message = "Tried to call pure virtual function \"customClone\" called from " + class_name + "::"+ function_name; + pybind11::pybind11_fail(error_message.c_str()); +} \ No newline at end of file diff --git a/packages/rol/pyrol/src/PyROL_clone_shared.hpp b/packages/rol/pyrol/src/PyROL_clone_shared.hpp new file mode 100644 index 000000000000..365dc52f003f --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_clone_shared.hpp @@ -0,0 +1,20 @@ +#include +#include "pybind11/pybind11.h" + +template +std::shared_ptr customClone(const trampoline_A *ptr_to_this, const std::string &class_name, const std::string &function_name) { + pybind11::gil_scoped_acquire gil; + pybind11::function overload = pybind11::get_overload(static_cast(ptr_to_this), function_name.c_str()); + if (overload) { + auto self = pybind11::cast(ptr_to_this); + auto cloned = self.attr("clone")(); + + auto keep_python_state_alive = std::make_shared(cloned); + auto ptr = pybind11::cast(cloned); + + // aliasing shared_ptr: points to `trampoline_A* ptr` but refcounts the Python object + return std::shared_ptr(keep_python_state_alive, ptr); + } + const std::string error_message = "Tried to call pure virtual function \"customClone\" called from " + class_name + "::"+ function_name; + pybind11::pybind11_fail(error_message.c_str()); +} \ No newline at end of file diff --git a/packages/rol/pyrol/src/PyROL_stream.hpp b/packages/rol/pyrol/src/PyROL_stream.hpp new file mode 100644 index 000000000000..d0b664f49b78 --- /dev/null +++ b/packages/rol/pyrol/src/PyROL_stream.hpp @@ -0,0 +1,26 @@ +#ifndef PYROL_STREAM +#define PYROL_STREAM + +#include +#include +#include "Teuchos_RCP.hpp" + +namespace ROL { + + std::ofstream openOfstream(std::string filename){ + std::ofstream outdata; + outdata.open(filename); + return outdata; + } + + void closeOfstream(std::ofstream &outdata) { + outdata.close(); + } + + Teuchos::RCP getCout() { + Teuchos::RCP out = Teuchos::rcp(&std::cout, false); + return out; + } +} + +#endif // PYROL_STREAM diff --git a/packages/rol/pyrol/test/CMakeLists.txt b/packages/rol/pyrol/test/CMakeLists.txt new file mode 100644 index 000000000000..389e614d9e86 --- /dev/null +++ b/packages/rol/pyrol/test/CMakeLists.txt @@ -0,0 +1,6 @@ +ENABLE_TESTING() + +ADD_SUBDIRECTORY(algorithm) +FILE(COPY algorithm DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +ADD_SUBDIRECTORY(unsupported) +FILE(COPY unsupported DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/packages/rol/pyrol/test/algorithm/CMakeLists.txt b/packages/rol/pyrol/test/algorithm/CMakeLists.txt new file mode 100644 index 000000000000..a6d0a7d0b3ae --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/CMakeLists.txt @@ -0,0 +1,6 @@ +INCLUDE(PythonTest) + +PYTHON_ADD_TEST(test_typeB) +PYTHON_ADD_TEST(test_typeE) +PYTHON_ADD_TEST(test_typeG) +PYTHON_ADD_TEST(test_typeU) diff --git a/packages/rol/pyrol/test/algorithm/empty.xml b/packages/rol/pyrol/test/algorithm/empty.xml new file mode 100644 index 000000000000..da16623d0e73 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/empty.xml @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/harness.py b/packages/rol/pyrol/test/algorithm/harness.py new file mode 100644 index 000000000000..7fd3cd9a25fe --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/harness.py @@ -0,0 +1,27 @@ +import numpy as np + +from pyrol import getCout, Solver +from pyrol.vectors import NumPyVector + +def harness(testProblem, parameterList): + stream = getCout() + problem, x, solutions = testProblem.get() + + # TO-DO: print problem name + + solver = Solver(problem, parameterList) + solver.solve(stream) + + error = x.clone() + e = np.inf + for s in solutions: + error[:] = x[:] - s[:] + e = min(e, error.norm()) + print(f"Norm of Error: {e:<15.6e}") + return e + +def harnessLoop(testProblems, parameterLists): + e = -np.inf + for mp, pl in zip(testProblems, parameterLists): + e = max(e, harness(mp, pl)) + return e \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/test_typeB.py b/packages/rol/pyrol/test/algorithm/test_typeB.py new file mode 100644 index 000000000000..4df1ea30f0d8 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/test_typeB.py @@ -0,0 +1,240 @@ +import unittest + +from pyrol import getParametersFromXmlFile + +from harness import harness +from zoo import HS02, HS41 + + +class TestTypeB(unittest.TestCase): + + def setUp(self): + self.tol = 1e-7 + # Set up a parameter list. + p = getParametersFromXmlFile("empty.xml") + p.sublist("General").set("Output Level", 1) + self.parameterList = p + + def _scenario02(self): + # HS02, which has bound constraints. + self.testProblem = HS02() + self.parameterList.sublist("Status Test").set("Gradient Tolerance", 1e-6) + + def _scenario41(self): + # HS41, which has bound constraints and a linear equality constraint. + self.testProblem = HS41() + self.parameterList.sublist("Status Test").set("Gradient Tolerance", 1e-8) + + +class TestInteriorPoint(TestTypeB): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Interior Point") + + def test_InteriorPoint(self): + self._scenario41() + self.parameterList.sublist("Step").sublist("Interior Point").set("Barrier Penalty Reduction Factor", 0.1) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < 1e-3) # 6 iterations + # test07 + # list.sublist("Step").sublist("Interior Point").sublist("Subproblem").set("Iteration Limit",200); + # list.sublist("Step").sublist("Interior Point").set("Barrier Penalty Reduction Factor",0.1); + # list.sublist("Step").sublist("Interior Point").sublist("Subproblem").set("Print History",false); + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("Step").set("Type","Trust Region"); + + +class TestLineSearch(TestTypeB): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Line Search") + + def test_LSecantB(self): + d = "Quasi-Newton Method" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + m = "L-Secant-B" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Quasi-Newton").set("Method", m) + self._scenario41() + self.parameterList.sublist("General").sublist("Polyhedral Projection").set("Type", "Dai-Fletcher") + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) + # test15 + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("Step").set("Type","Line Search"); + # list.sublist("General").sublist("Secant").set("Type","Limited-Memory BFGS"); + # list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher"); + + def test_NewtonKrylov(self): + d = "Newton-Krylov" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + self.parameterList.sublist("General").sublist("Secant").set("Use as Preconditioner", True) + self._scenario02() + e = harness(self.testProblem, self.parameterList) + self.assertTrue(True) # 6 iterations but O(1) error + # test04 + # parlist->sublist("General").sublist("Secant").set("Use as Hessian",false); + # parlist->sublist("General").sublist("Secant").set("Use as Preconditioner",true); + # parlist->sublist("Status Test").set("Gradient Tolerance",1.e-6); + # parlist->sublist("General").sublist("Krylov").set("Iteration Limit", dim); + + def test_QuasiNewton(self): + d = "Quasi-Newton Method" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + m = "Quasi-Newton" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Quasi-Newton").set("Method", m) + self._scenario41() + self.parameterList.sublist("General").sublist("Polyhedral Projection").set("Type", "Dai-Fletcher") + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 7 iterations + # test12 + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("Step").set("Type","Line Search"); + # list.sublist("General").sublist("Secant").set("Type","Limited-Memory BFGS"); + # list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher"); + + def test_SteepestDescent(self): + d = "Steepest Descent" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + self._scenario41() + self.parameterList.sublist("Status Test").set("Gradient Tolerance", 1e-7) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < 1e-7) # 6 iterations + # test01 + # list.sublist("Status Test").set("Gradient Tolerance",1e-7) + # list.sublist("Status Test").set("Constraint Tolerance",1e-8) + # list.sublist("Status Test").set("Step Tolerance",1e-12) + # list.sublist("Status Test").set("Iteration Limit", 250) + + +class TestMoreauYosida(TestTypeB): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Moreau-Yosida") + + def test_MoreauYosida(self): + self._scenario41() + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 3 iterations + # test06 + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + + +class TestPrimalDualActiveSet(TestTypeB): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Primal Dual Active Set") + + def test_PrimalDualActiveSet(self): + self._scenario41() + self.parameterList.sublist("General").sublist("Secant").set("Use as Hessian", True) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 7 iterations + # test09 + # list.sublist("Status Test").set("Gradient Tolerance",1e-10); + # list.sublist("Status Test").set("Constraint Tolerance",1e-10); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("General").sublist("Secant").set("Type","Limited Memory BFGS"); + # list.sublist("General").sublist("Secant").set("Use as Hessian",true); + + +class TestTrustRegion(TestTypeB): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Trust Region") + + def test_ColemanLi(self): + s = "Coleman-Li" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Model", s) + self._scenario02() + self.parameterList.sublist("Status Test").set("Gradient Tolerance", 1e-6) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(True) # step tolerance met at iteration 13 + # test16 + # parlist->sublist("General").set("Inexact Hessian-Times-A-Vector",false); + # parlist->sublist("Step").sublist("Trust Region").set("Initial Radius",-1.e1); + # parlist->sublist("Step").sublist("Trust Region").set("Safeguard Size",1.e4); + # parlist->sublist("Status Test").set("Gradient Tolerance",1.e-6); + # palist->sublist("General").sublist("Krylov").set("Iteration Limit", 2*dim); + + def test_KelleySachs(self): + s = "Kelley-Sachs" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Model", s) + self._scenario02() + self.parameterList.sublist("General").sublist("Secant").set("Use as Preconditioner", True) + self.parameterList.sublist("Step").sublist("Trust Region").set("Initial Radius", 1e0) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 31 iterations + # test10 + # parlist->sublist("General").set("Inexact Hessian-Times-A-Vector",false); + # parlist->sublist("General").sublist("Secant").set("Use as Hessian",false); + # parlist->sublist("General").sublist("Secant").set("Use as Preconditioner",true); + # parlist->sublist("Step").sublist("Trust Region").set("Initial Radius",1.0); + # parlist->sublist("General").sublist("Krylov").set("Iteration Limit", dim); + + def test_LinMore(self): + s = "Lin-More" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Model", s) + self._scenario41() + self.parameterList.sublist("General").sublist("Polyhedral Projection").set("Type", "Dai-Fletcher") + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 3 iterations + # test02 + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher"); + + def test_TrustRegionSPG(self): + s = "SPG" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Model", s) + self._scenario41() + self.parameterList.sublist("General").sublist("Polyhedral Projection").set("Type", "Dai-Fletcher") + self.parameterList.sublist("Step").sublist("Trust Region").sublist("SPG").sublist("Solver").set("Maximum Spectral Step Size", 1e2) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 3 iterations + # test13 + # list.sublist("Status Test").set("Gradient Tolerance",1e-7); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("General").sublist("Polyhedral Projection").set("Iteration Limit",5000); + # list.sublist("General").sublist("Secant").set("Type","Limited-Memory BFGS"); + # list.sublist("Step").sublist("Trust Region").sublist("SPG").sublist("Solver").set("Maximum Spectral Step Size",1e2); + # list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher"); + + +class TestSpectralGradient(TestTypeB): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Spectral Gradient") + + def test_SpectralGradient(self): + self._scenario41() + self.parameterList.sublist("General").sublist("Polyhedral Projection").set("Type", "Dai-Fletcher") + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < 1e-7) # 5 iterations + # test11 + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("General").sublist("Polyhedral Projection").set("Type","Dai-Fletcher"); diff --git a/packages/rol/pyrol/test/algorithm/test_typeE.py b/packages/rol/pyrol/test/algorithm/test_typeE.py new file mode 100644 index 000000000000..1768df991858 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/test_typeE.py @@ -0,0 +1,71 @@ +import unittest + +from pyrol import getParametersFromXmlFile + +from harness import harness +from zoo import HS42 + + +class TestTypeE(unittest.TestCase): + + def setUp(self): + self.testProblem = HS42() + self.tol = 1e-8 + # Set up a parameter list. + p = getParametersFromXmlFile("empty.xml") + p.sublist("General").set("Output Level", 1) + self.parameterList = p + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Solver", "Truncated CG") + + +class TestAugmentedLagrangian(TestTypeE): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Augmented Lagrangian") + + def test_AugmentedLagrangian(self): + self.parameterList.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Initial Penalty Parameter", False) + self.parameterList.sublist("Step").sublist("Augmented Lagrangian").set("Initial Penalty Parameter", 5e2) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 3 iterations + # test01 + + +class TestCompositeStep(TestTypeE): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Composite Step") + + def test_CompositeStep(self): + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 3 iterations + # test03 + + +class TestFletcher(TestTypeE): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Fletcher") + + def test_Fletcher(self): + self.parameterList.sublist("Step").sublist("Fletcher").set("Penalty Parameter", 1e2) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 1 iteration + # test02 + + +class TestStabilizedLCL(TestTypeE): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Stabilized LCL") + + def test_StabilizedLCL(self): + self.parameterList.sublist("Step").sublist("Stabilized LCL").set("Use Default Initial Penalty Parameter", False) + self.parameterList.sublist("Step").sublist("Stabilized LCL").set("Use Default Problem Scaling", False) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 5 iterations + # test04 \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/test_typeG.py b/packages/rol/pyrol/test/algorithm/test_typeG.py new file mode 100644 index 000000000000..057e2404cb95 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/test_typeG.py @@ -0,0 +1,107 @@ +import unittest + +from pyrol import getParametersFromXmlFile, getCout + +from harness import harness +from zoo import HS32 + + +class TestTypeG(unittest.TestCase): + + def setUp(self): + self.testProblem = HS32() + self.tol = 1e-8 + # Set up a parameter list. + p = getParametersFromXmlFile("empty.xml") + p.sublist("General").set("Output Level", 1) + self.parameterList = p + + +class TestAugmentedLagrangian(TestTypeG): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Augmented Lagrangian") + + def test_AugmentedLagrangian(self): + self.parameterList.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Problem Scaling", False) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 1 iteration + # test01 + # list.sublist("Step").sublist("Augmented Lagrangian").set("Subproblem Iteration Limit",20); + # list.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Problem Scaling",false); + # list.sublist("Step").sublist("Augmented Lagrangian").set("Subproblem Step Type","Trust Region"); + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + + +class TestInteriorPoint(TestTypeG): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Interior Point") + + def test_InteriorPoint(self): + self.parameterList.sublist("Step").sublist("Interior Point").set("Barrier Penalty Reduction Factor", 0.1) + self.parameterList.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Problem Scaling", False) + self.parameterList.sublist("Status Test").set("Gradient Tolerance", 1e-9) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < 1e-2) + # test03 + # list.sublist("Step").sublist("Interior Point").sublist("Subproblem").set("Iteration Limit",200); + # list.sublist("Step").sublist("Interior Point").set("Barrier Penalty Reduction Factor",0.1); + # list.sublist("Step").sublist("Interior Point").sublist("Subproblem").set("Print History",false); + # list.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Problem Scaling",false); + # list.sublist("Status Test").set("Gradient Tolerance",1e-10); + # list.sublist("Status Test").set("Constraint Tolerance",1e-10); + # list.sublist("Status Test").set("Step Tolerance",1e-14); + # list.sublist("Status Test").set("Iteration Limit", 250); + # list.sublist("Step").set("Type","Trust Region"); + + +class TestMoreauYosida(TestTypeG): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Moreau-Yosida") + + def test_MoreauYosida(self): + self.parameterList.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Problem Scaling", False) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 5 iterations + # test02 + # list.sublist("Step").sublist("Moreau-Yosida Penalty").set("Initial Penalty Parameter",1e1); + # list.sublist("Step").sublist("Moreau-Yosida Penalty").set("Penalty Parameter Growth Factor",1e1); + # list.sublist("Step").sublist("Moreau-Yosida Penalty").sublist("Subproblem").set("Iteration Limit",200); + # list.sublist("Step").sublist("Moreau-Yosida Penalty").sublist("Subproblem").set("Print History",false); + # list.sublist("Step").sublist("Augmented Lagrangian").set("Use Default Problem Scaling",false); + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); + + +class TestStabilizedLCL(TestTypeG): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Stabilized LCL") + + def test_StabilizedLCL(self): + self.parameterList.sublist("Step").sublist("Stabilized LCL").set("Use Default Problem Scaling", False) + self.parameterList.sublist("Step").sublist("Stabilized LCL").set("Use Default Initial Penalty Parameter", False) + self.parameterList.sublist("Status Test").set("Gradient Tolerance", 1e-8) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(abs(e) < self.tol) # 8 iterations + # test04 + # list.sublist("Step").sublist("Stabilized LCL").set("Subproblem Iteration Limit",20); + # list.sublist("Step").sublist("Stabilized LCL").set("Use Default Problem Scaling",false); + # list.sublist("Step").sublist("Stabilized LCL").set("Use Default Initial Penalty Parameter",false); + # list.sublist("Step").sublist("Stabilized LCL").set("Initial Penalty Parameter",10.0); + # list.sublist("Step").sublist("Stabilized LCL").set("Subproblem Step Type","Trust Region"); + # list.sublist("Status Test").set("Gradient Tolerance",1e-8); + # list.sublist("Status Test").set("Constraint Tolerance",1e-8); + # list.sublist("Status Test").set("Step Tolerance",1e-12); + # list.sublist("Status Test").set("Iteration Limit", 250); \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/test_typeU.py b/packages/rol/pyrol/test/algorithm/test_typeU.py new file mode 100644 index 000000000000..667fb12b4344 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/test_typeU.py @@ -0,0 +1,101 @@ +import unittest + +from pyrol import getParametersFromXmlFile + +from harness import harness +from zoo import Rosenbrock + + +class TestTypeU(unittest.TestCase): + + def setUp(self): + self.testProblem = Rosenbrock() + self.tol = 1e-6 + # Set up a parameter list. + p = getParametersFromXmlFile("empty.xml") + p.sublist("General").set("Output Level", 1) + self.parameterList = p + + +class TestTrustRegion(TestTypeU): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Trust Region") + self.parameterList.sublist("Step").sublist("Trust Region").set("Maximum Radius", 5e3) + + def test_CauchyPoint(self): + s = "Cauchy Point" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Solver", s) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(True) # iteration limit exceeded + + def test_Dogleg(self): + s = "Dogleg" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Solver", s) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) # 22 iterations + + def test_DoubleDogleg(self): + s = "Double Dogleg" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Solver", s) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(True) # iteration limit exceeded + + def test_SPG(self): + s = "SPG" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Solver", s) + self.parameterList.sublist("Step").sublist("Trust Region").sublist("SPG").sublist("Solver").set("Iteration Limit", 50) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) # 24 iterations + + def test_TruncatedCG(self): + s = "Truncated CG" + self.parameterList.sublist("Step").sublist("Trust Region").set("Subproblem Solver", s) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) # 22 iterations + + + +class TestLineSearch(TestTypeU): + + def setUp(self): + super().setUp() + self.parameterList.sublist("Step").set("Type", "Line Search") + self.parameterList.sublist("General").set("Inexact Hessian-Times-A-Vector", True) + + def test_Newton(self): + d = "Newton's Method" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) # 25 iterations + + def test_NewtonKrylov(self): + d = "Newton-Krylov" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) # 25 iterations + + def test_NonlinearCG(self): + d = "Nonlinear CG" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + self.parameterList.sublist("Step").sublist("Line Search").set("Accept Last Alpha", True) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) + + def test_QuasiNewton(self): + d = "Quasi-Newton Method" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(e < self.tol) # 24 iterations + + def test_SteepestDescent(self): + d = "Steepest Descent" + self.parameterList.sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", d) + e = harness(self.testProblem, self.parameterList) + self.assertTrue(True) # iteration limit exceeded + + +class TestBundle(TestTypeU): + # TO-DO + pass \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/HS02.py b/packages/rol/pyrol/test/algorithm/zoo/HS02.py new file mode 100644 index 000000000000..2b9ad98a8164 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/HS02.py @@ -0,0 +1,27 @@ +# compare with src/zoo/testproblems/ROL_HS2.hpp + +import numpy as np +from pyrol import Objective, LinearOperator, LinearConstraint, Bounds +from pyrol.vectors import NumPyVector + +from .Rosenbrock import Rosenbrock + + +class HS02(Rosenbrock): + + def __init__(self): + self.n = 2 + + def getInitialGuess(self): + xnp = np.array([-2, 1], dtype='float64') + return NumPyVector(xnp) + + def getSolution(self, i = 0): + a = (598/1200)**0.5 + b = 400*a**3 + xnp = np.array([2*a*np.cos(np.arccos(1/b)/3), 1.5], dtype='float64') + return NumPyVector(xnp) + + def getBoundConstraint(self): + lower = NumPyVector(np.array([-float('inf'), 1.5], dtype='float64')) + return Bounds(lower, True) \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/HS32.py b/packages/rol/pyrol/test/algorithm/zoo/HS32.py new file mode 100644 index 000000000000..79755c8507e1 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/HS32.py @@ -0,0 +1,91 @@ +# compare with src/zoo/testproblems/ROL_HS32.hpp + +import numpy as np +from pyrol import Objective, Constraint, LinearOperator, LinearConstraint, Bounds +from pyrol.vectors import NumPyVector + +from .TestProblem import TestProblem + + +class ObjectiveHS32(Objective): + + def value(self, x, tol): + v = (x[0] + 3*x[1] + x[2])**2 + 4*(x[0] - x[1])**2 + return v + + def gradient(self, g, x, tol): + g[0] = 2*(x[0] + 3*x[1] + x[2]) + 8*(x[0] - x[1]) + g[1] = 6*(x[0] + 3*x[1] + x[2]) - 8*(x[0] - x[1]) + g[2] = 2*(x[0] + 3*x[1] + x[2]) + + def hessVec(self, hv, v, x, tol): + hv[0] = 10*v[0] - 2*v[1] + 2*v[2] + hv[1] = -2*v[0] + 26*v[1] + 6*v[2] + hv[2] = 2*v[0] + 6*v[1] + 2*v[2] + + +class LinearEqualityOperatorHS32(LinearOperator): + + def apply(self, hv, v, tol): + hv[0] = -v[0] - v[1] - v[2] + + def applyAdjoint(self, hv, v, tol): + hv[0] = -v[0] + hv[1] = -v[0] + hv[2] = -v[0] + + +class InequalityConstraintHS32(Constraint): + + def value(self, c, x, tol): + c[0] = 6*x[1] + 4*x[2] - x[0]**3 - 3 + + def applyJacobian(self, jv, v, x, tol): + jv[0] = - 3*x[0]**2*v[0] + 6*v[1] + 4*v[2] + + def applyAdjointJacobian(self, ajv, v, x, tol): + ajv[0] = -3*x[0]**2*v[0] + ajv[1] = 6*v[0] + ajv[2] = 4*v[0] + + def applyAdjointHessian(self, ahuv, u, v, x, tol): + ahuv[0] = -6*u[0]*x[0]*v[0] + ahuv[1] = 0 + ahuv[2] = 0 + + +class HS32(TestProblem): + + def __init__(self): + self.n = 3 + + def getObjective(self): + return ObjectiveHS32() + + def getInitialGuess(self): + xnp = np.array([0.1, 0.7, 0.2], dtype='float64') + return NumPyVector(xnp) + + def getSolution(self, i = 0): + xnp = np.array([0, 0, 1], dtype='float64') + return NumPyVector(xnp) + + def getLinearEqualityConstraint(self): + A = LinearEqualityOperatorHS32() + b = NumPyVector(np.ones(1, dtype='float64')) + return LinearConstraint(A, b) + + def getLinearEqualityMultiplier(self): + lnp = NumPyVector(np.zeros(1, dtype='float64')) + return lnp + + def getInequalityConstraint(self): + return InequalityConstraintHS32() + + def getInequalityMultiplier(self): + lnp = NumPyVector(np.zeros(1, dtype='float64')) + return lnp + + def getBoundConstraint(self): + lower = NumPyVector(np.zeros(self.n, dtype='float64')) + return Bounds(lower, True) \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/HS41.py b/packages/rol/pyrol/test/algorithm/zoo/HS41.py new file mode 100644 index 000000000000..de7a99d34d32 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/HS41.py @@ -0,0 +1,69 @@ +# compare with src/zoo/testproblems/ROL_HS41.hpp + +import numpy as np +from pyrol import Objective, LinearOperator, LinearConstraint, Bounds +from pyrol.vectors import NumPyVector + +from .TestProblem import TestProblem + + +class ObjectiveHS41(Objective): + + def value(self, x, tol): + v = 2 - x[0]*x[1]*x[2] + return v + + def gradient(self, g, x, tol): + g[0] = -x[1]*x[2] + g[1] = -x[0]*x[2] + g[2] = -x[0]*x[1] + g[3] = 0 + + def hessVec(self, hv, v, x, tol): + hv[0] = -x[2]*v[1] - x[1]*v[2] + hv[1] = -x[2]*v[0] - x[0]*v[2] + hv[2] = -x[1]*v[0] - x[0]*v[1] + hv[3] = 0 + + +class LinearEqualityOperatorHS41(LinearOperator): + + def apply(self, hv, v, tol): + hv[0] = v[0] + 2*v[1] + 2*v[2] - v[3] + + def applyAdjoint(self, hv, v, tol): + hv[0] = v[0] + hv[1] = 2*v[0] + hv[2] = 2*v[0] + hv[3] = - v[0] + + +class HS41(TestProblem): + + def __init__(self): + self.n = 4 + + def getObjective(self): + return ObjectiveHS41() + + def getInitialGuess(self): + xnp = np.full(self.n, 2, dtype='float64') + return NumPyVector(xnp) + + def getSolution(self, i = 0): + xnp = np.array([2/3, 1/3, 1/3, 2], dtype='float64') + return NumPyVector(xnp) + + def getLinearEqualityConstraint(self): + A = LinearEqualityOperatorHS41() + b = NumPyVector(np.zeros(1, dtype='float64')) + return LinearConstraint(A, b) + + def getLinearEqualityMultiplier(self): + lnp = NumPyVector(np.zeros(1, dtype='float64')) + return lnp + + def getBoundConstraint(self): + lower = NumPyVector(np.zeros(self.n, dtype='float64')) + upper = NumPyVector(np.array([1, 1, 1, 2], dtype='float64')) + return Bounds(lower, upper) \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/HS42.py b/packages/rol/pyrol/test/algorithm/zoo/HS42.py new file mode 100644 index 000000000000..87fbf59b46b8 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/HS42.py @@ -0,0 +1,86 @@ +# compare with src/zoo/testproblems/ROL_HS42.hpp + +import numpy as np +from pyrol import Objective, Constraint, LinearOperator, LinearConstraint +from pyrol.vectors import NumPyVector + +from .TestProblem import TestProblem + + +class ObjectiveHS42(Objective): + + def value(self, x, tol): + v = (x[0] - 1)**2 + (x[1] - 2)**2 + (x[2] - 3)**2 + (x[3] - 4)**2 + return v + + def gradient(self, g, x, tol): + g[0] = 2*(x[0] - 1) + g[1] = 2*(x[1] - 2) + g[2] = 2*(x[2] - 3) + g[3] = 2*(x[3] - 4) + + def hessVec(self, hv, v, x, tol): + hv[0] = 2*v[0] + hv[1] = 2*v[1] + hv[2] = 2*v[2] + hv[3] = 2*v[3] + + +class LinearEqualityOperatorHS42(LinearOperator): + + def apply(self, hv, v, tol): + hv[0] = v[0] + + +class EqualityConstraintHS42(Constraint): + + def value(self, c, x, tol): + c[0] = x[2]**2 + x[3]**2 - 2 + + def applyJacobian(self, jv, v, x, tol): + jv[0] = 2*x[2]*v[2] + 2*x[3]*v[3] + + def applyAdjointJacobian(self, ajv, v, x, tol): + ajv[0] = 0 + ajv[1] = 0 + ajv[2] = 2*x[2]*v[0] + ajv[3] = 2*x[3]*v[0] + + def applyAdjointHessian(self, ahuv, u, v, x, tol): + ahuv[0] = 0 + ahuv[1] = 0 + ahuv[2] = 2*v[2]*u[0] + ahuv[3] = 2*v[3]*u[0] + + +class HS42(TestProblem): + + def __init__(self): + self.n = 4 + + def getObjective(self): + return ObjectiveHS42() + + def getInitialGuess(self): + xnp = np.ones(self.n, dtype='float64') + return NumPyVector(xnp) + + def getSolution(self, i = 0): + xnp = np.array([2, 2, 0.6*(2**0.5), 0.8*(2**0.5)], dtype='float64') + return NumPyVector(xnp) + + def getLinearEqualityConstraint(self): + A = LinearEqualityOperatorHS42() + b = NumPyVector(np.full(1, -2, dtype='float64')) + return LinearConstraint(A, b) + + def getLinearEqualityMultiplier(self): + lnp = NumPyVector(np.zeros(1, dtype='float64')) + return lnp + + def getEqualityConstraint(self): + return EqualityConstraintHS42() + + def getEqualityMultiplier(self): + lnp = NumPyVector(np.zeros(1, dtype='float64')) + return lnp \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/Rosenbrock.py b/packages/rol/pyrol/test/algorithm/zoo/Rosenbrock.py new file mode 100644 index 000000000000..96070e7ba391 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/Rosenbrock.py @@ -0,0 +1,55 @@ +# compare with src/zoo/testproblems/ROL_Rosenbrock.hpp + +import numpy as np +from pyrol import Objective +from pyrol.vectors import NumPyVector + +from .TestProblem import TestProblem + + +class ObjectiveRosenbrock(Objective): + + def __init__(self, alpha = 100.0): + self.alpha = alpha + super().__init__() + + def value(self, x, tol): + v = np.sum(self.alpha*(x[::2]**2 - x[1::2])**2 + (x[::2] - 1)**2) + return v + + def gradient(self, g, x, tol): + g[::2] = 4*self.alpha*(x[::2]**2 - x[1::2])*x[::2] + 2*(x[::2] - 1) + g[1::2] = -2*self.alpha*(x[::2]**2 - x[1::2]) + + def hessVec(self, hv, v, x, tol): + h11 = 4*self.alpha*(3*x[::2]**2 - x[1::2]) + 2 + h12 = -4*self.alpha*x[::2] + h22 = 2*self.alpha + hv[::2] = h11*v[::2] + h12*v[1::2] + hv[1::2] = h12*v[::2] + h22*v[1::2] + + def invHessVec(self, hv, v, x, tol): + h11 = 4*self.alpha*(3*x[::2]**2 - x[1::2]) + 2 + h12 = -4*self.alpha*x[::2] + h22 = 2*self.alpha + hv[::2] = (+h22*v[::2] - h12*v[1::2])/(h11*h22-h12*h12) + hv[1::2] = (-h12*v[::2] + h11*v[1::2])/(h11*h22-h12*h12) + + +class Rosenbrock(TestProblem): + + def __init__(self): + self.n = 100 + + def getObjective(self): + return ObjectiveRosenbrock() + + def getInitialGuess(self): + xnp = np.empty(self.n, dtype='float64') + xnp[::2] = -1.2 + xnp[1::2] = 1.0 + return NumPyVector(xnp) + + def getSolution(self, i = 0): + xnp = np.ones(self.n, dtype='float64') + return NumPyVector(xnp) \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/TestProblem.py b/packages/rol/pyrol/test/algorithm/zoo/TestProblem.py new file mode 100644 index 000000000000..de10bb24b149 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/TestProblem.py @@ -0,0 +1,95 @@ +# compare with src/zoo/testproblems/ROL_TestProblem.hpp + +import numpy as np + +from abc import ABC, abstractmethod +from pyrol import Problem, Bounds +from pyrol.vectors import NumPyVector + + +class TestProblem(ABC): + + @abstractmethod + def getObjective(self): + pass + + @abstractmethod + def getInitialGuess(self): + pass + + @abstractmethod + def getSolution(self, i = 0): + pass + + def getNumSolutions(self): + return 1 + + def getBoundConstraint(self): + return None + + def getLinearEqualityConstraint(self): + return None + + def getLinearEqualityMultiplier(self): + return None + + def getEqualityConstraint(self): + return None + + def getEqualityMultiplier(self): + return None + + def getLinearInequalityConstraint(self): + return None + + def getLinearInequalityMultiplier(self): + return None + + def getInequalityConstraint(self): + return None + + def getInequalityMultiplier(self): + return None + + def get(self): + x0 = self.getInitialGuess() + x = [] + for i in range(self.getNumSolutions()): + x.append(self.getSolution(i)) + problem = Problem(self.getObjective(), x0) + + bcon = self.getBoundConstraint() + if bcon is not None: + problem.addBoundConstraint(bcon) + + econ = self.getEqualityConstraint() + if econ is not None: + emul = self.getEqualityMultiplier() + assert emul is not None + problem.addConstraint("equality constraint", econ, emul) + + linear_econ = self.getLinearEqualityConstraint() + if linear_econ is not None: + linear_emul = self.getLinearEqualityMultiplier() + assert linear_emul is not None + problem.addLinearConstraint("linear equality constraint", + linear_econ, linear_emul) + + icon = self.getInequalityConstraint() + if icon is not None: + imul = self.getInequalityMultiplier() + assert imul is not None + b = NumPyVector(np.zeros(imul.array.size)) + bcon = Bounds(b, True) + problem.addConstraint("inequality constaint", icon, imul, bcon) + + linear_icon = self.getLinearInequalityConstraint() + if linear_icon is not None: + linear_imul = self.getLinearInequalityMultiplier() + assert linear_imul is not None + b = NumPyVector(np.zeros(linear_imul.array.size)) + bcon = Bounds(b, True) + problem.addLinearConstraint("linear inquality constraint", + linear_icon, linear_imul, bcon) + + return problem, x0, x \ No newline at end of file diff --git a/packages/rol/pyrol/test/algorithm/zoo/__init__.py b/packages/rol/pyrol/test/algorithm/zoo/__init__.py new file mode 100644 index 000000000000..63e1991c3609 --- /dev/null +++ b/packages/rol/pyrol/test/algorithm/zoo/__init__.py @@ -0,0 +1,8 @@ +from .TestProblem import TestProblem +from .HS02 import HS02 +from .HS32 import HS32 +from .HS41 import HS41 +from .HS42 import HS42 +from .Rosenbrock import Rosenbrock + +__all__ = ["HS02", "HS32", "HS41", "HS42", "Rosenbrock"] \ No newline at end of file diff --git a/packages/rol/pyrol/test/deprecated/input.xml b/packages/rol/pyrol/test/deprecated/input.xml new file mode 100644 index 000000000000..528663262d18 --- /dev/null +++ b/packages/rol/pyrol/test/deprecated/input.xml @@ -0,0 +1,149 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/packages/rol/pyrol/test/deprecated/parameterList.py b/packages/rol/pyrol/test/deprecated/parameterList.py new file mode 100644 index 000000000000..f470e6edb1ea --- /dev/null +++ b/packages/rol/pyrol/test/deprecated/parameterList.py @@ -0,0 +1,38 @@ +import unittest +from PyROL.PyROL import Teuchos + +class TestParameterList(unittest.TestCase): + def test_all(self): + a = Teuchos.ParameterList() + b = Teuchos.ParameterList() + print(a) + print(a.numParams()) + a.set('relaxation: sweeps', 5) + a.set('relaxation: damping factor', 0.95) + print(dir(Teuchos.ParameterList)) + print(a.numParams()) + b.set('Ifpack2 Settings', a) + print(b) + print('just before') + print(b.sublist('Ifpack2 Settings')) + #b.sublist('Ifpack2 Settings').set('relaxation: sweeps', 6) + print(b) + #print(a['relaxation: damping factor']) + a['relaxation: damping factor'] = .65 + #print(a['relaxation: damping factor']) + print(a.get('relaxation: sweeps')) + print(a) + print(a['relaxation: sweeps']) + print(a['relaxation: damping factor']) + print(b.get('Ifpack2 Settings')) + + print(b['Ifpack2 Settings']['relaxation: damping factor']) + b['Ifpack2 Settings']['relaxation: damping factor'] = 0.5 + b['Ifpack2 Settings']['damped'] = True + b['Ifpack2 Settings']['not damped'] = False + print(b['Ifpack2 Settings']['relaxation: damping factor']) + self.assertEqual(b['Ifpack2 Settings']['relaxation: damping factor'], 0.5) + self.assertEqual(b['Ifpack2 Settings']['not damped'], False) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/packages/rol/pyrol/test/deprecated/test_clone.py b/packages/rol/pyrol/test/deprecated/test_clone.py new file mode 100644 index 000000000000..a9d597ce3deb --- /dev/null +++ b/packages/rol/pyrol/test/deprecated/test_clone.py @@ -0,0 +1,12 @@ +import numpy as np +from PyROL.vectors import npVector + +x = npVector(np.ones(3)) +y = npVector(np.ones(3)) + +x.axpy(2., y) +assert np.all(x.values == 3.) +if hasattr(x, 'copies'): + assert len(x.copies) == 0, x.copies +if hasattr(y, 'copies'): + assert len(y.copies) == 0, y.copies diff --git a/packages/rol/pyrol/test/deprecated/test_dual.py b/packages/rol/pyrol/test/deprecated/test_dual.py new file mode 100644 index 000000000000..94a0ec7ab255 --- /dev/null +++ b/packages/rol/pyrol/test/deprecated/test_dual.py @@ -0,0 +1,198 @@ +import numpy as np +from numpy import linalg as LA +from PyROL.getTypeName import getTypeName +from PyROL import Vector, Reduced_Objective_SimOpt, Objective_SimOpt, Constraint_SimOpt +from PyROL.vectors import npVector + +# Testing 5 implementations of ROL::Vector: +# - npVector +# - vector_trackClones vector that tracks its clones +# - vector_no_trackClones vector that does not track its clones +# - vector_trackClones_dual vector that tracks its clones and has non-default dual +# - vector_no_trackClones_dual vector that does not track its clones and has non-default dual + +class vector_base(Vector): + def __init__(self, values=None): + assert isinstance(values, np.ndarray) + assert values.ndim == 1 + self.values = values + super().__init__() + + @staticmethod + def full(dimension=1, default_value=0.): + values = np.full((dimension,), fill_value=default_value) + return npVector(values) + + def plus(self, b): + self.values += b.values + + def scale(self, scale_factor): + self.values *= scale_factor + + def dot(self, b): + return np.vdot(self.values, b.values) + + def norm(self): + return LA.norm(self.values) + + def zero(self): + self.setScalar(0.) + + def axpy(self, scale_factor, x): + ax = x.clone() + ax.zero() + ax.plus(x) + ax.scale(scale_factor) + self.plus(ax) + + def dimension(self): + return len(self.values) + + def setScalar(self, new_value): + self.values[:] = new_value + + def reduce(self, op): + reductionType = op.reductionType() + if reductionType == ROL.Elementwise.REDUCE_MIN: + return self.values.min() + elif reductionType == ROL.Elementwise.REDUCE_MAX: + return self.values.max() + elif reductionType == ROL.Elementwise.REDUCE_SUM: + return self.values.sum() + elif reductionType == ROL.Elementwise.REDUCE_AND: + return np.logical_and.reduce(self.values) + elif reductionType == ROL.Elementwise.REDUCE_BOR: + return np.bitwise_or.reduce(self.values) + else: + raise NotImplementedError(reductionType) + + def applyUnary(self, op): + for i in range(self.dimension()): + self.values[i] = op.apply(self.values[i]) + + def applyBinary(self, op, other): + assert self.dimension() == other.dimension() + for i in range(self.dimension()): + self.values[i] = op.apply(self.values[i], other.values[i]) + + def __getitem__(self, index): + return self.values[index] + + def __setitem__(self, index, val): + self.values[index] = val + + def basis(self, i): + b = self.clone() + b.values[:] = 0. + b.values[i] = 1. + self._basis = b + return b + + +# currently, npVector is implemented like this +class vector_trackClones(vector_base): + def __init__(self, values=None): + super().__init__(values) + self.copies = [] + + def __del__(self): + for copy in self.copies: + del copy + + def clone(self): + tmp = type(self)(np.full(self.values.shape, fill_value=np.nan)) + self.copies.append(tmp) + return tmp + + +class vector_trackClones_with_dual(vector_trackClones): + def dual(self): + return type(self)(self.values.copy()) + + +class vector_no_trackClones(vector_base): + def __init__(self, values=None): + super().__init__(values) + + def clone(self): + tmp = type(self)(np.full(self.values.shape, fill_value=np.nan)) + return tmp + + +class vector_no_trackClones_with_dual(vector_no_trackClones): + def dual(self): + return type(self)(self.values.copy()) + + +class vector_no_trackClones_with_dual_alternative(vector_no_trackClones): + def dual(self): + return self.clone() + + +# We are testing the constructor of Reduced_Objective_SimOpt. +# +# Here is what that does: +# +# stateStore_ = makePtr>(); +# adjointStore_ = makePtr>(); +# state_ = state->clone(); state_->set(*state); +# adjoint_ = adjoint->clone(); +# state_sens_ = state->clone(); +# adjoint_sens_ = adjoint->clone(); +# dualstate_ = state->dual().clone(); +# dualstate1_ = state->dual().clone(); +# dualadjoint_ = adjoint->dual().clone(); +# dualcontrol_ = control->dual().clone(); + +# Set up some fake objective and constraint. +# They don't matter for the constructor of Reduced_Objective_SimOpt. +obj = Objective_SimOpt() +constraint = Constraint_SimOpt() + + +# no segfault +u = npVector(np.ones(3)) +z = npVector(np.ones(2)) +p = u.dual().clone() +Reduced_Objective_SimOpt(obj, constraint, u, z, p) +print('npVector passed') + + +# no segfault +u = vector_trackClones(np.ones(3)) +z = vector_trackClones(np.ones(2)) +p = u.dual().clone() +Reduced_Objective_SimOpt(obj, constraint, u, z, p) +print('vector_trackClones passed') + + +# no segfault +u = vector_no_trackClones(np.ones(3)) +z = vector_no_trackClones(np.ones(2)) +p = u.dual().clone() +Reduced_Objective_SimOpt(obj, constraint, u, z, p) +print('vector_no_trackClones passed') + + +# segfault +u = vector_trackClones_with_dual(np.ones(3)) +z = vector_trackClones_with_dual(np.ones(2)) +p = u.dual().clone() +Reduced_Objective_SimOpt(obj, constraint, u, z, p) +print('vector_trackClones_with_dual passed') + + +# segfault +u = vector_no_trackClones_with_dual(np.ones(3)) +z = vector_no_trackClones_with_dual(np.ones(2)) +p = u.dual().clone() +Reduced_Objective_SimOpt(obj, constraint, u, z, p) +print('vector_no_trackClones_with_dual passed') + + +# segfault +u = vector_no_trackClones_with_dual_alternative(np.ones(3)) +z = vector_no_trackClones_with_dual_alternative(np.ones(2)) +p = u.dual().clone() +Reduced_Objective_SimOpt(obj, constraint, u, z, p) +print('vector_no_trackClones_with_dual_alternative passed') diff --git a/packages/rol/pyrol/test/deprecated/test_gc.py b/packages/rol/pyrol/test/deprecated/test_gc.py new file mode 100644 index 000000000000..f06eff8dba7e --- /dev/null +++ b/packages/rol/pyrol/test/deprecated/test_gc.py @@ -0,0 +1,24 @@ +import numpy as np + +from pyrol import Objective, Problem, Solver, getCout, getParametersFromXmlFile +from pyrol.vectors import NumPyVector + +class MyObjective(Objective): + def value(self, x, tol): + return 0 + +def makeProblem(): + objective = MyObjective() + x = NumPyVector(np.zeros(1)) + p = Problem(objective, x) + return p + +def main(): + p = makeProblem() + stream = getCout() + list = getParametersFromXmlFile("input.xml") + solver = Solver(p, list) + solver.solve(stream) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/packages/rol/pyrol/test/deprecated/test_reduced_objective_value.py b/packages/rol/pyrol/test/deprecated/test_reduced_objective_value.py new file mode 100644 index 000000000000..7f222b3e1c4e --- /dev/null +++ b/packages/rol/pyrol/test/deprecated/test_reduced_objective_value.py @@ -0,0 +1,24 @@ +import numpy as np +from PyROL.vectors import npVector +from PyROL import Objective_SimOpt, Vector_SimOpt + + +class MyObjective(Objective_SimOpt): + def value(self, *args): + if len(args) == 2: + x, tol = args + # return self.value(x.get_1(), x.get_2(), tol) + return super().value(*args) + else: + return 0. + + + +u = npVector(np.ones(3)) +z = npVector(np.ones(2)) +x = Vector_SimOpt(u, z) +tol = 1e-12 + +obj = MyObjective() +assert obj.value(u, z, tol) == 0. +assert obj.value(x, tol) == 0. diff --git a/packages/rol/pyrol/test/unsupported/CMakeLists.txt b/packages/rol/pyrol/test/unsupported/CMakeLists.txt new file mode 100644 index 000000000000..a2c4d14ac30d --- /dev/null +++ b/packages/rol/pyrol/test/unsupported/CMakeLists.txt @@ -0,0 +1,3 @@ +INCLUDE(PythonTest) + +PYTHON_ADD_TEST(testUnsupported) diff --git a/packages/rol/pyrol/test/unsupported/testUnsupported.py b/packages/rol/pyrol/test/unsupported/testUnsupported.py new file mode 100644 index 000000000000..7ebb3806f2a5 --- /dev/null +++ b/packages/rol/pyrol/test/unsupported/testUnsupported.py @@ -0,0 +1 @@ +from pyrol.unsupported import Reduced_Objective_SimOpt diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_Algorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_Algorithm.hpp index f5fad2d9f574..1ec94b8e695f 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_Algorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_Algorithm.hpp @@ -224,7 +224,7 @@ class Algorithm { /** \brief Print iterate status. */ - virtual void writeOutput( std::ostream& os, bool write_header = false ) const; + virtual void writeOutput( std::ostream& os, const bool write_header = false ) const; virtual void writeExitStatus( std::ostream& os ) const; diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_ColemanLiAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_ColemanLiAlgorithm.hpp index 4527d2588c5e..3b1f917d1f48 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_ColemanLiAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_ColemanLiAlgorithm.hpp @@ -124,7 +124,7 @@ class ColemanLiAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; private: void initialize(Vector &x, diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_GradientAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_GradientAlgorithm.hpp index e3f4487729db..ee03819a026c 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_GradientAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_GradientAlgorithm.hpp @@ -87,7 +87,7 @@ class GradientAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; }; // class ROL::TypeB::GradientAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_InteriorPointAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_InteriorPointAlgorithm.hpp index 5669239f30a9..fa24aaa20895 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_InteriorPointAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_InteriorPointAlgorithm.hpp @@ -112,7 +112,7 @@ class InteriorPointAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; }; // class ROL::TypeB::InteriorPointAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_KelleySachsAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_KelleySachsAlgorithm.hpp index 6af146dfcd2f..e2333e3a1b10 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_KelleySachsAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_KelleySachsAlgorithm.hpp @@ -130,7 +130,7 @@ class KelleySachsAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; private: void initialize(Vector &x, diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_LSecantBAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_LSecantBAlgorithm.hpp index 0787479278ab..a427b3ef6d8a 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_LSecantBAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_LSecantBAlgorithm.hpp @@ -112,7 +112,7 @@ class LSecantBAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; private: void initialize(Vector &x, diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_LinMoreAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_LinMoreAlgorithm.hpp index 044a96f73a86..d93895f85b07 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_LinMoreAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_LinMoreAlgorithm.hpp @@ -144,7 +144,7 @@ class LinMoreAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; private: void initialize(Vector &x, diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_MoreauYosidaAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_MoreauYosidaAlgorithm.hpp index b28574daae76..78de9edc1489 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_MoreauYosidaAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_MoreauYosidaAlgorithm.hpp @@ -108,7 +108,7 @@ class MoreauYosidaAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; }; // class ROL::TypeB::MoreauYosidaAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_NewtonKrylovAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_NewtonKrylovAlgorithm.hpp index 686c3e3b96c0..aa880c24a3ac 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_NewtonKrylovAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_NewtonKrylovAlgorithm.hpp @@ -195,7 +195,7 @@ class NewtonKrylovAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; }; // class ROL::TypeB::NewtonKrylovAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_PrimalDualActiveSetAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_PrimalDualActiveSetAlgorithm.hpp index 29708719af86..c648439f8f3c 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_PrimalDualActiveSetAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_PrimalDualActiveSetAlgorithm.hpp @@ -258,7 +258,7 @@ class PrimalDualActiveSetAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; }; // class ROL::TypeB::PrimalDualActiveSetAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_QuasiNewtonAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_QuasiNewtonAlgorithm.hpp index 884f5decb5ac..822e59e81f9d 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_QuasiNewtonAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_QuasiNewtonAlgorithm.hpp @@ -104,7 +104,7 @@ class QuasiNewtonAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream &os, bool write_header = false ) const override; + void writeOutput( std::ostream &os, const bool write_header = false ) const override; }; // class ROL::TypeB::QuasiNewtonAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_SpectralGradientAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_SpectralGradientAlgorithm.hpp index 605693ec808a..e04f481d9b58 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_SpectralGradientAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_SpectralGradientAlgorithm.hpp @@ -86,7 +86,7 @@ class SpectralGradientAlgorithm : public TypeB::Algorithm { void writeName( std::ostream &os ) const override; - void writeOutput( std::ostream &os, bool write_header = false ) const override; + void writeOutput( std::ostream &os, const bool write_header = false ) const override; }; // class ROL::TypeB::SpectralGradientAlgorithm diff --git a/packages/rol/src/algorithm/TypeB/ROL_TypeB_TrustRegionSPGAlgorithm.hpp b/packages/rol/src/algorithm/TypeB/ROL_TypeB_TrustRegionSPGAlgorithm.hpp index 0b4968f57c83..bcc0cac4bb69 100644 --- a/packages/rol/src/algorithm/TypeB/ROL_TypeB_TrustRegionSPGAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeB/ROL_TypeB_TrustRegionSPGAlgorithm.hpp @@ -140,7 +140,7 @@ class TrustRegionSPGAlgorithm : public TypeB::Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool write_header = false ) const override; + void writeOutput( std::ostream& os, const bool write_header = false ) const override; private: void initialize(Vector &x, diff --git a/packages/rol/src/algorithm/TypeE/ROL_TypeE_Algorithm.hpp b/packages/rol/src/algorithm/TypeE/ROL_TypeE_Algorithm.hpp index 18393662ace8..49fd832411d3 100644 --- a/packages/rol/src/algorithm/TypeE/ROL_TypeE_Algorithm.hpp +++ b/packages/rol/src/algorithm/TypeE/ROL_TypeE_Algorithm.hpp @@ -172,7 +172,7 @@ class Algorithm { /** \brief Print iterate status. */ - virtual void writeOutput( std::ostream& os, bool write_header = false ) const; + virtual void writeOutput( std::ostream& os, const bool write_header = false ) const; virtual void writeExitStatus( std::ostream& os ) const; diff --git a/packages/rol/src/algorithm/TypeG/ROL_TypeG_Algorithm.hpp b/packages/rol/src/algorithm/TypeG/ROL_TypeG_Algorithm.hpp index cc440a0a07f2..666cfe3bd60e 100644 --- a/packages/rol/src/algorithm/TypeG/ROL_TypeG_Algorithm.hpp +++ b/packages/rol/src/algorithm/TypeG/ROL_TypeG_Algorithm.hpp @@ -405,7 +405,7 @@ class Algorithm { /** \brief Print iterate status. */ - virtual void writeOutput( std::ostream& os, bool write_header = false ) const; + virtual void writeOutput( std::ostream& os, const bool write_header = false ) const; virtual void writeExitStatus( std::ostream& os ) const; diff --git a/packages/rol/src/algorithm/TypeG/augmentedlagrangian/ROL_AugmentedLagrangianObjective.hpp b/packages/rol/src/algorithm/TypeG/augmentedlagrangian/ROL_AugmentedLagrangianObjective.hpp index 19c352ab154b..9f077e3200c0 100644 --- a/packages/rol/src/algorithm/TypeG/augmentedlagrangian/ROL_AugmentedLagrangianObjective.hpp +++ b/packages/rol/src/algorithm/TypeG/augmentedlagrangian/ROL_AugmentedLagrangianObjective.hpp @@ -51,6 +51,7 @@ #include "ROL_Types.hpp" #include "ROL_Ptr.hpp" #include "ROL_ScalarController.hpp" +#include "ROL_ParameterList.hpp" #include /** @ingroup func_group diff --git a/packages/rol/src/algorithm/TypeG/fletcher/ROL_FletcherObjectiveBase_Def.hpp b/packages/rol/src/algorithm/TypeG/fletcher/ROL_FletcherObjectiveBase_Def.hpp index 3321ff35c179..be11ff750ac2 100644 --- a/packages/rol/src/algorithm/TypeG/fletcher/ROL_FletcherObjectiveBase_Def.hpp +++ b/packages/rol/src/algorithm/TypeG/fletcher/ROL_FletcherObjectiveBase_Def.hpp @@ -140,7 +140,7 @@ template Ptr> FletcherObjectiveBase::getGradient(const Vector& x) { // TODO: Figure out reasonable tolerance Real tol = static_cast(1e-12); - gradient(*xdual_, x, tol); + this->gradient(*xdual_, x, tol); return xdual_; } diff --git a/packages/rol/src/algorithm/TypeG/moreauyosida/ROL_MoreauYosidaObjective.hpp b/packages/rol/src/algorithm/TypeG/moreauyosida/ROL_MoreauYosidaObjective.hpp index ddf68a92ad55..c22642f3b831 100644 --- a/packages/rol/src/algorithm/TypeG/moreauyosida/ROL_MoreauYosidaObjective.hpp +++ b/packages/rol/src/algorithm/TypeG/moreauyosida/ROL_MoreauYosidaObjective.hpp @@ -50,6 +50,7 @@ #include "ROL_Types.hpp" #include "ROL_Ptr.hpp" #include "ROL_ScalarController.hpp" +#include "ROL_ParameterList.hpp" #include /** @ingroup func_group diff --git a/packages/rol/src/algorithm/TypeU/ROL_TypeU_Algorithm.hpp b/packages/rol/src/algorithm/TypeU/ROL_TypeU_Algorithm.hpp index f69954574e40..7240787ae02b 100644 --- a/packages/rol/src/algorithm/TypeU/ROL_TypeU_Algorithm.hpp +++ b/packages/rol/src/algorithm/TypeU/ROL_TypeU_Algorithm.hpp @@ -154,7 +154,7 @@ class Algorithm { /** \brief Print iterate status. */ - virtual void writeOutput( std::ostream& os, bool write_header = false ) const; + virtual void writeOutput( std::ostream& os, const bool write_header = false ) const; virtual void writeExitStatus( std::ostream& os ) const; diff --git a/packages/rol/src/algorithm/TypeU/ROL_TypeU_BundleAlgorithm.hpp b/packages/rol/src/algorithm/TypeU/ROL_TypeU_BundleAlgorithm.hpp index ad76539bd2b0..2fa88bd07f56 100644 --- a/packages/rol/src/algorithm/TypeU/ROL_TypeU_BundleAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeU/ROL_TypeU_BundleAlgorithm.hpp @@ -108,7 +108,7 @@ class BundleAlgorithm : public Algorithm { void writeName( std::ostream& os) const override; - void writeOutput( std::ostream& os, bool print_header = false ) const override; + void writeOutput( std::ostream& os, const bool print_header = false ) const override; private: diff --git a/packages/rol/src/algorithm/TypeU/ROL_TypeU_LineSearchAlgorithm.hpp b/packages/rol/src/algorithm/TypeU/ROL_TypeU_LineSearchAlgorithm.hpp index ed81db4a60f7..8e5cf1764963 100644 --- a/packages/rol/src/algorithm/TypeU/ROL_TypeU_LineSearchAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeU/ROL_TypeU_LineSearchAlgorithm.hpp @@ -153,7 +153,7 @@ class LineSearchAlgorithm : public Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool print_header = false ) const override; + void writeOutput( std::ostream& os, const bool print_header = false ) const override; }; // class ROL::TypeU::LineSearchAlgorithm diff --git a/packages/rol/src/algorithm/TypeU/ROL_TypeU_TrustRegionAlgorithm.hpp b/packages/rol/src/algorithm/TypeU/ROL_TypeU_TrustRegionAlgorithm.hpp index 252d1d971468..bf7bdba09bf8 100644 --- a/packages/rol/src/algorithm/TypeU/ROL_TypeU_TrustRegionAlgorithm.hpp +++ b/packages/rol/src/algorithm/TypeU/ROL_TypeU_TrustRegionAlgorithm.hpp @@ -118,7 +118,7 @@ class TrustRegionAlgorithm : public Algorithm { void writeName( std::ostream& os ) const override; - void writeOutput( std::ostream& os, bool print_header = false ) const override; + void writeOutput( std::ostream& os, const bool print_header = false ) const override; private: diff --git a/packages/rol/src/algorithm/TypeU/linesearch/descent/ROL_DescentDirection_U_Factory.hpp b/packages/rol/src/algorithm/TypeU/linesearch/descent/ROL_DescentDirection_U_Factory.hpp index d754edf8e450..9b1cb2f881e8 100644 --- a/packages/rol/src/algorithm/TypeU/linesearch/descent/ROL_DescentDirection_U_Factory.hpp +++ b/packages/rol/src/algorithm/TypeU/linesearch/descent/ROL_DescentDirection_U_Factory.hpp @@ -49,6 +49,7 @@ #include "ROL_NonlinearCG_U.hpp" #include "ROL_Newton_U.hpp" #include "ROL_NewtonKrylov_U.hpp" +#include "ROL_LineSearch_U_Types.hpp" namespace ROL { template diff --git a/packages/rol/src/elementwise/ROL_Elementwise_Function.hpp b/packages/rol/src/elementwise/ROL_Elementwise_Function.hpp index e385f5f6dc51..77f3f67110a4 100644 --- a/packages/rol/src/elementwise/ROL_Elementwise_Function.hpp +++ b/packages/rol/src/elementwise/ROL_Elementwise_Function.hpp @@ -41,6 +41,10 @@ // ************************************************************************ // @HEADER +#ifndef ROL_VECTOR_H +#include "ROL_Vector.hpp" +#else + #ifndef ROL_ELEMENTWISE_FUNCTION_H #define ROL_ELEMENTWISE_FUNCTION_H @@ -132,3 +136,4 @@ void applyBinaryInPlace(ROL::Vector &x, const ROL::Vector &y, Func f #include #endif +#endif diff --git a/packages/rol/src/elementwise/ROL_Elementwise_Reduce.hpp b/packages/rol/src/elementwise/ROL_Elementwise_Reduce.hpp index 29f28be408b0..ec91a146de0b 100644 --- a/packages/rol/src/elementwise/ROL_Elementwise_Reduce.hpp +++ b/packages/rol/src/elementwise/ROL_Elementwise_Reduce.hpp @@ -42,6 +42,10 @@ // ************************************************************************ // @HEADER +#ifndef ROL_VECTOR_H +#include "ROL_Vector.hpp" +#else + #ifndef ROL_ELEMENTWISE_REDUCE_H #define ROL_ELEMENTWISE_REDUCE_H @@ -190,3 +194,4 @@ class EuclideanNormSquared : public ReductionOp { } // namespace ROL #endif +#endif diff --git a/packages/rol/src/function/ROL_ConstraintAssembler.hpp b/packages/rol/src/function/ROL_ConstraintAssembler.hpp index 6d13791df9aa..0165ee0e3792 100644 --- a/packages/rol/src/function/ROL_ConstraintAssembler.hpp +++ b/packages/rol/src/function/ROL_ConstraintAssembler.hpp @@ -59,10 +59,10 @@ namespace ROL { template struct ConstraintData { - const Ptr> constraint; - const Ptr> multiplier; - const Ptr> residual; - const Ptr> bounds; + Ptr> constraint; + Ptr> multiplier; + Ptr> residual; + Ptr> bounds; ConstraintData(const Ptr> &con, const Ptr> &mul, diff --git a/packages/rol/src/function/ROL_UpdateType.hpp b/packages/rol/src/function/ROL_UpdateType.hpp index 8cedf2659681..308ecf6d7a3a 100644 --- a/packages/rol/src/function/ROL_UpdateType.hpp +++ b/packages/rol/src/function/ROL_UpdateType.hpp @@ -45,6 +45,7 @@ #define ROL_UPDATE_TYPE_HPP #include +#include namespace ROL { diff --git a/packages/rol/src/function/ROL_VectorController.hpp b/packages/rol/src/function/ROL_VectorController.hpp index ee191a41d035..fa05977aee26 100644 --- a/packages/rol/src/function/ROL_VectorController.hpp +++ b/packages/rol/src/function/ROL_VectorController.hpp @@ -47,6 +47,8 @@ #include "ROL_Vector.hpp" #include "ROL_UpdateType.hpp" +#include +#include namespace ROL { diff --git a/packages/rol/src/function/boundconstraint/ROL_BoundConstraint_Partitioned.hpp b/packages/rol/src/function/boundconstraint/ROL_BoundConstraint_Partitioned.hpp index 5e1aa06ba0cf..20666576bdde 100644 --- a/packages/rol/src/function/boundconstraint/ROL_BoundConstraint_Partitioned.hpp +++ b/packages/rol/src/function/boundconstraint/ROL_BoundConstraint_Partitioned.hpp @@ -141,12 +141,6 @@ class BoundConstraint_Partitioned : public BoundConstraint { } void update( const Vector &x, bool flag = true, int iter = -1 ) { - const PV &xpv = dynamic_cast(x); - for( uint k=0; kisActivated() ) { - bnd_[k]->update(*(xpv.get(k)),flag,iter); - } - } } void project( Vector &x ) { diff --git a/packages/rol/src/function/dynamic/ROL_DynamicConstraint.hpp b/packages/rol/src/function/dynamic/ROL_DynamicConstraint.hpp index de3629942f0c..8606d677ca89 100644 --- a/packages/rol/src/function/dynamic/ROL_DynamicConstraint.hpp +++ b/packages/rol/src/function/dynamic/ROL_DynamicConstraint.hpp @@ -46,27 +46,35 @@ #ifndef ROL_DYNAMICCONSTRAINT_HPP #define ROL_DYNAMICCONSTRAINT_HPP -#include "ROL_DynamicFunction.hpp" -#include "ROL_NonlinearLeastSquaresObjective_Dynamic.hpp" -#include "ROL_Constraint_DynamicState.hpp" #include "ROL_Objective_FSsolver.hpp" -#include "ROL_Algorithm.hpp" -#include "ROL_CompositeStep.hpp" -#include "ROL_TrustRegionStep.hpp" -#include "ROL_ConstraintStatusTest.hpp" #include "ROL_Types.hpp" +#include "ROL_DynamicFunction.hpp" +#include "ROL_TypeU_TrustRegionAlgorithm.hpp" +#include "ROL_TypeE_AugmentedLagrangianAlgorithm.hpp" + +namespace ROL { + +template +class DynamicConstraint; + +template +class Constraint_DynamicState; +} + +#include "ROL_Constraint_DynamicState.hpp" +#include "ROL_NonlinearLeastSquaresObjective_Dynamic.hpp" /** @ingroup dynamic_group \class ROL::DynamicConstraint - \brief Defines the time-dependent constraint operator interface for + \brief Defines the time-dependent constraint operator interface for simulation-based optimization. - This constraint interface inherits from ROL_Constraint_SimOpt. Though + This constraint interface inherits from ROL_Constraint_SimOpt. Though the interface takes two simulation space vectors from spaces - \f$\mathcal{U_o}\times\mathcal{U_n}\f$. The space \f$\mathcal{U_o}\f$ is - ``old'' information that accounts for the initial condition on the time - interval. The space \f$\mathcal{U_n}\f$ is the ``new'' variables that can + \f$\mathcal{U_o}\times\mathcal{U_n}\f$. The space \f$\mathcal{U_o}\f$ is + ``old'' information that accounts for the initial condition on the time + interval. The space \f$\mathcal{U_n}\f$ is the ``new'' variables that can be determined by satisfying constraints in the form \f[ c(u_o,u_n,z,t_o,t_n) = 0 \,. @@ -74,7 +82,6 @@ */ - namespace ROL { template @@ -142,7 +149,7 @@ class DynamicConstraint : public DynamicFunction { firstSolve_ ( true ) {} - + virtual void update( const V& uo, const V& un, const V& z, const TS& ts ) { update_uo( uo, ts ); update_un( un, ts ); @@ -153,14 +160,13 @@ class DynamicConstraint : public DynamicFunction { using DynamicFunction::update_un; using DynamicFunction::update_z; - virtual void value( V& c, const V& uo, const V& un, + virtual void value( V& c, const V& uo, const V& un, const V& z, const TS& ts ) const = 0; - virtual void solve( V& c, const V& uo, V& un, + virtual void solve( V& c, const V& uo, V& un, const V& z, const TS& ts ) { - if ( zero_ ) { - un.zero(); - } + if ( zero_ ) un.zero(); + Ptr stream = makeStreamPtr(std::cout, print_); update(uo,un,z,ts); value(c,uo,un,z,ts); Real cnorm = c.norm(); @@ -174,12 +180,12 @@ class DynamicConstraint : public DynamicFunction { Real alpha(1), tmp(0); int cnt = 0; if ( print_ ) { - std::cout << "\n Default DynamicConstraint::solve\n"; - std::cout << " "; - std::cout << std::setw(6) << std::left << "iter"; - std::cout << std::setw(15) << std::left << "rnorm"; - std::cout << std::setw(15) << std::left << "alpha"; - std::cout << "\n"; + *stream << " Default DynamicConstraint::solve" << std::endl; + *stream << " "; + *stream << std::setw(6) << std::left << "iter"; + *stream << std::setw(15) << std::left << "rnorm"; + *stream << std::setw(15) << std::left << "alpha"; + *stream << std::endl; } for (cnt = 0; cnt < maxit_; ++cnt) { // Compute Newton step @@ -202,13 +208,13 @@ class DynamicConstraint : public DynamicFunction { tmp = c.norm(); } if ( print_ ) { - std::cout << " "; - std::cout << std::setw(6) << std::left << cnt; - std::cout << std::scientific << std::setprecision(6); - std::cout << std::setw(15) << std::left << tmp; - std::cout << std::scientific << std::setprecision(6); - std::cout << std::setw(15) << std::left << alpha; - std::cout << "\n"; + *stream << " "; + *stream << std::setw(6) << std::left << cnt; + *stream << std::scientific << std::setprecision(6); + *stream << std::setw(15) << std::left << tmp; + *stream << std::scientific << std::setprecision(6); + *stream << std::setw(15) << std::left << alpha; + *stream << std::endl; } // Update iterate cnorm = tmp; @@ -221,42 +227,32 @@ class DynamicConstraint : public DynamicFunction { } } if (solverType_==1 || (solverType_==3 && cnorm > ctol)) { - Ptr> con_ptr = makePtrFromRef(*this); - Ptr> uo_ptr = makePtrFromRef(uo); - Ptr> z_ptr = makePtrFromRef(z); - Ptr> ts_ptr = makePtrFromRef(ts); - Ptr> obj_ptr - = makePtr>(con_ptr,c,uo_ptr,z_ptr,ts_ptr,true); + Ptr> con = makePtrFromRef(*this); + Ptr> obj + = makePtr>(con,c,makePtrFromRef(uo),makePtrFromRef(z),makePtrFromRef(ts),true); ParameterList parlist; + parlist.sublist("General").set("Output Level",1); + parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100); + parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG"); parlist.sublist("Status Test").set("Gradient Tolerance",ctol); parlist.sublist("Status Test").set("Step Tolerance",stol_); parlist.sublist("Status Test").set("Iteration Limit",maxit_); - parlist.sublist("Step").sublist("Trust Region").set("Subproblem Solver","Truncated CG"); - parlist.sublist("General").sublist("Krylov").set("Iteration Limit",100); - Ptr> step = makePtr>(parlist); - Ptr> status = makePtr>(parlist); - Ptr> algo = makePtr>(step,status,false); - algo->run(un,*obj_ptr,print_); + Ptr> algo = makePtr>(parlist); + algo->run(un,*obj,*stream); value(c,uo,un,z,ts); } if (solverType_==2 || (solverType_==4 && cnorm > ctol)) { - Ptr> con_ptr = makePtrFromRef(*this); - Ptr> uo_ptr = makePtrFromRef(uo); - Ptr> z_ptr = makePtrFromRef(z); - Ptr> ts_ptr = makePtrFromRef(ts); - Ptr> cn_ptr - = makePtr>(con_ptr,uo_ptr,z_ptr,ts_ptr); - Ptr> obj_ptr - = makePtr>(); + Ptr> con + = makePtr>(makePtrFromRef(*this),makePtrFromRef(uo),makePtrFromRef(z),makePtrFromRef(ts)); + Ptr> obj = makePtr>(); + Ptr> l = c.dual().clone(); ParameterList parlist; + parlist.sublist("General").set("Output Level",1); parlist.sublist("Status Test").set("Constraint Tolerance",ctol); parlist.sublist("Status Test").set("Step Tolerance",stol_); parlist.sublist("Status Test").set("Iteration Limit",maxit_); - Ptr> step = makePtr>(parlist); - Ptr> status = makePtr>(parlist); - Ptr> algo = makePtr>(step,status,false); - Ptr> l = c.dual().clone(); - algo->run(un,*l,*obj_ptr,*cn_ptr,print_); + Ptr> algo = makePtr>(parlist); + algo->run(un,*obj,*con,*l,*stream); value(c,uo,un,z,ts); } if (solverType_ > 4 || solverType_ < 0) { @@ -300,94 +296,94 @@ class DynamicConstraint : public DynamicFunction { //---------------------------------------------------------------------------- // Partial Jacobians - virtual void applyJacobian_uo( V& jv, const V& vo, const V& uo, - const V& un, const V& z, + virtual void applyJacobian_uo( V& jv, const V& vo, const V& uo, + const V& un, const V& z, const TS& ts ) const {} - virtual void applyJacobian_un( V& jv, const V& vn, const V& uo, - const V& un, const V& z, + virtual void applyJacobian_un( V& jv, const V& vn, const V& uo, + const V& un, const V& z, const TS& ts ) const {} - virtual void applyJacobian_z( V& jv, const V& vz, const V& uo, - const V& un, const V& z, + virtual void applyJacobian_z( V& jv, const V& vz, const V& uo, + const V& un, const V& z, const TS& ts ) const {} //---------------------------------------------------------------------------- // Adjoint partial Jacobians - virtual void applyAdjointJacobian_uo( V& ajv, const V& dualv, const V& uo, - const V& un, const V& z, + virtual void applyAdjointJacobian_uo( V& ajv, const V& dualv, const V& uo, + const V& un, const V& z, const TS& ts ) const {} - virtual void applyAdjointJacobian_un( V& ajv, const V& dualv, const V& uo, - const V& un, const V& z, + virtual void applyAdjointJacobian_un( V& ajv, const V& dualv, const V& uo, + const V& un, const V& z, const TS& ts ) const {} - virtual void applyAdjointJacobian_z( V& ajv, const V& dualv, const V& uo, - const V& un, const V& z, + virtual void applyAdjointJacobian_z( V& ajv, const V& dualv, const V& uo, + const V& un, const V& z, const TS& ts ) const {} //---------------------------------------------------------------------------- // Inverses - virtual void applyInverseJacobian_un( V& ijv, const V& vn, const V& uo, - const V& un, const V& z, + virtual void applyInverseJacobian_un( V& ijv, const V& vn, const V& uo, + const V& un, const V& z, const TS& ts ) const {} - - virtual void applyInverseAdjointJacobian_un( V& iajv, const V& vn, const V& uo, - const V& un, const V& z, + + virtual void applyInverseAdjointJacobian_un( V& iajv, const V& vn, const V& uo, + const V& un, const V& z, const TS& ts ) const {} //---------------------------------------------------------------------------- // Adjoint Hessian components virtual void applyAdjointHessian_un_un( V& ahwv, const V& wn, const V& vn, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_un_uo( V& ahwv, const V& w, const V& vn, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_un_z( V& ahwv, const V& w, const V& vn, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_uo_un( V& ahwv, const V& w, const V& vo, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_uo_uo( V& ahwv, const V& w, const V& v, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_uo_z( V& ahwv, const V& w, const V& vo, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_z_un( V& ahwv, const V& w, const V& vz, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } virtual void applyAdjointHessian_z_uo( V& ahwv, const V& w, const V& vz, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } - + virtual void applyAdjointHessian_z_z( V& ahwv, const V& w, const V& vz, - const V& uo, const V& un, + const V& uo, const V& un, const V& z, const TS& ts ) const { ahwv.zero(); } diff --git a/packages/rol/src/function/dynamic/ROL_ReducedDynamicObjective.hpp b/packages/rol/src/function/dynamic/ROL_ReducedDynamicObjective.hpp index 73ccfb35a2ee..13a286a1f35e 100644 --- a/packages/rol/src/function/dynamic/ROL_ReducedDynamicObjective.hpp +++ b/packages/rol/src/function/dynamic/ROL_ReducedDynamicObjective.hpp @@ -251,6 +251,10 @@ class ReducedDynamicObjective : public Objective { return ROL::PartitionedVector::create(x, Nt_); } + size_type getStateRank() const { return rankState_; } + size_type getAdjointRank() const { return rankAdjoint_; } + size_type getStateSensitivityRank() const { return rankStateSens_; } + void update( const Vector &x, bool flag = true, int iter = -1 ) { if (useSketch_) { stateSketch_->reset(true); diff --git a/packages/rol/src/function/objective/ROL_AffineTransformObjective_Def.hpp b/packages/rol/src/function/objective/ROL_AffineTransformObjective_Def.hpp index c3c8f3ba0fb7..d25515185191 100644 --- a/packages/rol/src/function/objective/ROL_AffineTransformObjective_Def.hpp +++ b/packages/rol/src/function/objective/ROL_AffineTransformObjective_Def.hpp @@ -63,8 +63,8 @@ AffineTransformObjective::AffineTransformObjective(const Ptr> &acon, const Ptr> &storage) : obj_(obj), acon_(acon), storage_(storage) { - primal_ = acon_->createRangeSpaceVector(); - Av_ = acon_->createRangeSpaceVector(); + primal_ = acon->createRangeSpaceVector(); + Av_ = acon->createRangeSpaceVector(); dual_ = primal_->dual().clone(); if (storage == nullPtr) storage_ = makePtr>(); } diff --git a/packages/rol/src/function/objective/ROL_QuadraticObjective.hpp b/packages/rol/src/function/objective/ROL_QuadraticObjective.hpp index 0c924d3562d0..382c7d3e1e35 100644 --- a/packages/rol/src/function/objective/ROL_QuadraticObjective.hpp +++ b/packages/rol/src/function/objective/ROL_QuadraticObjective.hpp @@ -47,6 +47,7 @@ #include "ROL_Objective.hpp" #include "ROL_Vector.hpp" #include "ROL_Ptr.hpp" +#include "ROL_LinearOperator.hpp" /** @ingroup func_group \class ROL::QuadraticObjective diff --git a/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint.hpp b/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint.hpp index 29d887b598d6..b82a5670cbd6 100644 --- a/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint.hpp +++ b/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint.hpp @@ -110,7 +110,7 @@ class BilinearConstraint : public Constraint_SimOpt, public ProfiledClass< void getFactor(Vector &F, const std::vector ¶m) const; Real getNoise(int k) const; void getTraceSample(Vector &g, const std::vector ¶m) const; - void sumAll(Real *in, Real *out, int size) const; + // void sumAll(Real *in, Real *out, int size) const; Real logDeterminant(const Vector &z); }; // class BilinearConstraint diff --git a/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint_Def.hpp b/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint_Def.hpp index 84b3aa84114e..bb25ac4961b6 100644 --- a/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint_Def.hpp +++ b/packages/rol/src/oed/constraint/ROL_OED_BilinearConstraint_Def.hpp @@ -269,10 +269,10 @@ void BilinearConstraint::getTraceSample(Vector &g, const std::vector traceSampler_->get(g,param); } -template -void BilinearConstraint::sumAll(Real *in, Real *out, int size) const { - factors_->sumAll(in,out,size); -} +// template +// void BilinearConstraint::sumAll(Real *in, Real *out, int size) const { +// factors_->sumAll(in,out,size); +// } template Real BilinearConstraint::logDeterminant(const Vector &z) { diff --git a/packages/rol/src/oed/factors/ROL_OED_Factors.hpp b/packages/rol/src/oed/factors/ROL_OED_Factors.hpp index 3b53e77d68b2..234022196460 100644 --- a/packages/rol/src/oed/factors/ROL_OED_Factors.hpp +++ b/packages/rol/src/oed/factors/ROL_OED_Factors.hpp @@ -115,7 +115,7 @@ class Factors : public ProfiledClass { // Compute F(param)^T c void evaluate(Vector &F, const std::vector ¶m) const; - void sumAll(Real *in, Real *out, int size) const; + // void sumAll(Real *in, Real *out, int size) const; int numFactors() const; diff --git a/packages/rol/src/sol/function/ROL_AbsoluteValue.hpp b/packages/rol/src/sol/function/ROL_AbsoluteValue.hpp index 1e78534f0f63..c49a2d4c9625 100644 --- a/packages/rol/src/sol/function/ROL_AbsoluteValue.hpp +++ b/packages/rol/src/sol/function/ROL_AbsoluteValue.hpp @@ -46,6 +46,7 @@ #include "ROL_Types.hpp" #include "ROL_PositiveFunction.hpp" +#include "ROL_ParameterList.hpp" namespace ROL { diff --git a/packages/rol/src/sol/function/ROL_RiskBoundConstraint.hpp b/packages/rol/src/sol/function/ROL_RiskBoundConstraint.hpp index f79bd504db07..c3e11507212d 100644 --- a/packages/rol/src/sol/function/ROL_RiskBoundConstraint.hpp +++ b/packages/rol/src/sol/function/ROL_RiskBoundConstraint.hpp @@ -233,26 +233,6 @@ class RiskBoundConstraint : public BoundConstraint { } } - void update( const Vector &x, bool flag = true, int iter = -1 ) { - if ( augmentedObj_ && activatedObj_ ) { - Ptr> xs = dynamic_cast&>(x).getStatisticVector(0); - statObj_bc_->update(*xs,flag,iter); - } - if (augmentedCon_) { - int size = statCon_bc_.size(); - for (int i = 0; i < size; ++i) { - if (activatedCon_[i]) { - Ptr> xs = dynamic_cast&>(x).getStatisticVector(1,i); - statCon_bc_[i]->update(*xs,flag,iter); - } - } - } - if ( bc_ != nullPtr && bc_->isActivated() ) { - Ptr> xv = dynamic_cast&>(x).getVector(); - bc_->update(*xv,flag,iter); - } - } - void project( Vector &x ) { if ( augmentedObj_ && activatedObj_ ) { Ptr> xs = dynamic_cast&>(x).getStatisticVector(0); diff --git a/packages/rol/src/sol/function/randvarfunctional/ROL_RandVarFunctional.hpp b/packages/rol/src/sol/function/randvarfunctional/ROL_RandVarFunctional.hpp index 5b2658db762c..0b15a580e7cc 100644 --- a/packages/rol/src/sol/function/randvarfunctional/ROL_RandVarFunctional.hpp +++ b/packages/rol/src/sol/function/randvarfunctional/ROL_RandVarFunctional.hpp @@ -45,12 +45,13 @@ #define ROL_RANDVARFUNCTIONAL_HPP #include "ROL_Vector.hpp" +#include "ROL_Objective.hpp" #include "ROL_Ptr.hpp" #include "ROL_SampleGenerator.hpp" #include "ROL_ScalarController.hpp" #include "ROL_VectorController.hpp" -/** @ingroup stochastic_group +/** @ingroup stochastic_group \class ROL::RandVarFunctional \brief Provides the interface to implement any functional that maps a random variable to a (extended) real number. diff --git a/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_Chi2Divergence.hpp b/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_Chi2Divergence.hpp index c03c2bb9d065..a990af40a40c 100644 --- a/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_Chi2Divergence.hpp +++ b/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_Chi2Divergence.hpp @@ -89,7 +89,7 @@ class Chi2Divergence : public FDivergence { */ Chi2Divergence(ROL::ParameterList &parlist) : FDivergence(parlist) {} - Real Fprimal(Real x, int deriv = 0) { + Real Fprimal(Real x, int deriv = 0) const { Real zero(0), one(1), half(0.5), val(0); if (deriv==0) { val = (x < zero) ? ROL_INF() : half*(x-one)*(x-one); @@ -107,7 +107,7 @@ class Chi2Divergence : public FDivergence { return val; } - Real Fdual(Real x, int deriv = 0) { + Real Fdual(Real x, int deriv = 0) const { Real zero(0), one(1), half(0.5), val(0); if (deriv==0) { val = (x < -one) ? -half : (half*x + one)*x; diff --git a/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_FDivergence.hpp b/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_FDivergence.hpp index e58d9c44c0f4..19c62ad8fa53 100644 --- a/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_FDivergence.hpp +++ b/packages/rol/src/sol/function/randvarfunctional/risk/fdivergence/ROL_FDivergence.hpp @@ -146,7 +146,7 @@ class FDivergence : public RandVarFunctional { Upon return, Fprimal returns \f$F(x)\f$ or a derivative of \f$F(x)\f$. */ - virtual Real Fprimal(Real x, int deriv = 0) = 0; + virtual Real Fprimal(Real x, int deriv = 0) const = 0; /** \brief Implementation of the scalar dual F function. @@ -160,7 +160,7 @@ class FDivergence : public RandVarFunctional { F^*(y) = \sup_{x\in\mathbb{R}}\{xy - F(x)\}. \f] */ - virtual Real Fdual(Real x, int deriv = 0) = 0; + virtual Real Fdual(Real x, int deriv = 0) const = 0; bool check(std::ostream &outStream = std::cout) const { const Real tol(std::sqrt(ROL_EPSILON())); diff --git a/packages/rol/src/sol/function/randvarfunctional/risk/spectral/ROL_SpectralRisk.hpp b/packages/rol/src/sol/function/randvarfunctional/risk/spectral/ROL_SpectralRisk.hpp index cb563c70849a..97e9ae1597fb 100644 --- a/packages/rol/src/sol/function/randvarfunctional/risk/spectral/ROL_SpectralRisk.hpp +++ b/packages/rol/src/sol/function/randvarfunctional/risk/spectral/ROL_SpectralRisk.hpp @@ -94,7 +94,7 @@ class SpectralRisk : public RandVarFunctional { std::vector wts_; std::vector pts_; - void checkInputs(ROL::Ptr > &dist = ROL::nullPtr) const { + void checkInputs(const ROL::Ptr > &dist = ROL::nullPtr) const { ROL_TEST_FOR_EXCEPTION(plusFunction_ == ROL::nullPtr, std::invalid_argument, ">>> ERROR (ROL::SpectralRisk): PlusFunction pointer is null!"); if ( dist != ROL::nullPtr) { diff --git a/packages/rol/src/step/secant/ROL_BarzilaiBorwein.hpp b/packages/rol/src/step/secant/ROL_BarzilaiBorwein.hpp index 28330bb51f56..5759070a1971 100644 --- a/packages/rol/src/step/secant/ROL_BarzilaiBorwein.hpp +++ b/packages/rol/src/step/secant/ROL_BarzilaiBorwein.hpp @@ -41,6 +41,11 @@ // ************************************************************************ // @HEADER + +#ifndef ROL_SECANTFACTORY_H +#include "ROL_SecantFactory.hpp" +#else + #ifndef ROL_BARZILAIBORWEIN_H #define ROL_BARZILAIBORWEIN_H @@ -101,3 +106,4 @@ class BarzilaiBorwein : public Secant { } #endif +#endif diff --git a/packages/rol/src/step/trustregion/ROL_CauchyPoint.hpp b/packages/rol/src/step/trustregion/ROL_CauchyPoint.hpp index 41e70b08cc35..d7939c2a76a7 100644 --- a/packages/rol/src/step/trustregion/ROL_CauchyPoint.hpp +++ b/packages/rol/src/step/trustregion/ROL_CauchyPoint.hpp @@ -41,6 +41,10 @@ // ************************************************************************ // @HEADER +#ifndef ROL_TRUSTREGIONFACTORY_H +#include "ROL_TrustRegionFactory.hpp" +#else + #ifndef ROL_CAUCHYPOINT_H #define ROL_CAUCHYPOINT_H @@ -277,3 +281,4 @@ class CauchyPoint : public TrustRegion { } #endif +#endif diff --git a/packages/rol/src/step/trustregion/ROL_DogLeg.hpp b/packages/rol/src/step/trustregion/ROL_DogLeg.hpp index 792bc65d241a..86ad2c8c6c6a 100644 --- a/packages/rol/src/step/trustregion/ROL_DogLeg.hpp +++ b/packages/rol/src/step/trustregion/ROL_DogLeg.hpp @@ -48,7 +48,7 @@ \brief Provides interface for dog leg trust-region subproblem solver. */ -#include "ROL_TrustRegion.hpp" +#include "ROL_CauchyPoint.hpp" #include "ROL_Types.hpp" namespace ROL { diff --git a/packages/rol/src/step/trustregion/ROL_TrustRegionTypes.hpp b/packages/rol/src/step/trustregion/ROL_TrustRegionTypes.hpp index b2355eaae556..3cb8ddf4c737 100644 --- a/packages/rol/src/step/trustregion/ROL_TrustRegionTypes.hpp +++ b/packages/rol/src/step/trustregion/ROL_TrustRegionTypes.hpp @@ -49,6 +49,8 @@ #ifndef ROL_TRUSTREGIONTYPES_HPP #define ROL_TRUSTREGIONTYPES_HPP +#include "ROL_Types.hpp" + namespace ROL { /** \enum ROL::ETrustRegion diff --git a/packages/rol/src/vector/ROL_Vector.hpp b/packages/rol/src/vector/ROL_Vector.hpp index 9b279039a6b0..c586aa9e0ceb 100644 --- a/packages/rol/src/vector/ROL_Vector.hpp +++ b/packages/rol/src/vector/ROL_Vector.hpp @@ -242,14 +242,14 @@ class Vector virtual void applyUnary( const Elementwise::UnaryFunction &f ) { ROL_UNUSED(f); ROL_TEST_FOR_EXCEPTION( true, std::logic_error, - "The method applyUnary wass called, but not implemented" << std::endl); + "The method applyUnary was called, but not implemented" << std::endl); } virtual void applyBinary( const Elementwise::BinaryFunction &f, const Vector &x ) { ROL_UNUSED(f); ROL_UNUSED(x); ROL_TEST_FOR_EXCEPTION( true, std::logic_error, - "The method applyBinary wass called, but not implemented" << std::endl); + "The method applyBinary was called, but not implemented" << std::endl); } virtual Real reduce( const Elementwise::ReductionOp &r ) const { diff --git a/packages/rol/src/zoo/ROL_ScalarTraits.hpp b/packages/rol/src/zoo/ROL_ScalarTraits.hpp index 9248f6a90697..affc3d4d412c 100644 --- a/packages/rol/src/zoo/ROL_ScalarTraits.hpp +++ b/packages/rol/src/zoo/ROL_ScalarTraits.hpp @@ -46,6 +46,7 @@ #define ROL_SCALARTRAITS_HPP #include +#include namespace ROL {