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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "NLPInterfacePack_NLPDirectThyraModelEvaluator.hpp"
00043 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00044 #include "AbstractLinAlgPack_VectorOut.hpp"
00045 #include "AbstractLinAlgPack_ThyraAccessors.hpp"
00046 #include "AbstractLinAlgPack_VectorSpaceThyra.hpp"
00047 #include "AbstractLinAlgPack_VectorMutableThyra.hpp"
00048 #include "AbstractLinAlgPack_MultiVectorMutableThyra.hpp"
00049 #include "AbstractLinAlgPack_MatrixOpNonsingThyra.hpp"
00050 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00051 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00052 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00053 #include "AbstractLinAlgPack_BasisSystem.hpp"
00054 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00055 #include "Thyra_ModelEvaluatorHelpers.hpp"
00056 #include "Thyra_DetachedVectorView.hpp"
00057 #include "Thyra_VectorStdOps.hpp"
00058 #include "Thyra_MultiVectorStdOps.hpp"
00059 #include "Teuchos_AbstractFactoryStd.hpp"
00060 #include "Teuchos_Assert.hpp"
00061 #include "Teuchos_dyn_cast.hpp"
00062 
00063 namespace NLPInterfacePack {
00064 
00065 NLPDirectThyraModelEvaluator::NLPDirectThyraModelEvaluator()
00066   :DfDp_is_const_(false)
00067 {}
00068 
00069 NLPDirectThyraModelEvaluator::NLPDirectThyraModelEvaluator(
00070   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> > &model 
00071   ,const int p_idx 
00072   ,const int g_idx 
00073   ,const objDirecFiniteDiffCalculator_ptr_t objDirecFiniteDiffCalculator
00074   ,const conDirecFiniteDiffCalculator_ptr_t conDirecFiniteDiffCalculator
00075   )
00076   :DfDp_is_const_(false)
00077 {
00078   initialize(
00079     model,p_idx,g_idx
00080     ,objDirecFiniteDiffCalculator,conDirecFiniteDiffCalculator
00081     );
00082 }
00083 
00084 void NLPDirectThyraModelEvaluator::initialize(
00085   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> > &model
00086   ,const int p_idx
00087   ,const int g_idx
00088   ,const objDirecFiniteDiffCalculator_ptr_t objDirecFiniteDiffCalculator
00089   ,const conDirecFiniteDiffCalculator_ptr_t conDirecFiniteDiffCalculator
00090   )
00091 {
00092   typedef Thyra::ModelEvaluatorBase MEB;
00093   if(objDirecFiniteDiffCalculator.get())
00094     this->set_objDirecFiniteDiffCalculator(objDirecFiniteDiffCalculator);
00095   if(conDirecFiniteDiffCalculator.get())
00096     this->set_conDirecFiniteDiffCalculator(conDirecFiniteDiffCalculator);
00097   initializeBase(model,p_idx,g_idx);
00098   Thyra::ModelEvaluatorBase::OutArgs<double> model_outArgs = model->createOutArgs();
00099   MEB::DerivativeProperties model_W_properties = model_outArgs.get_W_properties();
00100   if( p_idx >= 0 ) {
00101     TEUCHOS_TEST_FOR_EXCEPTION(
00102       (conDirecFiniteDiffCalculator_.get()==0 && !model_outArgs.supports(MEB::OUT_ARG_DfDp,p_idx).supports(MEB::DERIV_MV_BY_COL))
00103       ,std::invalid_argument
00104       ,"Error, model must support computing DfDp("<<p_idx<<") as a"
00105       " column-oriented multi-vector if not using finite differences!"
00106       );
00107   }
00108   DfDp_is_const_  = false;
00109 }
00110 
00111 // Overridden public members from NLP
00112 
00113 void NLPDirectThyraModelEvaluator::initialize(bool test_setup)
00114 {
00115   if(initialized_) {
00116     NLPDirect::initialize(test_setup);
00117     return;
00118   }
00119   NLPThyraModelEvaluatorBase::initialize(test_setup);
00120   NLPDirect::initialize(test_setup);
00121 }
00122 
00123 void NLPDirectThyraModelEvaluator::unset_quantities()
00124 {
00125   NLPDirect::unset_quantities();
00126 }
00127 
00128 // Overridden public members from NLPObjGrad
00129 
00130 bool NLPDirectThyraModelEvaluator::supports_Gf() const
00131 {
00132   if(objDirecFiniteDiffCalculator_.get())
00133     return false;
00134   return true;
00135   // ToDo: Change this to false when only the operator for model_DgDx is
00136   // supported!
00137 }
00138 
00139 bool NLPDirectThyraModelEvaluator::supports_Gf_prod() const
00140 {
00141   return ( objDirecFiniteDiffCalculator_.get() != NULL );
00142 }
00143 
00144 value_type NLPDirectThyraModelEvaluator::calc_Gf_prod(
00145   const Vector& x, const Vector& d, bool newx
00146   ) const
00147 {
00148   TEUCHOS_TEST_FOR_EXCEPT(objDirecFiniteDiffCalculator_.get()==0);
00149   typedef Thyra::ModelEvaluatorBase MEB;
00150   MEB::InArgs<value_type> basePoint = model_->createInArgs();
00151   MEB::InArgs<value_type> directions = model_->createInArgs();
00152   MEB::OutArgs<value_type> baseFunc = model_->createOutArgs();
00153   MEB::OutArgs<value_type> variations = model_->createOutArgs();
00154   set_x(x,&basePoint);
00155   set_x(d,&directions);
00156   variations.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00157   objDirecFiniteDiffCalculator_->calcVariations(
00158     *model_,basePoint,directions,baseFunc,variations
00159     );
00160   return Thyra::get_ele(*variations.get_g(g_idx_),0);
00161 }
00162 
00163 // Overridden public members from NLPDirect
00164 
00165 Range1D NLPDirectThyraModelEvaluator::var_dep() const
00166 {
00167   return basis_sys_->var_dep();
00168 }
00169 
00170 Range1D NLPDirectThyraModelEvaluator::var_indep() const
00171 {
00172   return basis_sys_->var_indep();
00173 }
00174 
00175 const NLPDirect::mat_fcty_ptr_t
00176 NLPDirectThyraModelEvaluator::factory_D() const
00177 {
00178   return basis_sys_->factory_D();
00179 }
00180 
00181 const NLPDirect::mat_sym_nonsing_fcty_ptr_t
00182 NLPDirectThyraModelEvaluator::factory_S() const
00183 {
00184   return basis_sys_->factory_S();
00185 }
00186 
00187 void NLPDirectThyraModelEvaluator::calc_point(
00188   const Vector &x
00189   ,value_type *f
00190   ,VectorMutable *c
00191   ,bool recalc_c
00192   ,VectorMutable *Gf
00193   ,VectorMutable *py
00194   ,VectorMutable *rGf
00195   ,MatrixOp *GcU
00196   ,MatrixOp *D
00197   ,MatrixOp *Uz
00198   ) const
00199 {
00200   using Teuchos::FancyOStream;
00201   using Teuchos::OSTab;
00202   using Teuchos::dyn_cast;
00203   using Teuchos::RCP;
00204   using Teuchos::rcp;
00205   using Teuchos::rcp_const_cast;
00206   using Teuchos::rcp_dynamic_cast;
00207   using AbstractLinAlgPack::VectorSpaceThyra;
00208   using AbstractLinAlgPack::VectorMutableThyra;
00209   using AbstractLinAlgPack::MultiVectorMutableThyra;
00210   using AbstractLinAlgPack::MatrixOpThyra;
00211   using AbstractLinAlgPack::MatrixOpNonsingThyra;
00212   typedef Thyra::ModelEvaluatorBase MEB;
00213   typedef Teuchos::RCP<Thyra::VectorBase<value_type> > TVecPtr;
00214   typedef MEB::DerivativeMultiVector<value_type> DerivMV;
00215   typedef MEB::Derivative<value_type> Deriv;
00216   //
00217   // Get output and verbosity
00218   //
00219   const Teuchos::RCP<Teuchos::FancyOStream>
00220     out = this->getOStream();
00221   const Teuchos::EVerbosityLevel
00222     verbLevel = ( showModelEvaluatorTrace() ? this->getVerbLevel() : Teuchos::VERB_NONE );
00223   Teuchos::OSTab tab(out);
00224   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00225   VOTSME modelOutputTempState(model_,out,verbLevel);
00226   typedef Teuchos::VerboseObjectTempState<Thyra::DirectionalFiniteDiffCalculator<value_type> > VOTSDFDC;
00227   VOTSDFDC objDirecFiniteDiffCalculatorOutputTempState(objDirecFiniteDiffCalculator_,out,verbLevel);
00228   VOTSDFDC conDirecFiniteDiffCalculatorOutputTempState(conDirecFiniteDiffCalculator_,out,verbLevel);
00229   const bool trace = static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW);
00230   if(out.get() && trace)
00231     *out << "\nEntering MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00232   //
00233   // Validate input
00234   //
00235   TEUCHOS_TEST_FOR_EXCEPT(GcU!=NULL); // Can't handle these yet!
00236   TEUCHOS_TEST_FOR_EXCEPT(Uz!=NULL);
00237   TEUCHOS_TEST_FOR_EXCEPTION(
00238     objDirecFiniteDiffCalculator_.get()!=NULL && Gf!=NULL, std::logic_error
00239     ,"Error, can not compute full gradient vector Gf when using directional finite differences!"
00240     );
00241   //
00242   // Compute using derivative objects that come from the underlying model.
00243   //
00244   // Set the input and output arguments
00245   //
00246   // ToDo: Disallow computation Gf and instead just compute
00247   // the operators for this!
00248   //
00249   MEB::InArgs<value_type> model_inArgs = model_->createInArgs();
00250   MEB::OutArgs<value_type> model_outArgs = model_->createOutArgs();
00251   NLPObjGrad::ObjGradInfo obj_grad_info;
00252   obj_grad_info.Gf = Gf;
00253   obj_grad_info.f = f;
00254   if(recalc_c) obj_grad_info.c = c;
00255   preprocessBaseInOutArgs(
00256     x,true,NULL,&obj_grad_info,NULL
00257     ,&model_inArgs,&model_outArgs,NULL,NULL,NULL,NULL
00258     );
00259   if( py || rGf || D ) {
00260     if(thyra_C_.get()==NULL)
00261       thyra_C_ = model_->create_W();
00262     model_outArgs.set_W(thyra_C_.assert_not_null());
00263   }
00264   bool new_thyra_N = false;
00265   if( rGf || D ) {
00266     if(thyra_N_.get()==NULL) {
00267       thyra_N_ = Thyra::create_DfDp_mv(*model_,p_idx_,MEB::DERIV_MV_BY_COL).getMultiVector();
00268       new_thyra_N = true;
00269     }
00270     if(
00271       !conDirecFiniteDiffCalculator_.get()
00272       &&
00273       ( new_thyra_N
00274         || model_outArgs.get_DfDp_properties(p_idx_).linearity!=MEB::DERIV_LINEARITY_CONST )
00275       )
00276     {
00277       model_outArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00278     }
00279   }
00280   if ( !is_null(model_outArgs.get_W()) || !is_null(model_outArgs.get_W_op()) ) {
00281     if (model_inArgs.supports(MEB::IN_ARG_alpha))
00282       model_inArgs.set_alpha(0.0);
00283     if (model_inArgs.supports(MEB::IN_ARG_beta))
00284       model_inArgs.set_beta(1.0);
00285   }
00286   //
00287   // Evaluate the functions
00288   //
00289   model_->evalModel(model_inArgs,model_outArgs);
00290   //
00291   // Postprocess the evaluation
00292   //
00293   postprocessBaseOutArgs(&model_outArgs,Gf,f,recalc_c?c:NULL);
00294   // Setup solve components
00295   const VectorSpaceThyra *space_c;
00296   const VectorSpaceThyra *space_xD;
00297   Teuchos::RCP<const Thyra::VectorBase<value_type> > thyra_c;
00298   Teuchos::RCP<Thyra::VectorBase<value_type> > thyra_py;
00299   RCP<MatrixOp> D_used = rcp(D,false);
00300   RCP<Thyra::MultiVectorBase<value_type> > thyra_D;
00301 
00302   if( py || ( ( rGf || D ) && conDirecFiniteDiffCalculator_.get() ) ) {
00303     space_c = &dyn_cast<const VectorSpaceThyra>(c->space()),
00304       space_xD = &dyn_cast<const VectorSpaceThyra>(py->space());
00305     get_thyra_vector(*space_c,*c,&thyra_c);
00306     get_thyra_vector(*space_xD,py,&thyra_py);
00307   }
00308 
00309   if( D || rGf ) {
00310     if(!D) D_used = this->factory_D()->create();
00311     thyra_D =
00312       rcp_const_cast<Thyra::MultiVectorBase<value_type> >(
00313         rcp_dynamic_cast<const Thyra::MultiVectorBase<value_type> >(
00314           dyn_cast<MultiVectorMutableThyra>(*D_used).thyra_multi_vec()
00315           )
00316         );
00317   }
00318 
00319   if(
00320     ( rGf || D )
00321     &&
00322     !is_null(conDirecFiniteDiffCalculator_)
00323     &&
00324     ( new_thyra_N || !DfDp_is_const() )
00325     )
00326   {
00327     if(out.get() && trace)
00328       *out << "\nComputing thyra_N using directional finite differences ...\n";
00329     // !
00330     if(thyra_c.get())
00331       model_outArgs.set_f(rcp_const_cast<Thyra::VectorBase<value_type> >(thyra_c)); // Okay, will not be changed!
00332     typedef Thyra::DirectionalFiniteDiffCalculator<value_type>::SelectedDerivatives SelectedDerivatives; 
00333     MEB::OutArgs<value_type>
00334       model_fdOutArgs = conDirecFiniteDiffCalculator_->createOutArgs(
00335         *model_,SelectedDerivatives().supports(MEB::OUT_ARG_DfDp,0)
00336         );
00337     model_fdOutArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00338     conDirecFiniteDiffCalculator_->calcDerivatives(
00339       *model_
00340       ,model_inArgs // The base point to compute the derivatives at
00341       ,model_outArgs // The base function values
00342       ,model_fdOutArgs // The derivatives that we want to compute
00343       );
00344   }
00345   // Perform solve
00346   if( ( D || rGf ) && py ) {
00347     // Solve for py and D together
00348     if(out.get() && trace)
00349       *out << "\nSolving C*[py,D] = -[c,N] simultaneously ...\n";
00350     const int nind = thyra_N_->domain()->dim();
00351     RCP<Thyra::MultiVectorBase<value_type> > 
00352       thyra_cN = Thyra::createMembers(thyra_N_->range(),nind+1);
00353     Thyra::assign(thyra_cN->col(0).ptr(), *thyra_c);
00354     Thyra::assign(thyra_cN->subView(Teuchos::Range1D(1,nind)).ptr(), *thyra_N_);
00355     RCP<Thyra::MultiVectorBase<value_type> > 
00356       thyra_pyD = Thyra::createMembers(thyra_D->range(),nind+1);
00357     Thyra::assign(thyra_pyD.ptr(), 0.0);
00358     Thyra::SolveStatus<value_type>
00359       solveStatus = Thyra::solve(*thyra_C_, Thyra::NOTRANS, *thyra_cN, thyra_pyD.ptr());
00360     if(out.get() && trace) {
00361       *out
00362         << "\nsolve status:\n";
00363       OSTab(out).o() << solveStatus;
00364     }
00365     Thyra::scale(-1.0, thyra_pyD.ptr());
00366     Thyra::assign(thyra_py.ptr(), *thyra_pyD->col(0));
00367     Thyra::assign(thyra_D.ptr(), *thyra_pyD->subView(Teuchos::Range1D(1,nind)));
00368   }
00369   else {
00370     // Solve for py or D
00371     if( py ) {
00372       if(out.get() && trace)
00373         *out << "\nSolving C*py = -c ...\n";
00374       Thyra::assign(thyra_py.ptr(), 0.0);
00375       Thyra::SolveStatus<value_type> solveStatus =
00376         Thyra::solve<value_type>(*thyra_C_, Thyra::NOTRANS, *thyra_c, thyra_py.ptr());
00377       if (nonnull(out) && trace) {
00378         *out << "\nsolve status:\n";
00379         OSTab(out).o() << solveStatus;
00380       }
00381       Thyra::Vt_S(thyra_py.ptr(), -1.0);
00382     }
00383     if( D || rGf ) {
00384       if(out.get() && trace)
00385         *out << "\nSolving C*D = -N ...\n";
00386       Thyra::assign(thyra_D.ptr(), 0.0);
00387       Thyra::SolveStatus<value_type>
00388         solveStatus = Thyra::solve(*thyra_C_, Thyra::NOTRANS, *thyra_N_, thyra_D.ptr());
00389       if(out.get() && trace) {
00390         *out
00391           << "\nsolve status:\n";
00392         OSTab(out).o() << solveStatus;
00393       }
00394       Thyra::scale(-1.0, thyra_D.ptr());
00395     }
00396   }
00397   if(thyra_py.get()) {
00398     free_thyra_vector(*space_c,*c,&thyra_c);
00399     commit_thyra_vector(*space_xD,py,&thyra_py);
00400   }
00401   // Compute reduced gradient
00402   if(rGf) {
00403     // rGf = D' * Gf_xD + Gf_xI
00404     const Range1D
00405       var_dep = basis_sys_->var_dep(),
00406       var_indep = basis_sys_->var_indep();
00407     if(objDirecFiniteDiffCalculator_.get()) {
00408       if(out.get() && trace)
00409         *out << "\nComputing rGf using directional finite differences ...\n";
00410       //
00411       // Compute each component of the reduced gradient as
00412       //
00413       // rGf(i) = fd_product( model.g, [ D(:,i); e(i) ] )
00414       //
00415       AbstractLinAlgPack::VectorDenseMutableEncap d_rGf(*rGf);
00416       TVecPtr e_i = createMember(model_->get_p_space(p_idx_));
00417       Thyra::assign(e_i.ptr(), 0.0);
00418       MEB::InArgs<value_type> dir = model_->createInArgs();
00419       dir.set_p(p_idx_,e_i);
00420       MEB::OutArgs<value_type> bfunc = model_->createOutArgs();
00421       if(f) {
00422         bfunc.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00423         Thyra::set_ele(0, *f, bfunc.get_g(g_idx_).ptr());
00424       }
00425       MEB::OutArgs<value_type> var = model_->createOutArgs();
00426       var.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00427       const int np = var_indep.size();
00428       for( int i = 0; i < np; ++i ) {
00429         set_ele(i, 1.0, e_i.ptr());
00430         dir.set_x(thyra_D->col(i));
00431         objDirecFiniteDiffCalculator_->calcVariations(
00432           *model_,model_inArgs,dir,bfunc,var
00433           );
00434         dir.set_x(Teuchos::null);
00435         Thyra::set_ele(i, 0.0, e_i.ptr());
00436         d_rGf()[i] = Thyra::get_ele(*var.get_g(g_idx_),0);
00437       }
00438     }
00439     else {
00440       LinAlgOpPack::V_MtV( rGf, *D_used, BLAS_Cpp::trans, *Gf->sub_view(var_dep) );
00441       LinAlgOpPack::Vp_V( rGf, *Gf->sub_view(var_indep) );
00442       // ToDo: Just compute the operators associated with Gf and not Gf directly!
00443     }
00444   }
00445   // * ToDo: Add specialized algorithm for computing D using an inexact Jacobian
00446   // * ToDo: Add in logic for inexact solves
00447   if(out.get() && trace)
00448     *out << "\nLeaving MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00449 }
00450 
00451 void NLPDirectThyraModelEvaluator::calc_semi_newton_step(
00452   const Vector &x
00453   ,VectorMutable *c
00454   ,bool recalc_c
00455   ,VectorMutable *py
00456   ) const
00457 {
00458   if(recalc_c) {
00459     //
00460     // Recompute c
00461     //
00462     TEUCHOS_TEST_FOR_EXCEPT(true);
00463   }
00464   // Compute py = - inv(C)*c
00465   TEUCHOS_TEST_FOR_EXCEPT(true);
00466 }
00467 
00468 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends