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