MOOCHO/Thyra Adapter Software Version of the Day
NLPInterfacePack_NLPDirectThyraModelEvaluator.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "NLPInterfacePack_NLPDirectThyraModelEvaluator.hpp"
00030 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00031 #include "AbstractLinAlgPack_VectorOut.hpp"
00032 #include "AbstractLinAlgPack_ThyraAccessors.hpp"
00033 #include "AbstractLinAlgPack_VectorSpaceThyra.hpp"
00034 #include "AbstractLinAlgPack_VectorMutableThyra.hpp"
00035 #include "AbstractLinAlgPack_MultiVectorMutableThyra.hpp"
00036 #include "AbstractLinAlgPack_MatrixOpNonsingThyra.hpp"
00037 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00038 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00039 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00040 #include "AbstractLinAlgPack_BasisSystem.hpp"
00041 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00042 #include "Thyra_ModelEvaluatorHelpers.hpp"
00043 #include "Thyra_DetachedVectorView.hpp"
00044 #include "Thyra_VectorStdOps.hpp"
00045 #include "Thyra_MultiVectorStdOps.hpp"
00046 #include "Teuchos_AbstractFactoryStd.hpp"
00047 #include "Teuchos_TestForException.hpp"
00048 #include "Teuchos_dyn_cast.hpp"
00049 
00050 namespace NLPInterfacePack {
00051 
00052 NLPDirectThyraModelEvaluator::NLPDirectThyraModelEvaluator()
00053   :DfDp_is_const_(false)
00054 {}
00055 
00056 NLPDirectThyraModelEvaluator::NLPDirectThyraModelEvaluator(
00057   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> > &model 
00058   ,const int p_idx 
00059   ,const int g_idx 
00060   ,const objDirecFiniteDiffCalculator_ptr_t objDirecFiniteDiffCalculator
00061   ,const conDirecFiniteDiffCalculator_ptr_t conDirecFiniteDiffCalculator
00062   )
00063   :DfDp_is_const_(false)
00064 {
00065   initialize(
00066     model,p_idx,g_idx
00067     ,objDirecFiniteDiffCalculator,conDirecFiniteDiffCalculator
00068     );
00069 }
00070 
00071 void NLPDirectThyraModelEvaluator::initialize(
00072   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> > &model
00073   ,const int p_idx
00074   ,const int g_idx
00075   ,const objDirecFiniteDiffCalculator_ptr_t objDirecFiniteDiffCalculator
00076   ,const conDirecFiniteDiffCalculator_ptr_t conDirecFiniteDiffCalculator
00077   )
00078 {
00079   typedef Thyra::ModelEvaluatorBase MEB;
00080   if(objDirecFiniteDiffCalculator.get())
00081     this->set_objDirecFiniteDiffCalculator(objDirecFiniteDiffCalculator);
00082   if(conDirecFiniteDiffCalculator.get())
00083     this->set_conDirecFiniteDiffCalculator(conDirecFiniteDiffCalculator);
00084   initializeBase(model,p_idx,g_idx);
00085   Thyra::ModelEvaluatorBase::OutArgs<double> model_outArgs = model->createOutArgs();
00086   MEB::DerivativeProperties model_W_properties = model_outArgs.get_W_properties();
00087   if( p_idx >= 0 ) {
00088     TEST_FOR_EXCEPTION(
00089       (conDirecFiniteDiffCalculator_.get()==0 && !model_outArgs.supports(MEB::OUT_ARG_DfDp,p_idx).supports(MEB::DERIV_MV_BY_COL))
00090       ,std::invalid_argument
00091       ,"Error, model must support computing DfDp("<<p_idx<<") as a"
00092       " column-oriented multi-vector if not using finite differences!"
00093       );
00094   }
00095   DfDp_is_const_  = false;
00096 }
00097 
00098 // Overridden public members from NLP
00099 
00100 void NLPDirectThyraModelEvaluator::initialize(bool test_setup)
00101 {
00102   if(initialized_) {
00103     NLPDirect::initialize(test_setup);
00104     return;
00105   }
00106   NLPThyraModelEvaluatorBase::initialize(test_setup);
00107   NLPDirect::initialize(test_setup);
00108 }
00109 
00110 void NLPDirectThyraModelEvaluator::unset_quantities()
00111 {
00112   NLPDirect::unset_quantities();
00113 }
00114 
00115 // Overridden public members from NLPObjGrad
00116 
00117 bool NLPDirectThyraModelEvaluator::supports_Gf() const
00118 {
00119   if(objDirecFiniteDiffCalculator_.get())
00120     return false;
00121   return true;
00122   // ToDo: Change this to false when only the operator for model_DgDx is
00123   // supported!
00124 }
00125 
00126 bool NLPDirectThyraModelEvaluator::supports_Gf_prod() const
00127 {
00128   return ( objDirecFiniteDiffCalculator_.get() != NULL );
00129 }
00130 
00131 value_type NLPDirectThyraModelEvaluator::calc_Gf_prod(
00132   const Vector& x, const Vector& d, bool newx
00133   ) const
00134 {
00135   TEST_FOR_EXCEPT(objDirecFiniteDiffCalculator_.get()==0);
00136   typedef Thyra::ModelEvaluatorBase MEB;
00137   MEB::InArgs<value_type> basePoint = model_->createInArgs();
00138   MEB::InArgs<value_type> directions = model_->createInArgs();
00139   MEB::OutArgs<value_type> baseFunc = model_->createOutArgs();
00140   MEB::OutArgs<value_type> variations = model_->createOutArgs();
00141   set_x(x,&basePoint);
00142   set_x(d,&directions);
00143   variations.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00144   objDirecFiniteDiffCalculator_->calcVariations(
00145     *model_,basePoint,directions,baseFunc,variations
00146     );
00147   return Thyra::get_ele(*variations.get_g(g_idx_),0);
00148 }
00149 
00150 // Overridden public members from NLPDirect
00151 
00152 Range1D NLPDirectThyraModelEvaluator::var_dep() const
00153 {
00154   return basis_sys_->var_dep();
00155 }
00156 
00157 Range1D NLPDirectThyraModelEvaluator::var_indep() const
00158 {
00159   return basis_sys_->var_indep();
00160 }
00161 
00162 const NLPDirect::mat_fcty_ptr_t
00163 NLPDirectThyraModelEvaluator::factory_D() const
00164 {
00165   return basis_sys_->factory_D();
00166 }
00167 
00168 const NLPDirect::mat_sym_nonsing_fcty_ptr_t
00169 NLPDirectThyraModelEvaluator::factory_S() const
00170 {
00171   return basis_sys_->factory_S();
00172 }
00173 
00174 void NLPDirectThyraModelEvaluator::calc_point(
00175   const Vector &x
00176   ,value_type *f
00177   ,VectorMutable *c
00178   ,bool recalc_c
00179   ,VectorMutable *Gf
00180   ,VectorMutable *py
00181   ,VectorMutable *rGf
00182   ,MatrixOp *GcU
00183   ,MatrixOp *D
00184   ,MatrixOp *Uz
00185   ) const
00186 {
00187   using Teuchos::FancyOStream;
00188   using Teuchos::OSTab;
00189   using Teuchos::dyn_cast;
00190   using Teuchos::RCP;
00191   using Teuchos::rcp;
00192   using Teuchos::rcp_const_cast;
00193   using Teuchos::rcp_dynamic_cast;
00194   using AbstractLinAlgPack::VectorSpaceThyra;
00195   using AbstractLinAlgPack::VectorMutableThyra;
00196   using AbstractLinAlgPack::MultiVectorMutableThyra;
00197   using AbstractLinAlgPack::MatrixOpThyra;
00198   using AbstractLinAlgPack::MatrixOpNonsingThyra;
00199   typedef Thyra::ModelEvaluatorBase MEB;
00200   typedef Teuchos::RCP<Thyra::VectorBase<value_type> > TVecPtr;
00201   typedef MEB::DerivativeMultiVector<value_type> DerivMV;
00202   typedef MEB::Derivative<value_type> Deriv;
00203   //
00204   // Get output and verbosity
00205   //
00206   const Teuchos::RCP<Teuchos::FancyOStream>
00207     out = this->getOStream();
00208   const Teuchos::EVerbosityLevel
00209     verbLevel = ( showModelEvaluatorTrace() ? this->getVerbLevel() : Teuchos::VERB_NONE );
00210   Teuchos::OSTab tab(out);
00211   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00212   VOTSME modelOutputTempState(model_,out,verbLevel);
00213   typedef Teuchos::VerboseObjectTempState<Thyra::DirectionalFiniteDiffCalculator<value_type> > VOTSDFDC;
00214   VOTSDFDC objDirecFiniteDiffCalculatorOutputTempState(objDirecFiniteDiffCalculator_,out,verbLevel);
00215   VOTSDFDC conDirecFiniteDiffCalculatorOutputTempState(conDirecFiniteDiffCalculator_,out,verbLevel);
00216   const bool trace = static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW);
00217   if(out.get() && trace)
00218     *out << "\nEntering MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00219   //
00220   // Validate input
00221   //
00222   TEST_FOR_EXCEPT(GcU!=NULL); // Can't handle these yet!
00223   TEST_FOR_EXCEPT(Uz!=NULL);
00224   TEST_FOR_EXCEPTION(
00225     objDirecFiniteDiffCalculator_.get()!=NULL && Gf!=NULL, std::logic_error
00226     ,"Error, can not compute full gradient vector Gf when using directional finite differences!"
00227     );
00228   //
00229   // Compute using derivative objects that come from the underlying model.
00230   //
00231   // Set the input and output arguments
00232   //
00233   // ToDo: Disallow computation Gf and instead just compute
00234   // the operators for this!
00235   //
00236   MEB::InArgs<value_type> model_inArgs = model_->createInArgs();
00237   MEB::OutArgs<value_type> model_outArgs = model_->createOutArgs();
00238   NLPObjGrad::ObjGradInfo obj_grad_info;
00239   obj_grad_info.Gf = Gf;
00240   obj_grad_info.f = f;
00241   if(recalc_c) obj_grad_info.c = c;
00242   preprocessBaseInOutArgs(
00243     x,true,NULL,&obj_grad_info,NULL
00244     ,&model_inArgs,&model_outArgs,NULL,NULL,NULL,NULL
00245     );
00246   if( py || rGf || D ) {
00247     if(thyra_C_.get()==NULL)
00248       thyra_C_ = model_->create_W();
00249     model_outArgs.set_W(thyra_C_.assert_not_null());
00250   }
00251   bool new_thyra_N = false;
00252   if( rGf || D ) {
00253     if(thyra_N_.get()==NULL) {
00254       thyra_N_ = Thyra::create_DfDp_mv(*model_,p_idx_,MEB::DERIV_MV_BY_COL).getMultiVector();
00255       new_thyra_N = true;
00256     }
00257     if(
00258       !conDirecFiniteDiffCalculator_.get()
00259       &&
00260       ( new_thyra_N
00261         || model_outArgs.get_DfDp_properties(p_idx_).linearity!=MEB::DERIV_LINEARITY_CONST )
00262       )
00263     {
00264       model_outArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00265     }
00266   }
00267   if ( !is_null(model_outArgs.get_W()) || !is_null(model_outArgs.get_W_op()) ) {
00268     if (model_inArgs.supports(MEB::IN_ARG_alpha))
00269       model_inArgs.set_alpha(0.0);
00270     if (model_inArgs.supports(MEB::IN_ARG_beta))
00271       model_inArgs.set_beta(1.0);
00272   }
00273   //
00274   // Evaluate the functions
00275   //
00276   model_->evalModel(model_inArgs,model_outArgs);
00277   //
00278   // Postprocess the evaluation
00279   //
00280   postprocessBaseOutArgs(&model_outArgs,Gf,f,recalc_c?c:NULL);
00281   // Setup solve components
00282   const VectorSpaceThyra *space_c;
00283   const VectorSpaceThyra *space_xD;
00284   Teuchos::RCP<const Thyra::VectorBase<value_type> > thyra_c;
00285   Teuchos::RCP<Thyra::VectorBase<value_type> > thyra_py;
00286   RCP<MatrixOp> D_used = rcp(D,false);
00287   RCP<Thyra::MultiVectorBase<value_type> > thyra_D;
00288 
00289   if( py || ( ( rGf || D ) && conDirecFiniteDiffCalculator_.get() ) ) {
00290     space_c = &dyn_cast<const VectorSpaceThyra>(c->space()),
00291       space_xD = &dyn_cast<const VectorSpaceThyra>(py->space());
00292     get_thyra_vector(*space_c,*c,&thyra_c);
00293     get_thyra_vector(*space_xD,py,&thyra_py);
00294   }
00295 
00296   if( D || rGf ) {
00297     if(!D) D_used = this->factory_D()->create();
00298     thyra_D =
00299       rcp_const_cast<Thyra::MultiVectorBase<value_type> >(
00300         rcp_dynamic_cast<const Thyra::MultiVectorBase<value_type> >(
00301           dyn_cast<MultiVectorMutableThyra>(*D_used).thyra_multi_vec()
00302           )
00303         );
00304   }
00305 
00306   if(
00307     ( rGf || D )
00308     &&
00309     !is_null(conDirecFiniteDiffCalculator_)
00310     &&
00311     ( new_thyra_N || !DfDp_is_const() )
00312     )
00313   {
00314     if(out.get() && trace)
00315       *out << "\nComputing thyra_N using directional finite differences ...\n";
00316     // !
00317     if(thyra_c.get())
00318       model_outArgs.set_f(rcp_const_cast<Thyra::VectorBase<value_type> >(thyra_c)); // Okay, will not be changed!
00319     typedef Thyra::DirectionalFiniteDiffCalculator<value_type>::SelectedDerivatives SelectedDerivatives; 
00320     MEB::OutArgs<value_type>
00321       model_fdOutArgs = conDirecFiniteDiffCalculator_->createOutArgs(
00322         *model_,SelectedDerivatives().supports(MEB::OUT_ARG_DfDp,0)
00323         );
00324     model_fdOutArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00325     conDirecFiniteDiffCalculator_->calcDerivatives(
00326       *model_
00327       ,model_inArgs // The base point to compute the derivatives at
00328       ,model_outArgs // The base function values
00329       ,model_fdOutArgs // The derivatives that we want to compute
00330       );
00331   }
00332   // Perform solve
00333   if( ( D || rGf ) && py ) {
00334     // Solve for py and D together
00335     if(out.get() && trace)
00336       *out << "\nSolving C*[py,D] = -[c,N] simultaneously ...\n";
00337     const int nind = thyra_N_->domain()->dim();
00338     RCP<Thyra::MultiVectorBase<value_type> > 
00339       thyra_cN = Thyra::createMembers(thyra_N_->range(),nind+1);
00340     Thyra::assign(thyra_cN->col(0).ptr(), *thyra_c);
00341     Thyra::assign(thyra_cN->subView(Teuchos::Range1D(1,nind)).ptr(), *thyra_N_);
00342     RCP<Thyra::MultiVectorBase<value_type> > 
00343       thyra_pyD = Thyra::createMembers(thyra_D->range(),nind+1);
00344     Thyra::assign(thyra_pyD.ptr(), 0.0);
00345     Thyra::SolveStatus<value_type>
00346       solveStatus = Thyra::solve(*thyra_C_, Thyra::NOTRANS, *thyra_cN, thyra_pyD.ptr());
00347     if(out.get() && trace) {
00348       *out
00349         << "\nsolve status:\n";
00350       OSTab(out).o() << solveStatus;
00351     }
00352     Thyra::scale(-1.0, thyra_pyD.ptr());
00353     Thyra::assign(thyra_py.ptr(), *thyra_pyD->col(0));
00354     Thyra::assign(thyra_D.ptr(), *thyra_pyD->subView(Teuchos::Range1D(1,nind)));
00355   }
00356   else {
00357     // Solve for py or D
00358     if( py ) {
00359       if(out.get() && trace)
00360         *out << "\nSolving C*py = -c ...\n";
00361       Thyra::assign(thyra_py.ptr(), 0.0);
00362       Thyra::SolveStatus<value_type> solveStatus =
00363         Thyra::solve<value_type>(*thyra_C_, Thyra::NOTRANS, *thyra_c, thyra_py.ptr());
00364       if (nonnull(out) && trace) {
00365         *out << "\nsolve status:\n";
00366         OSTab(out).o() << solveStatus;
00367       }
00368       Thyra::Vt_S(thyra_py.ptr(), -1.0);
00369     }
00370     if( D || rGf ) {
00371       if(out.get() && trace)
00372         *out << "\nSolving C*D = -N ...\n";
00373       Thyra::assign(thyra_D.ptr(), 0.0);
00374       Thyra::SolveStatus<value_type>
00375         solveStatus = Thyra::solve(*thyra_C_, Thyra::NOTRANS, *thyra_N_, thyra_D.ptr());
00376       if(out.get() && trace) {
00377         *out
00378           << "\nsolve status:\n";
00379         OSTab(out).o() << solveStatus;
00380       }
00381       Thyra::scale(-1.0, thyra_D.ptr());
00382     }
00383   }
00384   if(thyra_py.get()) {
00385     free_thyra_vector(*space_c,*c,&thyra_c);
00386     commit_thyra_vector(*space_xD,py,&thyra_py);
00387   }
00388   // Compute reduced gradient
00389   if(rGf) {
00390     // rGf = D' * Gf_xD + Gf_xI
00391     const Range1D
00392       var_dep = basis_sys_->var_dep(),
00393       var_indep = basis_sys_->var_indep();
00394     if(objDirecFiniteDiffCalculator_.get()) {
00395       if(out.get() && trace)
00396         *out << "\nComputing rGf using directional finite differences ...\n";
00397       //
00398       // Compute each component of the reduced gradient as
00399       //
00400       // rGf(i) = fd_product( model.g, [ D(:,i); e(i) ] )
00401       //
00402       AbstractLinAlgPack::VectorDenseMutableEncap d_rGf(*rGf);
00403       TVecPtr e_i = createMember(model_->get_p_space(p_idx_));
00404       Thyra::assign(e_i.ptr(), 0.0);
00405       MEB::InArgs<value_type> dir = model_->createInArgs();
00406       dir.set_p(p_idx_,e_i);
00407       MEB::OutArgs<value_type> bfunc = model_->createOutArgs();
00408       if(f) {
00409         bfunc.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00410         Thyra::set_ele(0, *f, bfunc.get_g(g_idx_).ptr());
00411       }
00412       MEB::OutArgs<value_type> var = model_->createOutArgs();
00413       var.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00414       const int np = var_indep.size();
00415       for( int i = 0; i < np; ++i ) {
00416         set_ele(i, 1.0, e_i.ptr());
00417         dir.set_x(thyra_D->col(i));
00418         objDirecFiniteDiffCalculator_->calcVariations(
00419           *model_,model_inArgs,dir,bfunc,var
00420           );
00421         dir.set_x(Teuchos::null);
00422         Thyra::set_ele(i, 0.0, e_i.ptr());
00423         d_rGf()[i] = Thyra::get_ele(*var.get_g(g_idx_),0);
00424       }
00425     }
00426     else {
00427       LinAlgOpPack::V_MtV( rGf, *D_used, BLAS_Cpp::trans, *Gf->sub_view(var_dep) );
00428       LinAlgOpPack::Vp_V( rGf, *Gf->sub_view(var_indep) );
00429       // ToDo: Just compute the operators associated with Gf and not Gf directly!
00430     }
00431   }
00432   // * ToDo: Add specialized algorithm for computing D using an inexact Jacobian
00433   // * ToDo: Add in logic for inexact solves
00434   if(out.get() && trace)
00435     *out << "\nLeaving MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00436 }
00437 
00438 void NLPDirectThyraModelEvaluator::calc_semi_newton_step(
00439   const Vector &x
00440   ,VectorMutable *c
00441   ,bool recalc_c
00442   ,VectorMutable *py
00443   ) const
00444 {
00445   if(recalc_c) {
00446     //
00447     // Recompute c
00448     //
00449     TEST_FOR_EXCEPT(true);
00450   }
00451   // Compute py = - inv(C)*c
00452   TEST_FOR_EXCEPT(true);
00453 }
00454 
00455 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends