Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultInverseModelEvaluator.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP
00030 #define THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP
00031 
00032 
00033 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00034 #include "Thyra_ModelEvaluatorHelpers.hpp"
00035 #include "Thyra_DetachedVectorView.hpp"
00036 #include "Thyra_ParameterDrivenMultiVectorInput.hpp"
00037 #include "Thyra_VectorSpaceFactoryBase.hpp"
00038 #include "Thyra_MultiVectorStdOps.hpp"
00039 #include "Thyra_AssertOp.hpp"
00040 #include "Teuchos_StandardMemberCompositionMacros.hpp"
00041 #include "Teuchos_StandardCompositionMacros.hpp"
00042 #include "Teuchos_ParameterListAcceptor.hpp"
00043 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00044 #include "Teuchos_StandardParameterEntryValidators.hpp"
00045 #include "Teuchos_Time.hpp"
00046 
00047 
00048 namespace Thyra {
00049 
00050 
00249 template<class Scalar>
00250 class DefaultInverseModelEvaluator
00251   : virtual public ModelEvaluatorDelegatorBase<Scalar>
00252   , virtual public Teuchos::ParameterListAcceptor
00253 {
00254 public:
00255 
00257   STANDARD_CONST_COMPOSITION_MEMBERS( VectorBase<Scalar>, observationTarget );
00258 
00260   STANDARD_CONST_COMPOSITION_MEMBERS( VectorBase<Scalar>, parameterBase );
00261 
00263   STANDARD_CONST_COMPOSITION_MEMBERS( LinearOpBase<Scalar>, observationMatchWeightingOp );
00264 
00266   STANDARD_CONST_COMPOSITION_MEMBERS( LinearOpBase<Scalar>, parameterRegularizationWeightingOp );
00267 
00270   STANDARD_NONCONST_COMPOSITION_MEMBERS( MultiVectorFileIOBase<Scalar>, observationTargetIO );
00271 
00274   STANDARD_NONCONST_COMPOSITION_MEMBERS( MultiVectorFileIOBase<Scalar>, parameterBaseIO );
00275 
00278 
00280   DefaultInverseModelEvaluator();
00281 
00283   void initialize(
00284     const RCP<ModelEvaluator<Scalar> > &thyraModel
00285     );
00286 
00288   void uninitialize(
00289     RCP<ModelEvaluator<Scalar> > *thyraModel
00290     );
00291 
00293 
00296 
00304   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00306   RCP<Teuchos::ParameterList> getNonconstParameterList();
00308   RCP<Teuchos::ParameterList> unsetParameterList();
00310   RCP<const Teuchos::ParameterList> getParameterList() const;
00319   RCP<const Teuchos::ParameterList> getValidParameters() const;
00320 
00322 
00325 
00327   std::string description() const;
00328 
00330 
00333 
00335   RCP<const VectorSpaceBase<Scalar> > get_p_space(int l) const;
00337   RCP<const VectorSpaceBase<Scalar> > get_g_space(int j) const;
00339   ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00340 
00342 
00343 private:
00344 
00347 
00349   ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00351   void evalModelImpl(
00352     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00353     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00354     ) const;
00355 
00357 
00358 private:
00359 
00360   // ////////////////////////////////
00361   // Private data members
00362 
00363   mutable RCP<const Teuchos::ParameterList> validParamList_;
00364   RCP<Teuchos::ParameterList>  paramList_;
00365 
00366   RCP<const VectorSpaceBase<Scalar> > inv_g_space_;
00367 
00368   mutable ModelEvaluatorBase::InArgs<Scalar> prototypeInArgs_;
00369   mutable ModelEvaluatorBase::OutArgs<Scalar> prototypeOutArgs_;
00370   mutable bool usingObservationTargetAsParameter_;
00371 
00372   int obs_idx_;
00373   int p_idx_;
00374 
00375 
00376   double observationMultiplier_;
00377   double parameterMultiplier_; 
00378 
00379   bool observationTargetAsParameter_;
00380   
00381   bool observationPassThrough_;
00382 
00383   Teuchos::EVerbosityLevel localVerbLevel_;
00384 
00385   mutable ParameterDrivenMultiVectorInput<Scalar> observationTargetReader_;
00386   mutable ParameterDrivenMultiVectorInput<Scalar> parameterBaseReader_;
00387 
00388   static const std::string ObservationIndex_name_;
00389   static const int ObservationIndex_default_;
00390 
00391   static const std::string ParameterSubvectorIndex_name_;
00392   static const int ParameterSubvectorIndex_default_;
00393 
00394   static const std::string ObservationMultiplier_name_;
00395   static const double ObservationMultiplier_default_;
00396 
00397   static const std::string ObservationTargetVector_name_;
00398 
00399   static const std::string ObservationTargetAsParameter_name_;
00400   static const bool ObservationTargetAsParameter_default_;
00401 
00402   static const std::string ObservationPassThrough_name_;
00403   static const bool ObservationPassThrough_default_;
00404 
00405   static const std::string ParameterMultiplier_name_;
00406   static const double ParameterMultiplier_default_;
00407 
00408   static const std::string ParameterBaseVector_name_;
00409 
00410   // ////////////////////////////////
00411   // Private member functions
00412 
00413   void initializeDefaults();
00414 
00415   void initializeInArgsOutArgs() const;
00416 
00417   RCP<const VectorSpaceBase<Scalar> > get_obs_space() const;
00418 
00419 };
00420 
00421 
00426 template<class Scalar>
00427 RCP<DefaultInverseModelEvaluator<Scalar> >
00428 defaultInverseModelEvaluator(
00429   const RCP<ModelEvaluator<Scalar> > &thyraModel
00430   )
00431 {
00432   RCP<DefaultInverseModelEvaluator<Scalar> >
00433     inverseModel = Teuchos::rcp(new DefaultInverseModelEvaluator<Scalar>);
00434   inverseModel->initialize(thyraModel);
00435   return inverseModel;
00436 }
00437 
00438 
00439 // /////////////////////////////////
00440 // Implementations
00441 
00442 
00443 // Static data members
00444 
00445 
00446 template<class Scalar>
00447 const std::string
00448 DefaultInverseModelEvaluator<Scalar>::ObservationIndex_name_
00449 = "Observation Index";
00450 
00451 template<class Scalar>
00452 const int
00453 DefaultInverseModelEvaluator<Scalar>::ObservationIndex_default_
00454 = -1;
00455 
00456 
00457 template<class Scalar>
00458 const std::string
00459 DefaultInverseModelEvaluator<Scalar>::ParameterSubvectorIndex_name_
00460 = "Parameter Subvector Ordinal";
00461 
00462 template<class Scalar>
00463 const int
00464 DefaultInverseModelEvaluator<Scalar>::ParameterSubvectorIndex_default_
00465 = 0;
00466 
00467 
00468 template<class Scalar>
00469 const std::string
00470 DefaultInverseModelEvaluator<Scalar>::ObservationMultiplier_name_
00471 = "Observation Multiplier";
00472 
00473 template<class Scalar>
00474 const double
00475 DefaultInverseModelEvaluator<Scalar>::ObservationMultiplier_default_
00476 = 1.0;
00477 
00478 
00479 template<class Scalar>
00480 const std::string
00481 DefaultInverseModelEvaluator<Scalar>::ObservationTargetVector_name_
00482 = "Observation Target Vector";
00483 
00484 
00485 template<class Scalar>
00486 const std::string
00487 DefaultInverseModelEvaluator<Scalar>::ObservationTargetAsParameter_name_
00488 = "Observation Target as Parameter";
00489 
00490 template<class Scalar>
00491 const bool
00492 DefaultInverseModelEvaluator<Scalar>::ObservationTargetAsParameter_default_
00493 = false;
00494 
00495 
00496 template<class Scalar>
00497 const std::string
00498 DefaultInverseModelEvaluator<Scalar>::ObservationPassThrough_name_
00499 = "Observation Pass Through";
00500 
00501 template<class Scalar>
00502 const bool
00503 DefaultInverseModelEvaluator<Scalar>::ObservationPassThrough_default_
00504 = false;
00505 
00506 
00507 template<class Scalar>
00508 const std::string
00509 DefaultInverseModelEvaluator<Scalar>::ParameterMultiplier_name_
00510 = "Parameter Multiplier";
00511 
00512 template<class Scalar>
00513 const double
00514 DefaultInverseModelEvaluator<Scalar>::ParameterMultiplier_default_
00515 = 1e-6;
00516 
00517 
00518 template<class Scalar>
00519 const std::string
00520 DefaultInverseModelEvaluator<Scalar>::ParameterBaseVector_name_
00521 = "Parameter Base Vector";
00522 
00523 
00524 // Constructors/initializers/accessors/utilities
00525 
00526 
00527 template<class Scalar>
00528 DefaultInverseModelEvaluator<Scalar>::DefaultInverseModelEvaluator()
00529   :usingObservationTargetAsParameter_(false), obs_idx_(-1),p_idx_(0),
00530    observationTargetAsParameter_(false),
00531    observationPassThrough_(ObservationPassThrough_default_),
00532    localVerbLevel_(Teuchos::VERB_DEFAULT)
00533 {}
00534 
00535 
00536 template<class Scalar>
00537 void DefaultInverseModelEvaluator<Scalar>::initialize(
00538   const RCP<ModelEvaluator<Scalar> > &thyraModel
00539   )
00540 {
00541   this->ModelEvaluatorDelegatorBase<Scalar>::initialize(thyraModel);
00542   inv_g_space_= thyraModel->get_x_space()->smallVecSpcFcty()->createVecSpc(1);
00543   // Get ready for reinitalization
00544   prototypeInArgs_ = ModelEvaluatorBase::InArgs<Scalar>();
00545   prototypeOutArgs_ = ModelEvaluatorBase::OutArgs<Scalar>();
00546 }
00547 
00548 
00549 template<class Scalar>
00550 void DefaultInverseModelEvaluator<Scalar>::uninitialize(
00551   RCP<ModelEvaluator<Scalar> >  *thyraModel
00552   )
00553 {
00554   if(thyraModel) *thyraModel = this->getUnderlyingModel();
00555   this->ModelEvaluatorDelegatorBase<Scalar>::uninitialize();
00556 }
00557 
00558 
00559 // Overridden from Teuchos::ParameterListAcceptor
00560 
00561 
00562 template<class Scalar>
00563 void DefaultInverseModelEvaluator<Scalar>::setParameterList(
00564   RCP<Teuchos::ParameterList> const& paramList
00565   )
00566 {
00567 
00568   using Teuchos::Array;
00569   using Teuchos::getParameterPtr;
00570   using Teuchos::rcp;
00571   using Teuchos::sublist;
00572 
00573   // Validate and set the parameter list
00574   TEST_FOR_EXCEPT(0==paramList.get());
00575   paramList->validateParameters(*getValidParameters(),0);
00576   paramList_ = paramList;
00577 
00578   // Parameters for observation matching term
00579   obs_idx_ = paramList_->get(
00580     ObservationIndex_name_,ObservationIndex_default_);
00581   observationPassThrough_ = paramList_->get(
00582     ObservationPassThrough_name_, ObservationPassThrough_default_ );
00583 #ifdef TEUCHOS_DEBUG
00584     TEST_FOR_EXCEPTION(
00585       ( obs_idx_ < 0 &&  observationPassThrough_ ), std::logic_error,
00586       "Error, the observation function index obs_idx = " << obs_idx_ << " is not\n"
00587       "allowed when the observation is simply passed through!"
00588       );
00589 #endif
00590   observationMultiplier_ = paramList_->get(
00591     ObservationMultiplier_name_,ObservationMultiplier_default_);
00592   if (!ObservationPassThrough_default_) {
00593     observationTargetAsParameter_ = paramList_->get(
00594       ObservationTargetAsParameter_name_, ObservationTargetAsParameter_default_ );
00595     if(get_observationTargetIO().get()) {
00596       observationTargetReader_.set_vecSpc(get_obs_space());
00597       Teuchos::VerboseObjectTempState<ParameterDrivenMultiVectorInput<Scalar> >
00598         vots_observationTargetReader(
00599           rcp(&observationTargetReader_,false)
00600           ,this->getOStream(),this->getVerbLevel()
00601           );
00602       observationTargetReader_.setParameterList(
00603         sublist(paramList_,ObservationTargetVector_name_)
00604         );
00605       RCP<VectorBase<Scalar> >
00606         observationTarget;
00607       observationTargetReader_.readVector(
00608         "observation target vector",&observationTarget);
00609       observationTarget_ = observationTarget;
00610     }
00611   }
00612   else {
00613     observationTargetAsParameter_ = false;
00614     observationTarget_ = Teuchos::null;
00615   }
00616   
00617   // Parameters for parameter matching term
00618   p_idx_ = paramList_->get(
00619     ParameterSubvectorIndex_name_,ParameterSubvectorIndex_default_);
00620   parameterMultiplier_ = paramList_->get(
00621     ParameterMultiplier_name_,ParameterMultiplier_default_);
00622   if(get_parameterBaseIO().get()) {
00623     parameterBaseReader_.set_vecSpc(this->get_p_space(p_idx_));
00624     Teuchos::VerboseObjectTempState<ParameterDrivenMultiVectorInput<Scalar> >
00625       vots_parameterBaseReader(
00626         rcp(&parameterBaseReader_,false)
00627         ,this->getOStream(),this->getVerbLevel()
00628         );
00629     parameterBaseReader_.setParameterList(
00630       sublist(paramList_,ParameterBaseVector_name_)
00631       );
00632     RCP<VectorBase<Scalar> >
00633       parameterBase;
00634     parameterBaseReader_.readVector(
00635       "parameter base vector",&parameterBase);
00636     parameterBase_ = parameterBase;
00637   }
00638 
00639   // Verbosity settings
00640   localVerbLevel_ = this->readLocalVerbosityLevelValidatedParameter(*paramList_);
00641   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00642 
00643 #ifdef TEUCHOS_DEBUG
00644   paramList_->validateParameters(*getValidParameters(),0);
00645 #endif // TEUCHOS_DEBUG
00646 
00647   // Get ready for reinitalization
00648   prototypeInArgs_ = ModelEvaluatorBase::InArgs<Scalar>();
00649   prototypeOutArgs_ = ModelEvaluatorBase::OutArgs<Scalar>();
00650 
00651 }
00652 
00653 
00654 template<class Scalar>
00655 RCP<Teuchos::ParameterList>
00656 DefaultInverseModelEvaluator<Scalar>::getNonconstParameterList()
00657 {
00658   return paramList_;
00659 }
00660 
00661 
00662 template<class Scalar>
00663 RCP<Teuchos::ParameterList>
00664 DefaultInverseModelEvaluator<Scalar>::unsetParameterList()
00665 {
00666   RCP<Teuchos::ParameterList> _paramList = paramList_;
00667   paramList_ = Teuchos::null;
00668   return _paramList;
00669 }
00670 
00671 
00672 template<class Scalar>
00673 RCP<const Teuchos::ParameterList>
00674 DefaultInverseModelEvaluator<Scalar>::getParameterList() const
00675 {
00676   return paramList_;
00677 }
00678 
00679 
00680 template<class Scalar>
00681 RCP<const Teuchos::ParameterList>
00682 DefaultInverseModelEvaluator<Scalar>::getValidParameters() const
00683 {
00684   if(validParamList_.get()==NULL) {
00685     RCP<Teuchos::ParameterList>
00686       pl = Teuchos::rcp(new Teuchos::ParameterList());
00687     pl->set( ObservationIndex_name_,ObservationIndex_default_,
00688       "The index of the observation function, obs_idx.\n"
00689       "If obs_idx < 0, then the observation will be the state vector x.\n"
00690       "If obs_idx >= 0, then the observation will be the response function g(obs_idx)."
00691       );
00692     pl->set( ParameterSubvectorIndex_name_,ParameterSubvectorIndex_default_,
00693       "The index of the parameter subvector that will be used in the\n"
00694       "regularization term."
00695       );
00696     pl->set( ObservationMultiplier_name_,ObservationMultiplier_default_,
00697       "observationMultiplier"
00698       );
00699     if(this->get_observationTargetIO().get())
00700       observationTargetReader_.set_fileIO(this->get_observationTargetIO());
00701     pl->sublist(ObservationTargetVector_name_).setParameters(
00702       *observationTargetReader_.getValidParameters()
00703       );
00704     pl->set( ObservationPassThrough_name_, ObservationPassThrough_default_,
00705       "If true, then the observation will just be used instead of the least-squares\n"
00706       "function.  This allows you to add a parameter regularization term to any existing\n"
00707       "response function!"
00708       );
00709     pl->set( ObservationTargetAsParameter_name_, ObservationTargetAsParameter_default_,
00710       "If true, then a parameter will be accepted for the state observation vector\n"
00711       "to allow it to be set by an external client through the InArgs object."
00712       );
00713     pl->set( ParameterMultiplier_name_,ParameterMultiplier_default_,
00714       "parameterMultiplier" );
00715     if(this->get_parameterBaseIO().get())
00716       parameterBaseReader_.set_fileIO(this->get_parameterBaseIO());
00717     pl->sublist(ParameterBaseVector_name_).setParameters(
00718       *parameterBaseReader_.getValidParameters()
00719       );
00720     this->setLocalVerbosityLevelValidatedParameter(&*pl);
00721     Teuchos::setupVerboseObjectSublist(&*pl);
00722     validParamList_ = pl;
00723   }
00724   return validParamList_;
00725 }
00726 
00727 
00728 // Overridden from ModelEvaulator.
00729 
00730 
00731 template<class Scalar>
00732 RCP<const VectorSpaceBase<Scalar> >
00733 DefaultInverseModelEvaluator<Scalar>::get_p_space(int l) const
00734 {
00735   if (prototypeInArgs_.Np()==0)
00736     initializeInArgsOutArgs();
00737   if ( l == prototypeInArgs_.Np()-1 && usingObservationTargetAsParameter_ )
00738     return get_obs_space();
00739   return this->getUnderlyingModel()->get_p_space(l);
00740 }
00741 
00742 
00743 template<class Scalar>
00744 RCP<const VectorSpaceBase<Scalar> >
00745 DefaultInverseModelEvaluator<Scalar>::get_g_space(int j) const
00746 {
00747   if (prototypeOutArgs_.Np()==0)
00748     initializeInArgsOutArgs();
00749   if (j==prototypeOutArgs_.Ng()-1)
00750     return inv_g_space_;
00751   return this->getUnderlyingModel()->get_g_space(j);
00752 }
00753 
00754 
00755 template<class Scalar>
00756 ModelEvaluatorBase::InArgs<Scalar>
00757 DefaultInverseModelEvaluator<Scalar>::createInArgs() const
00758 {
00759   if (prototypeInArgs_.Np()==0)
00760     initializeInArgsOutArgs();
00761   return prototypeInArgs_;
00762 }
00763 
00764 
00765 // Public functions overridden from Teuchos::Describable
00766 
00767 
00768 template<class Scalar>
00769 std::string DefaultInverseModelEvaluator<Scalar>::description() const
00770 {
00771   const RCP<const ModelEvaluator<Scalar> >
00772     thyraModel = this->getUnderlyingModel();
00773   std::ostringstream oss;
00774   oss << "Thyra::DefaultInverseModelEvaluator{";
00775   oss << "thyraModel=";
00776   if(thyraModel.get())
00777     oss << "\'"<<thyraModel->description()<<"\'";
00778   else
00779     oss << "NULL";
00780   oss << "}";
00781   return oss.str();
00782 }
00783 
00784 
00785 // Private functions overridden from ModelEvaulatorDefaultBase
00786 
00787 
00788 template<class Scalar>
00789 ModelEvaluatorBase::OutArgs<Scalar>
00790 DefaultInverseModelEvaluator<Scalar>::createOutArgsImpl() const
00791 {
00792   if (prototypeOutArgs_.Np()==0)
00793     initializeInArgsOutArgs();
00794   return prototypeOutArgs_;
00795 }
00796 
00797 
00798 template<class Scalar>
00799 void DefaultInverseModelEvaluator<Scalar>::evalModelImpl(
00800   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00801   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00802   ) const
00803 {
00804 
00805   using std::endl;
00806   using Teuchos::rcp;
00807   using Teuchos::rcp_const_cast;
00808   using Teuchos::rcp_dynamic_cast;
00809   using Teuchos::OSTab;
00810   typedef Teuchos::ScalarTraits<Scalar>  ST;
00811   typedef typename ST::magnitudeType ScalarMag;
00812   typedef ModelEvaluatorBase MEB;
00813 
00814   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_LOCALVERBLEVEL_BEGIN(
00815     "Thyra::DefaultInverseModelEvaluator",inArgs,outArgs,localVerbLevel_
00816     );
00817 
00818   const bool trace = out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_LOW);
00819   const bool print_p = out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_MEDIUM);
00820   const bool print_x = out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_EXTREME);
00821   const bool print_o = print_x;
00822 
00823   //
00824   // A) See what needs to be computed
00825   //
00826 
00827   VectorBase<Scalar>
00828     *g_inv_out = outArgs.get_g(outArgs.Ng()-1).get();
00829   MultiVectorBase<Scalar>
00830     *DgDx_inv_trans_out = get_mv(
00831       outArgs.get_DgDx(outArgs.Ng()-1),"DgDx",MEB::DERIV_TRANS_MV_BY_ROW
00832       ).get();
00833   MultiVectorBase<Scalar>
00834     *DgDp_inv_trans_out = get_mv(
00835       outArgs.get_DgDp(outArgs.Ng()-1,p_idx_),"DgDp",MEB::DERIV_TRANS_MV_BY_ROW
00836       ).get();
00837   const bool computeInverseFunction = ( g_inv_out || DgDx_inv_trans_out || DgDp_inv_trans_out );
00838   
00839   //
00840   // B) Compute all of the needed functions from the base model
00841   //
00842 
00843   if(trace)
00844     *out << "\nComputing the base point and the observation(s) ...\n";
00845 
00846   MEB::InArgs<Scalar>  wrappedInArgs = thyraModel->createInArgs();
00847   wrappedInArgs.setArgs(inArgs,true);
00848   MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
00849   wrappedOutArgs.setArgs(outArgs,true);
00850   RCP<VectorBase<Scalar> > wrapped_o;
00851   MEB::Derivative<Scalar> wrapped_DoDx;
00852   MEB::Derivative<Scalar> wrapped_DoDp_trans;
00853   if( obs_idx_ >= 0 && computeInverseFunction )
00854   {
00855     wrapped_o = createMember(thyraModel->get_g_space(obs_idx_));
00856     wrappedOutArgs.set_g(obs_idx_,wrapped_o);
00857     if( DgDx_inv_trans_out ) {
00858       if (!observationPassThrough_)
00859         wrapped_DoDx = thyraModel->create_DgDx_op(obs_idx_);
00860       else
00861         wrapped_DoDx = Thyra::create_DgDx_mv(
00862           *thyraModel, obs_idx_, MEB::DERIV_TRANS_MV_BY_ROW );
00863       wrappedOutArgs.set_DgDx(obs_idx_,wrapped_DoDx);
00864     }
00865     if( DgDp_inv_trans_out ) {
00866       wrapped_DoDp_trans = create_DgDp_mv(
00867         *thyraModel, obs_idx_, p_idx_, MEB::DERIV_TRANS_MV_BY_ROW
00868         );
00869       wrappedOutArgs.set_DgDp(obs_idx_,p_idx_,wrapped_DoDp_trans);
00870     }
00871     // 2007/07/28: rabartl: Above, we really should check if these output
00872     // arguments have already been set by the client.  If they are, then we
00873     // need to make sure that they are of the correct form or we need to throw
00874     // an exception!
00875   }
00876 
00877   if (!wrappedOutArgs.isEmpty()) {
00878     thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
00879   }
00880   else {
00881     if(trace)
00882       *out << "\nSkipping the evaluation of the underlying model since "
00883            << "there is nothing to compute ...\n";
00884   }
00885   
00886   bool failed = wrappedOutArgs.isFailed();
00887 
00888   //
00889   // C) Assemble the final observation and paramter terms
00890   //
00891   
00892   if ( !failed && computeInverseFunction ) {
00893 
00894     //
00895     // Compute the inverse response function and its derivatives
00896     //
00897 
00898     RCP<const VectorBase<Scalar> >
00899       x_in = inArgs.get_x(),
00900       p_in = inArgs.get_p(p_idx_);
00901 
00902     const MEB::InArgs<Scalar> nominalValues = this->getNominalValues();
00903     RCP<const VectorBase<Scalar> >
00904       x = ( !is_null(x_in) ? x_in : nominalValues.get_x().assert_not_null() ),
00905       p = ( !is_null(p_in) ? p_in : nominalValues.get_p(p_idx_).assert_not_null() );
00906 
00907     const RCP<const VectorSpaceBase<Scalar> >
00908       o_space = get_obs_space(),
00909       p_space = this->get_p_space(p_idx_);
00910 
00911     const Ordinal
00912       no = o_space->dim(),
00913       np = p_space->dim();
00914     
00915     if (trace)
00916       *out << "\nno = " << no
00917            << "\nnp = " << np
00918            << endl;
00919 
00920 #ifdef TEUCHOS_DEBUG
00921     TEST_FOR_EXCEPTION(
00922       observationPassThrough_ && no != 1, std::logic_error,
00923       "Error, the observation function dimension no="<<no<<" > 1 is not allowed"
00924       " when the observation is passed through as the observation matching term!"
00925       );
00926 #endif
00927 
00928     // Compute diff_o if needed
00929     RCP<const VectorBase<Scalar> > o;
00930     RCP<VectorBase<Scalar> > diff_o;
00931     if( !observationPassThrough_ && ( g_inv_out || DgDx_inv_trans_out )  ) {
00932       if (obs_idx_ < 0 ) o = x; else o = wrapped_o; // can't use ( test ? x : wrapped_o )!
00933       if(trace) *out << "\n||o||inf = " << norm_inf(*o) << endl;
00934       if (print_o) *out << "\no = " << *o;
00935       diff_o = createMember(o_space);
00936       RCP<const VectorBase<Scalar> >
00937         observationTarget
00938         = ( observationTargetAsParameter_
00939           ? inArgs.get_p(inArgs.Np()-1)
00940           : Teuchos::null
00941           );
00942       if (is_null(observationTarget) ) {
00943         observationTarget = observationTarget_;
00944         if (trace)
00945           *out << "\n||ot||inf = " << norm_inf(*observationTarget) << endl;
00946         if (print_o)
00947           *out << "\not = " << *observationTarget;
00948       }
00949       if (!is_null(observationTarget)) {
00950         V_VmV( &*diff_o, *o, *observationTarget );
00951       }
00952       else {
00953         assign( &*diff_o, *o );
00954       }
00955       if(trace)
00956         *out << "\n||diff_o||inf = " << norm_inf(*diff_o) << endl;
00957       if (print_o)
00958         *out << "\ndiff_o = " << *diff_o;
00959     }
00960   
00961     // Compute diff_p if needed
00962     RCP<VectorBase<Scalar> > diff_p;
00963     if( g_inv_out || DgDp_inv_trans_out ) {
00964       if(trace) *out << "\n||p||inf = " << norm_inf(*p) << endl;
00965       if(print_p) *out << "\np = " << Teuchos::describe(*p,Teuchos::VERB_EXTREME);
00966       diff_p = createMember(p_space);
00967       if (!is_null(parameterBase_) ) {
00968         if(trace) *out << "\n||pt||inf = " << norm_inf(*parameterBase_) << endl;
00969         if(print_p) *out << "\npt = " << Teuchos::describe(*parameterBase_,Teuchos::VERB_EXTREME);
00970         V_VmV( &*diff_p, *p, *parameterBase_ );
00971       }
00972       else {
00973         assign( &*diff_p, *p );
00974       }
00975       if(trace) *out << "\n||diff_p|| = " << norm(*diff_p) << endl;
00976       if(print_p) *out << "\ndiff_p = " << Teuchos::describe(*diff_p,Teuchos::VERB_EXTREME);
00977     }
00978     
00979 
00980     // Get and check Q_o and Q_p
00981 
00982     RCP<const LinearOpBase<Scalar> >
00983       Q_o = this->get_observationMatchWeightingOp(),
00984       Q_p = this->get_parameterRegularizationWeightingOp();
00985 
00986 #ifdef TEUCHOS_DEBUG
00987     if (!is_null(Q_o)) {
00988       THYRA_ASSERT_VEC_SPACES(
00989         "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
00990         *Q_o->range(), *o_space
00991         );
00992       THYRA_ASSERT_VEC_SPACES(
00993         "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
00994         *Q_o->domain(), *o_space
00995         );
00996     }
00997     if (!is_null(Q_p)) {
00998       THYRA_ASSERT_VEC_SPACES(
00999         "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
01000         *Q_p->range(), *p_space
01001         );
01002       THYRA_ASSERT_VEC_SPACES(
01003         "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
01004         *Q_p->domain(), *p_space
01005         );
01006     }
01007     // Note, we have not proved that Q_o and Q_p are s.p.d. but at least we
01008     // have established that that have the right range and domain spaces!
01009 #endif
01010 
01011     // Compute Q_o * diff_o
01012     RCP<VectorBase<Scalar> > Q_o_diff_o;
01013     if ( !is_null(Q_o) && !is_null(diff_o) ) {
01014       Q_o_diff_o = createMember(Q_o->range()); // Should be same as domain!
01015       apply( *Q_o, NOTRANS, *diff_o, &*Q_o_diff_o );
01016     }
01017     
01018     // Compute Q_p * diff_p
01019     RCP<VectorBase<Scalar> > Q_p_diff_p;
01020     if ( !is_null(Q_p)  && !is_null(diff_p)  ) {
01021       Q_p_diff_p = createMember(Q_p->range()); // Should be same as domain!
01022       apply( *Q_p, NOTRANS, *diff_p, &*Q_p_diff_p );
01023     }
01024 
01025     // Compute g_inv(x,p)
01026     if(g_inv_out) {
01027       if(trace)
01028         *out << "\nComputing inverse response function ginv = g(Np-1) ...\n";
01029       const Scalar observationTerm
01030         = ( observationPassThrough_
01031           ? get_ele(*wrapped_o,0) // ToDo; Verify that this is already a scalar
01032           : ( observationMultiplier_ != ST::zero()
01033             ? ( !is_null(Q_o)
01034               ?  observationMultiplier_*0.5*dot(*diff_o,*Q_o_diff_o)
01035               : observationMultiplier_*(0.5/no)*dot(*diff_o,*diff_o)
01036               )
01037             : ST::zero()
01038             )
01039           );
01040       const Scalar parameterTerm
01041         = ( parameterMultiplier_ != ST::zero()
01042           ? ( !is_null(Q_p)
01043             ?  parameterMultiplier_*0.5*dot(*diff_p,*Q_p_diff_p)
01044             : parameterMultiplier_*(0.5/np)*dot(*diff_p,*diff_p)
01045             )
01046           : ST::zero()
01047           );
01048       const Scalar g_inv_val = observationTerm+parameterTerm;
01049       if(trace)
01050         *out
01051           << "\nObservation matching term of ginv = g(Np-1):"
01052           << "\n  observationMultiplier = " << observationMultiplier_
01053           << "\n  observationMultiplier*observationMatch(x,p) = " << observationTerm
01054           << "\nParameter regularization term of ginv = g(Np-1):"
01055           << "\n  parameterMultiplier = " << parameterMultiplier_
01056           << "\n  parameterMultiplier*parameterRegularization(p) = " << parameterTerm
01057           << "\nginv = " << g_inv_val
01058           << "\n";
01059       set_ele(0,observationTerm+parameterTerm,g_inv_out);
01060     }
01061 
01062     // Compute d(g_inv)/d(x)^T
01063     if(DgDx_inv_trans_out) {
01064       if(trace)
01065         *out << "\nComputing inverse response function derivative DginvDx^T:\n";
01066       if (!observationPassThrough_) {
01067         if( obs_idx_ < 0 ) {
01068           if (!is_null(Q_o)) {
01069             if (trace)
01070               *out << "\nDginvDx^T = observationMultiplier * Q_o * diff_o ...\n";
01071             V_StV(
01072               &*DgDx_inv_trans_out->col(0),
01073               observationMultiplier_,
01074               *Q_o_diff_o
01075               );
01076           }
01077           else {
01078             if (trace)
01079               *out << "\nDginvDx^T = observationMultiplier * (1/no) * diff_o ...\n";
01080             V_StV(
01081               &*DgDx_inv_trans_out->col(0),
01082               Scalar(observationMultiplier_*(1.0/no)),
01083               *diff_o
01084               );
01085           }
01086         }
01087         else {
01088           //if (trace)
01089           //  *out << "\n||DoDx^T||inf = " << norms_inf(*wrapped_DoDx.getMultiVector()) << endl;
01090           if (print_o && print_x)
01091             *out << "\nDoDx = " << *wrapped_DoDx.getLinearOp();
01092           if (!is_null(Q_o)) {
01093             if (trace)
01094               *out << "\nDginvDx^T = observationMultiplier * DoDx^T * Q_o * diff_o ...\n";
01095             apply(
01096               *wrapped_DoDx.getLinearOp(), CONJTRANS,
01097               *Q_o_diff_o,
01098               &*DgDx_inv_trans_out->col(0),
01099               observationMultiplier_
01100               );
01101           }
01102           else {
01103             if (trace)
01104               *out << "\nDginvDx^T = (observationMultiplier*(1/no)) * DoDx^T * diff_o ...\n";
01105             apply(
01106               *wrapped_DoDx.getLinearOp(), CONJTRANS,
01107               *diff_o,
01108               &*DgDx_inv_trans_out->col(0),
01109               Scalar(observationMultiplier_*(1.0/no))
01110               );
01111           }
01112         }
01113       }
01114       else {
01115         if (trace)
01116           *out << "\nDginvDx^T = observationMultiplier * DoDx^T ...\n";
01117         V_StV(
01118           &*DgDx_inv_trans_out->col(0), observationMultiplier_,
01119           *wrapped_DoDx.getMultiVector()->col(0)
01120           );
01121       }
01122       if(trace)
01123         *out << "\n||DginvDx^T||inf = " << norms_inf(*DgDx_inv_trans_out) << "\n";
01124       if (print_x)
01125         *out << "\nDginvDx^T = " << *DgDx_inv_trans_out;
01126     }
01127 
01128     // Compute d(g_inv)/d(p)^T
01129     if(DgDp_inv_trans_out) {
01130       if(trace)
01131         *out << "\nComputing inverse response function derivative DginvDp^T ...\n";
01132       if (obs_idx_ >= 0) {
01133         if (trace)
01134           *out << "\n||DoDp^T|| = " << norms_inf(*wrapped_DoDp_trans.getMultiVector()) << endl;
01135         if (print_p)
01136           *out << "\nDoDp^T = " << Teuchos::describe(*wrapped_DoDp_trans.getMultiVector(),Teuchos::VERB_EXTREME);
01137       }
01138       if(trace)
01139         *out << "\nDginvDp^T = 0 ...\n";
01140       assign( &*DgDp_inv_trans_out->col(0), ST::zero() );
01141       // DgDp^T += observationMultiplier * d(observationMatch)/d(p)^T
01142       if (!observationPassThrough_) {
01143         if ( obs_idx_ >= 0 ) {
01144           if ( !is_null(Q_o) ) {
01145             if(trace)
01146               *out << "\nDginvDp^T += observationMultiplier* * (DoDp^T) * Q_o * diff_o ...\n";
01147             apply(
01148               *wrapped_DoDp_trans.getMultiVector(), NOTRANS,
01149               *Q_o_diff_o,
01150               &*DgDp_inv_trans_out->col(0),
01151               Scalar(observationMultiplier_*(1.0/no)),
01152               ST::one()
01153               );
01154           }
01155           else {
01156             if(trace)
01157               *out << "\nDgDp^T += observationMultiplier* * (DoDp^T) * Q_o * diff_o ...\n";
01158             apply(
01159               *wrapped_DoDp_trans.getMultiVector(), NOTRANS,
01160               *diff_o,
01161               &*DgDp_inv_trans_out->col(0),
01162               Scalar(observationMultiplier_*(1.0/no)),
01163               ST::one()
01164               );
01165           }
01166           if(trace)
01167             *out << "\n||DginvDp^T||inf = " << norms_inf(*DgDp_inv_trans_out) << "\n";
01168           if (print_p)
01169             *out << "\nDginvDp^T = " << *DgDp_inv_trans_out;
01170         }
01171         else {
01172           // d(observationMatch)/d(p)^T = 0, nothing to do!
01173         }
01174       }
01175       else {
01176         if(trace)
01177           *out << "\nDginvDp^T += (observationMultiplier*(1/no)) * (DoDp^T) * diff_o ...\n";
01178         Vp_StV(
01179           &*DgDp_inv_trans_out->col(0), observationMultiplier_,
01180           *wrapped_DoDp_trans.getMultiVector()->col(0)
01181           );
01182         
01183       }
01184       // DgDp^T += parameterMultiplier * d(parameterRegularization)/d(p)^T
01185       if( parameterMultiplier_ != ST::zero() ) {
01186         if ( !is_null(Q_p) ) {
01187           if(trace)
01188             *out << "\nDginvDp^T += parameterMultiplier * Q_p * diff_p ...\n";
01189           Vp_StV(
01190             &*DgDp_inv_trans_out->col(0),
01191             parameterMultiplier_,
01192             *Q_p_diff_p
01193             );
01194         }
01195         else {
01196           if(trace)
01197             *out << "\nDginvDp^T += (parameterMultiplier*(1.0/np)) * diff_p ...\n";
01198           Vp_StV(
01199             &*DgDp_inv_trans_out->col(0),
01200             Scalar(parameterMultiplier_*(1.0/np)),
01201             *diff_p
01202             );
01203         }
01204         if(trace)
01205           *out << "\n||DginvDp^T||inf = " << norms_inf(*DgDp_inv_trans_out) << "\n";
01206         if (print_p)
01207         *out << "\nDginvDp^T = " << *DgDp_inv_trans_out;
01208       }
01209       else {
01210         // This term is zero so there is nothing to do!
01211       }
01212     }
01213 
01214   }
01215   
01216   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
01217   
01218 }
01219 
01220 
01221 // private
01222 
01223 
01224 template<class Scalar>
01225 void DefaultInverseModelEvaluator<Scalar>::initializeDefaults()
01226 {
01227   obs_idx_ = ObservationIndex_default_;
01228   p_idx_ = ParameterSubvectorIndex_default_;
01229   observationMultiplier_ = ObservationMultiplier_default_;
01230   parameterMultiplier_ = ParameterMultiplier_default_; 
01231 }
01232 
01233 
01234 template<class Scalar>
01235 void DefaultInverseModelEvaluator<Scalar>::initializeInArgsOutArgs() const
01236 {
01237 
01238   typedef ModelEvaluatorBase MEB;
01239 
01240   const RCP<const ModelEvaluator<Scalar> >
01241     thyraModel = this->getUnderlyingModel();
01242 
01243   const MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
01244   const int wrapped_Np = wrappedInArgs.Np();
01245 
01246   MEB::InArgsSetup<Scalar> inArgs;
01247   inArgs.setModelEvalDescription(this->description());
01248   const bool supports_x = wrappedInArgs.supports(MEB::IN_ARG_x);
01249   usingObservationTargetAsParameter_ = ( supports_x && observationTargetAsParameter_ );
01250   inArgs.setSupports(
01251     wrappedInArgs,
01252     wrapped_Np + ( usingObservationTargetAsParameter_ ? 1 : 0 )
01253     );
01254   prototypeInArgs_ = inArgs;
01255 
01256   const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
01257   const int wrapped_Ng = wrappedOutArgs.Ng();
01258 
01259   MEB::OutArgsSetup<Scalar> outArgs;
01260   outArgs.setModelEvalDescription(inArgs.modelEvalDescription());
01261   outArgs.set_Np_Ng( prototypeInArgs_.Np(), wrapped_Ng+1 );
01262   outArgs.setSupports(wrappedOutArgs);
01263   outArgs.setSupports(MEB::OUT_ARG_DgDx,wrapped_Ng,MEB::DERIV_TRANS_MV_BY_ROW);
01264   outArgs.setSupports(MEB::OUT_ARG_DgDp,wrapped_Ng,p_idx_,MEB::DERIV_TRANS_MV_BY_ROW);
01265   prototypeOutArgs_ = outArgs;
01266   
01267 }
01268 
01269 
01270 template<class Scalar>
01271 RCP<const VectorSpaceBase<Scalar> >
01272 DefaultInverseModelEvaluator<Scalar>::get_obs_space() const
01273 {
01274   return ( obs_idx_ < 0 ? this->get_x_space() : this->get_g_space(obs_idx_) );
01275 }
01276 
01277 
01278 } // namespace Thyra
01279 
01280 
01281 #endif // THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines