Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultLumpedParameterModelEvaluator.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_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
00030 #define THYRA_DEFAUL_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
00031 
00032 
00033 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00034 #include "Thyra_ModelEvaluatorHelpers.hpp"
00035 #include "Thyra_DetachedVectorView.hpp"
00036 #include "Teuchos_ParameterListAcceptor.hpp"
00037 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00038 #include "Teuchos_Time.hpp"
00039 #include "Teuchos_Assert.hpp"
00040 #include "Teuchos_as.hpp"
00041 
00042 #include "sillyModifiedGramSchmidt.hpp" // This is just an example!
00043 
00044 
00045 namespace Thyra {
00046 
00047 
00226 template<class Scalar>
00227 class DefaultLumpedParameterModelEvaluator
00228   : virtual public ModelEvaluatorDelegatorBase<Scalar>
00229   , virtual public Teuchos::ParameterListAcceptor
00230 {
00231 public:
00232 
00235 
00237   DefaultLumpedParameterModelEvaluator();
00238 
00240   void initialize(
00241     const RCP<ModelEvaluator<Scalar> > &thyraModel
00242     );
00243 
00245   void uninitialize(
00246     RCP<ModelEvaluator<Scalar> > *thyraModel
00247     );
00248 
00249   // 2007/07/30: rabartl: ToDo: Add functions to set and get the underlying
00250   // basis matrix!
00251 
00253 
00256 
00258   std::string description() const;
00259 
00261 
00264 
00266   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00268   RCP<Teuchos::ParameterList> getNonconstParameterList();
00270   RCP<Teuchos::ParameterList> unsetParameterList();
00272   RCP<const Teuchos::ParameterList> getParameterList() const;
00274   RCP<const Teuchos::ParameterList> getValidParameters() const;
00275 
00277 
00280 
00282   RCP<const VectorSpaceBase<Scalar> > get_p_space(int l) const;
00284   RCP<const Array<std::string> > get_p_names(int l) const;
00286   ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00288   ModelEvaluatorBase::InArgs<Scalar> getLowerBounds() const;
00290   ModelEvaluatorBase::InArgs<Scalar> getUpperBounds() const;
00292   void reportFinalPoint(
00293     const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
00294     const bool wasSolved
00295     );
00296 
00298 
00299 private:
00300 
00303 
00305   ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00307   void evalModelImpl(
00308     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00309     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00310     ) const;
00311 
00313 
00314 private:
00315 
00316   // ////////////////////////////////
00317   // Private data members
00318 
00319   mutable bool isInitialized_;
00320   mutable bool nominalValuesAndBoundsUpdated_;
00321 
00322   mutable RCP<const Teuchos::ParameterList> validParamList_;
00323   RCP<Teuchos::ParameterList> paramList_;
00324 
00325   // Parameters read from the parameter list
00326   int p_idx_;
00327   bool autogenerateBasisMatrix_;
00328   int numberOfBasisColumns_;
00329   bool nominalValueIsParameterBase_;
00330   bool ignoreParameterBounds_;
00331   Teuchos::EVerbosityLevel localVerbLevel_;
00332   bool dumpBasisMatrix_;
00333 
00334   // Reduced affine parameter model
00335   mutable RCP<const MultiVectorBase<Scalar> > B_;
00336   mutable RCP<const VectorBase<Scalar> > p_orig_base_;
00337 
00338   // Nominal values and bounds
00339   mutable ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00340   mutable ModelEvaluatorBase::InArgs<Scalar> lowerBounds_;
00341   mutable ModelEvaluatorBase::InArgs<Scalar> upperBounds_;
00342 
00343   // Static
00344 
00345   static const std::string ParameterSubvectorIndex_name_;
00346   static const int ParameterSubvectorIndex_default_;
00347 
00348   static const std::string AutogenerateBasisMatrix_name_;
00349   static const bool AutogenerateBasisMatrix_default_;
00350 
00351   static const std::string NumberOfBasisColumns_name_;
00352   static const int NumberOfBasisColumns_default_;
00353 
00354   static const std::string NominalValueIsParameterBase_name_;
00355   static const bool NominalValueIsParameterBase_default_;
00356 
00357   static const std::string ParameterBaseVector_name_;
00358 
00359   static const std::string IgnoreParameterBounds_name_;
00360   static const bool IgnoreParameterBounds_default_;
00361 
00362   static const std::string DumpBasisMatrix_name_;
00363   static const bool DumpBasisMatrix_default_;
00364 
00365   // ////////////////////////////////
00366   // Private member functions
00367 
00368   // These functions are used to implement late initialization so that the
00369   // need for clients to order function calls is reduced.
00370 
00371   // Finish enough initialization to defined spaces etc.
00372   void finishInitialization() const;
00373 
00374   // Generate the parameter basis matrix B.
00375   void generateParameterBasisMatrix() const;
00376 
00377   // Finish all of initialization needed to define nominal values, bounds, and
00378   // p_orig_base.  This calls finishInitialization().
00379   void updateNominalValuesAndBounds() const;
00380 
00381   // Map from p -> p_orig.
00382   RCP<VectorBase<Scalar> >
00383   map_from_p_to_p_orig( const VectorBase<Scalar> &p ) const;
00384 
00385   // Set up the arguments for DhDp_orig to be computed by the underlying model.
00386   void setupWrappedParamDerivOutArgs(
00387     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs, // in
00388     ModelEvaluatorBase::OutArgs<Scalar> *wrappedOutArgs // in/out
00389     ) const;
00390 
00391   // Create DhDp_orig needed to assembled DhDp
00392   ModelEvaluatorBase::Derivative<Scalar>
00393   create_deriv_wrt_p_orig(
00394     const ModelEvaluatorBase::Derivative<Scalar> &DhDp,
00395     const ModelEvaluatorBase::EDerivativeMultiVectorOrientation requiredOrientation
00396     ) const;
00397 
00398   // After DhDp_orig has been computed, assemble DhDp or DhDp^T for all deriv
00399   // output arguments.
00400   void assembleParamDerivOutArgs(
00401     const ModelEvaluatorBase::OutArgs<Scalar> &wrappedOutArgs, // in
00402     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs // in/out
00403     ) const;
00404 
00405   // Given a single DhDp_orig, assemble DhDp
00406   void assembleParamDeriv(
00407     const ModelEvaluatorBase::Derivative<Scalar> &DhDp_orig, // in
00408     const ModelEvaluatorBase::Derivative<Scalar> &DhDp // in/out
00409     ) const;
00410 
00411 };
00412 
00413 
00418 template<class Scalar>
00419 RCP<DefaultLumpedParameterModelEvaluator<Scalar> >
00420 defaultLumpedParameterModelEvaluator(
00421   const RCP<ModelEvaluator<Scalar> > &thyraModel
00422   )
00423 {
00424   RCP<DefaultLumpedParameterModelEvaluator<Scalar> >
00425     paramLumpedModel = Teuchos::rcp(new DefaultLumpedParameterModelEvaluator<Scalar>);
00426   paramLumpedModel->initialize(thyraModel);
00427   return paramLumpedModel;
00428 }
00429 
00430 
00431 // /////////////////////////////////
00432 // Implementations
00433 
00434 
00435 // Static data members
00436 
00437 
00438 template<class Scalar>
00439 const std::string
00440 DefaultLumpedParameterModelEvaluator<Scalar>::ParameterSubvectorIndex_name_
00441 = "Parameter Subvector Index";
00442 
00443 template<class Scalar>
00444 const int
00445 DefaultLumpedParameterModelEvaluator<Scalar>::ParameterSubvectorIndex_default_
00446 = 0;
00447 
00448 
00449 template<class Scalar>
00450 const std::string
00451 DefaultLumpedParameterModelEvaluator<Scalar>::AutogenerateBasisMatrix_name_
00452 = "Auto-generate Basis Matrix";
00453 
00454 template<class Scalar>
00455 const bool
00456 DefaultLumpedParameterModelEvaluator<Scalar>::AutogenerateBasisMatrix_default_
00457 = true;
00458 
00459 
00460 template<class Scalar>
00461 const std::string
00462 DefaultLumpedParameterModelEvaluator<Scalar>::NumberOfBasisColumns_name_
00463 = "Number of Basis Columns";
00464 
00465 template<class Scalar>
00466 const int
00467 DefaultLumpedParameterModelEvaluator<Scalar>::NumberOfBasisColumns_default_
00468 = 1;
00469 
00470 
00471 template<class Scalar>
00472 const std::string
00473 DefaultLumpedParameterModelEvaluator<Scalar>::NominalValueIsParameterBase_name_
00474 = "Nominal Value is Parameter Base";
00475 
00476 template<class Scalar>
00477 const bool
00478 DefaultLumpedParameterModelEvaluator<Scalar>::NominalValueIsParameterBase_default_
00479 = true;
00480 
00481 
00482 template<class Scalar>
00483 const std::string
00484 DefaultLumpedParameterModelEvaluator<Scalar>::ParameterBaseVector_name_
00485 = "Parameter Base Vector";
00486 
00487 
00488 template<class Scalar>
00489 const std::string
00490 DefaultLumpedParameterModelEvaluator<Scalar>::IgnoreParameterBounds_name_
00491 = "Ignore Parameter Bounds";
00492 
00493 template<class Scalar>
00494 const bool
00495 DefaultLumpedParameterModelEvaluator<Scalar>::IgnoreParameterBounds_default_
00496 = false;
00497 
00498 
00499 template<class Scalar>
00500 const std::string
00501 DefaultLumpedParameterModelEvaluator<Scalar>::DumpBasisMatrix_name_
00502 = "Dump Basis Matrix";
00503 
00504 template<class Scalar>
00505 const bool
00506 DefaultLumpedParameterModelEvaluator<Scalar>::DumpBasisMatrix_default_
00507 = false;
00508 
00509 
00510 // Constructors/initializers/accessors/utilities
00511 
00512 
00513 template<class Scalar>
00514 DefaultLumpedParameterModelEvaluator<Scalar>::DefaultLumpedParameterModelEvaluator()
00515   :isInitialized_(false),
00516    nominalValuesAndBoundsUpdated_(false),
00517    p_idx_(ParameterSubvectorIndex_default_),
00518    autogenerateBasisMatrix_(AutogenerateBasisMatrix_default_),
00519    numberOfBasisColumns_(NumberOfBasisColumns_default_),
00520    nominalValueIsParameterBase_(NominalValueIsParameterBase_default_),
00521    ignoreParameterBounds_(IgnoreParameterBounds_default_),
00522    localVerbLevel_(Teuchos::VERB_DEFAULT),
00523    dumpBasisMatrix_(DumpBasisMatrix_default_)
00524 {}
00525 
00526 
00527 template<class Scalar>
00528 void DefaultLumpedParameterModelEvaluator<Scalar>::initialize(
00529   const RCP<ModelEvaluator<Scalar> > &thyraModel
00530   )
00531 {
00532   isInitialized_ = false;
00533   nominalValuesAndBoundsUpdated_ = false;
00534   this->ModelEvaluatorDelegatorBase<Scalar>::initialize(thyraModel);
00535 }
00536 
00537 
00538 template<class Scalar>
00539 void DefaultLumpedParameterModelEvaluator<Scalar>::uninitialize(
00540   RCP<ModelEvaluator<Scalar> > *thyraModel
00541   )
00542 {
00543   isInitialized_ = false;
00544   if(thyraModel) *thyraModel = this->getUnderlyingModel();
00545   this->ModelEvaluatorDelegatorBase<Scalar>::uninitialize();
00546 }
00547 
00548 
00549 // Public functions overridden from Teuchos::Describable
00550 
00551 
00552 template<class Scalar>
00553 std::string
00554 DefaultLumpedParameterModelEvaluator<Scalar>::description() const
00555 {
00556   const RCP<const ModelEvaluator<Scalar> >
00557     thyraModel = this->getUnderlyingModel();
00558   std::ostringstream oss;
00559   oss << "Thyra::DefaultLumpedParameterModelEvaluator{";
00560   oss << "thyraModel=";
00561   if(thyraModel.get())
00562     oss << "\'"<<thyraModel->description()<<"\'";
00563   else
00564     oss << "NULL";
00565   oss << "}";
00566   return oss.str();
00567 }
00568 
00569 
00570 // Overridden from Teuchos::ParameterListAcceptor
00571 
00572 
00573 template<class Scalar>
00574 void DefaultLumpedParameterModelEvaluator<Scalar>::setParameterList(
00575   RCP<Teuchos::ParameterList> const& paramList
00576   )
00577 {
00578 
00579   using Teuchos::getParameterPtr;
00580   using Teuchos::rcp;
00581   using Teuchos::sublist;
00582 
00583   isInitialized_ = false;
00584   nominalValuesAndBoundsUpdated_ = false;
00585 
00586   // Validate and set the parameter list
00587   TEST_FOR_EXCEPT(is_null(paramList));
00588   paramList->validateParameters(*getValidParameters(),0);
00589   paramList_ = paramList;
00590 
00591   // Read in parameters
00592   p_idx_ = paramList_->get(
00593     ParameterSubvectorIndex_name_, ParameterSubvectorIndex_default_ );
00594   autogenerateBasisMatrix_ = paramList_->get(
00595     AutogenerateBasisMatrix_name_, AutogenerateBasisMatrix_default_ );
00596   if (autogenerateBasisMatrix_) {
00597     numberOfBasisColumns_ = paramList_->get(
00598       NumberOfBasisColumns_name_, NumberOfBasisColumns_default_ );
00599   }
00600   nominalValueIsParameterBase_ = paramList_->get(
00601     NominalValueIsParameterBase_name_, NominalValueIsParameterBase_default_ );
00602   if (!nominalValueIsParameterBase_) {
00603     TEST_FOR_EXCEPT("ToDo: Implement reading parameter base vector from file!");
00604   }
00605   ignoreParameterBounds_ = paramList_->get(
00606     IgnoreParameterBounds_name_, IgnoreParameterBounds_default_ );
00607   dumpBasisMatrix_ = paramList_->get(
00608     DumpBasisMatrix_name_, DumpBasisMatrix_default_ );
00609 
00610   // Verbosity settings
00611   localVerbLevel_ = this->readLocalVerbosityLevelValidatedParameter(*paramList_);
00612   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00613 
00614 #ifdef TEUCHOS_DEBUG
00615   paramList_->validateParameters(*getValidParameters(),0);
00616 #endif
00617 
00618 }
00619 
00620 
00621 template<class Scalar>
00622 RCP<Teuchos::ParameterList>
00623 DefaultLumpedParameterModelEvaluator<Scalar>::getNonconstParameterList()
00624 {
00625   return paramList_;
00626 }
00627 
00628 
00629 template<class Scalar>
00630 RCP<Teuchos::ParameterList>
00631 DefaultLumpedParameterModelEvaluator<Scalar>::unsetParameterList()
00632 {
00633   RCP<Teuchos::ParameterList> _paramList = paramList_;
00634   paramList_ = Teuchos::null;
00635   return _paramList;
00636 }
00637 
00638 
00639 template<class Scalar>
00640 RCP<const Teuchos::ParameterList>
00641 DefaultLumpedParameterModelEvaluator<Scalar>::getParameterList() const
00642 {
00643   return paramList_;
00644 }
00645 
00646 
00647 template<class Scalar>
00648 RCP<const Teuchos::ParameterList>
00649 DefaultLumpedParameterModelEvaluator<Scalar>::getValidParameters() const
00650 {
00651   if(validParamList_.get()==NULL) {
00652     RCP<Teuchos::ParameterList>
00653       pl = Teuchos::rcp(new Teuchos::ParameterList());
00654     pl->set( ParameterSubvectorIndex_name_, ParameterSubvectorIndex_default_,
00655       "Determines the index of the parameter subvector in the underlying model\n"
00656       "for which the reduced basis representation will be determined." );
00657     pl->set( AutogenerateBasisMatrix_name_, AutogenerateBasisMatrix_default_,
00658       "If true, then a basis matrix will be auto-generated for a given number\n"
00659       " of basis vectors." );
00660     pl->set( NumberOfBasisColumns_name_, NumberOfBasisColumns_default_,
00661       "If a basis is auto-generated, then this parameter gives the number\n"
00662       "of columns in the basis matrix that will be created.  Warning!  This\n"
00663       "number must be less than or equal to the number of original parameters\n"
00664       "or an exception will be thrown!" );
00665     pl->set( NominalValueIsParameterBase_name_, NominalValueIsParameterBase_default_,
00666       "If true, then the nominal values for the full parameter subvector from the\n"
00667       "underlying model will be used for p_orig_base.  This allows p==0 to give\n"
00668       "the nominal values for the parameters." );
00669     /*
00670     if(this->get_parameterBaseIO().get())
00671       parameterBaseReader_.set_fileIO(this->get_parameterBaseIO());
00672     pl->sublist(ParameterBaseVector_name_).setParameters(
00673       *parameterBaseReader_.getValidParameters()
00674       );
00675     */
00676     pl->set( IgnoreParameterBounds_name_, IgnoreParameterBounds_default_,
00677       "If true, then any bounds on the parameter subvector will be ignored." );
00678     pl->set( DumpBasisMatrix_name_, DumpBasisMatrix_default_,
00679       "If true, then the basis matrix will be printed the first time it is created\n"
00680       "as part of the verbose output and as part of the Describable::describe(...)\n"
00681       "output for any verbositiy level >= \"low\"." );
00682     this->setLocalVerbosityLevelValidatedParameter(&*pl);
00683     Teuchos::setupVerboseObjectSublist(&*pl);
00684     validParamList_ = pl;
00685   }
00686   return validParamList_;
00687 }
00688 
00689 
00690 // Overridden from ModelEvaulator.
00691 
00692 
00693 template<class Scalar>
00694 RCP<const VectorSpaceBase<Scalar> >
00695 DefaultLumpedParameterModelEvaluator<Scalar>::get_p_space(int l) const
00696 {
00697   finishInitialization();
00698   if (l == p_idx_)
00699     return B_->domain();
00700   return this->getUnderlyingModel()->get_p_space(l);
00701 }
00702 
00703 
00704 template<class Scalar>
00705 RCP<const Array<std::string> >
00706 DefaultLumpedParameterModelEvaluator<Scalar>::get_p_names(int l) const
00707 {
00708   finishInitialization();
00709   if (l == p_idx_)
00710     return Teuchos::null; // Names for these parameters would be meaningless!
00711   return this->getUnderlyingModel()->get_p_names(l);
00712 }
00713 
00714 
00715 template<class Scalar>
00716 ModelEvaluatorBase::InArgs<Scalar>
00717 DefaultLumpedParameterModelEvaluator<Scalar>::getNominalValues() const
00718 {
00719   updateNominalValuesAndBounds();
00720   return nominalValues_;
00721 }
00722 
00723 
00724 template<class Scalar>
00725 ModelEvaluatorBase::InArgs<Scalar>
00726 DefaultLumpedParameterModelEvaluator<Scalar>::getLowerBounds() const
00727 {
00728   updateNominalValuesAndBounds();
00729   return lowerBounds_;
00730 }
00731 
00732 
00733 template<class Scalar>
00734 ModelEvaluatorBase::InArgs<Scalar>
00735 DefaultLumpedParameterModelEvaluator<Scalar>::getUpperBounds() const
00736 {
00737   updateNominalValuesAndBounds();
00738   return upperBounds_;
00739 }
00740 
00741 
00742 template<class Scalar>
00743 void DefaultLumpedParameterModelEvaluator<Scalar>::reportFinalPoint(
00744   const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
00745   const bool wasSolved
00746   )
00747 {
00748 
00749   typedef ModelEvaluatorBase MEB;
00750 
00751   // Make sure that everything has been initialized
00752   updateNominalValuesAndBounds();
00753 
00754   const RCP<ModelEvaluator<Scalar> >
00755     thyraModel = this->getNonconstUnderlyingModel();
00756 
00757   // By default, copy all input arguments since they will all be the same
00758   // except for the given reduced p.  We will then replace the reduced p with
00759   // p_orig below.
00760   MEB::InArgs<Scalar> wrappedFinalPoint = thyraModel->createInArgs();
00761   wrappedFinalPoint.setArgs(finalPoint);
00762 
00763   // Replace p with p_orig.
00764   RCP<const VectorBase<Scalar> > p;
00765   if (!is_null(p=finalPoint.get_p(p_idx_))) {
00766     wrappedFinalPoint.set_p(p_idx_,map_from_p_to_p_orig(*p));
00767   }
00768 
00769   thyraModel->reportFinalPoint(wrappedFinalPoint,wasSolved);
00770 
00771 }
00772 
00773 
00774 // Private functions overridden from ModelEvaulatorDefaultBase
00775 
00776 
00777 template<class Scalar>
00778 ModelEvaluatorBase::OutArgs<Scalar>
00779 DefaultLumpedParameterModelEvaluator<Scalar>::createOutArgsImpl() const
00780 {
00781   ModelEvaluatorBase::OutArgsSetup<Scalar>
00782     outArgs = this->getUnderlyingModel()->createOutArgs();
00783   outArgs.setModelEvalDescription(this->description());
00784   return outArgs;
00785   // 2007/07/31: rabartl: ToDo: We need to manually set the forms of the
00786   // derivatives that this class object will support!  This needs to be based
00787   // on tests of what the forms of derivatives the underlying model supports.
00788 }
00789 
00790 
00791 template<class Scalar>
00792 void DefaultLumpedParameterModelEvaluator<Scalar>::evalModelImpl(
00793   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00794   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00795   ) const
00796 {
00797 
00798   // This routine is pretty simple for the most part.  By default, we just
00799   // pass everything through to the underlying model evaluator except for
00800   // arguments reated to the parameter subvector with index
00801   // p_idx_.
00802 
00803   using Teuchos::rcp;
00804   using Teuchos::rcp_const_cast;
00805   using Teuchos::rcp_dynamic_cast;
00806   using Teuchos::OSTab;
00807   typedef Teuchos::ScalarTraits<Scalar>  ST;
00808   typedef typename ST::magnitudeType ScalarMag;
00809   typedef ModelEvaluatorBase MEB;
00810 
00811   // Make sure that everything has been initialized
00812   updateNominalValuesAndBounds();
00813 
00814   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_LOCALVERBLEVEL_BEGIN(
00815     "Thyra::DefaultLumpedParameterModelEvaluator",inArgs,outArgs,localVerbLevel_
00816     );
00817 
00818   //
00819   // A) Setup InArgs
00820   //
00821 
00822   // By default, copy all input arguments since they will all be the same
00823   // except for the given reduced p.  We will then replace the reduced p with
00824   // p_orig below.
00825   MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
00826   wrappedInArgs.setArgs(inArgs);
00827   
00828   // Replace p with p_orig.
00829   RCP<const VectorBase<Scalar> > p;
00830   if (!is_null(p=wrappedInArgs.get_p(p_idx_))) {
00831     if (
00832       dumpBasisMatrix_
00833       && includesVerbLevel(localVerbLevel,Teuchos::VERB_MEDIUM)
00834       )
00835     {
00836       *out << "\nB = " << Teuchos::describe(*B_,Teuchos::VERB_EXTREME);
00837     }
00838     wrappedInArgs.set_p(p_idx_,map_from_p_to_p_orig(*p));
00839   }
00840 
00841   //
00842   // B) Setup OutArgs
00843   //
00844 
00845   // By default, copy all output arguments since they will all be the same
00846   // except for those derivatives w.r.t. p(p_idx).  We will then replace the
00847   // derivative objects w.r.t. given reduced p with the derivarive objects
00848   // w.r.t. p_orig below.
00849   MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
00850   wrappedOutArgs.setArgs(outArgs);
00851 
00852   // Set derivative output arguments for p_orig if derivatives for p are
00853   // reqeusted in outArgs
00854   setupWrappedParamDerivOutArgs(outArgs,&wrappedOutArgs);
00855 
00856   //
00857   // C) Evaluate the underlying model functions
00858   //
00859 
00860   if (includesVerbLevel(localVerbLevel,Teuchos::VERB_LOW))
00861     *out << "\nEvaluating the fully parameterized underlying model ...\n";
00862   // Compute the underlying functions in terms of p_orig, including
00863   // derivatives w.r.t. p_orig.
00864   thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
00865 
00866   //
00867   // D) Postprocess the output arguments
00868   //
00869 
00870   // Assemble the derivatives for p given derivatives for p_orig computed
00871   // above.
00872   assembleParamDerivOutArgs(wrappedOutArgs,outArgs);
00873   
00874   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00875   
00876 }
00877 
00878 
00879 // private
00880 
00881 
00882 template<class Scalar>
00883 void DefaultLumpedParameterModelEvaluator<Scalar>::finishInitialization() const
00884 {
00885 
00886   typedef ScalarTraits<Scalar> ST;
00887   typedef ModelEvaluatorBase MEB;
00888 
00889   if (isInitialized_)
00890     return;
00891 
00892   //
00893   // A) Get the underlying model
00894   //
00895 
00896   const RCP<const ModelEvaluator<Scalar> >
00897     thyraModel = this->getUnderlyingModel();
00898 
00899   TEST_FOR_EXCEPTION(
00900     is_null(thyraModel), std::logic_error,
00901     "Error, the underlying model evaluator must be set!" );
00902 
00903   //
00904   // B) Create B for the reduced affine model for the given parameter subvector
00905   //
00906 
00907   if (autogenerateBasisMatrix_) {
00908     generateParameterBasisMatrix();
00909   }
00910   else {
00911     TEST_FOR_EXCEPTION(
00912       true, std::logic_error,
00913       "Error, we don't handle a client-set parameter basis matrix yet!" );
00914   }
00915 
00916   isInitialized_ = true;
00917 
00918 }
00919 
00920 
00921 template<class Scalar>
00922 void DefaultLumpedParameterModelEvaluator<Scalar>::generateParameterBasisMatrix() const
00923 {
00924 
00925   using Teuchos::as;
00926   typedef ScalarTraits<Scalar> ST;
00927 
00928   const RCP<const ModelEvaluator<Scalar> >
00929     thyraModel = this->getUnderlyingModel();
00930 
00931   const RCP<const VectorSpaceBase<Scalar> >
00932     p_orig_space = thyraModel->get_p_space(p_idx_);
00933 
00934   const Ordinal p_orig_dim = p_orig_space->dim();
00935 
00936   TEST_FOR_EXCEPTION(
00937     !( 1 <= numberOfBasisColumns_ && numberOfBasisColumns_ <= p_orig_dim ),
00938     std::logic_error,
00939     "Error, the number of basis columns = " << numberOfBasisColumns_ << " does not\n"
00940     "fall in the range [1,"<<p_orig_dim<<"]!" );
00941 
00942   //
00943   // Create and randomize B
00944   //
00945   // Here we make the first column all ones and then randomize columns 1
00946   // through numberOfBasisColumns_-1 so that the average entry is 1.0 with a
00947   // spread of 1.0.  This is just to give as well a scaled starting matrix as
00948   // possible that will hopefully yeild a well scaled orthonomal B after we
00949   // are finished.
00950 
00951   RCP<MultiVectorBase<Scalar> >
00952     B = createMembers(p_orig_space,numberOfBasisColumns_);
00953   assign( &*B->col(0), ST::one() );
00954   if (numberOfBasisColumns_ > 1) {
00955     seed_randomize<double>(0);
00956     Thyra::randomize( as<Scalar>(0.5*ST::one()), as<Scalar>(1.5*ST::one()), &*B );
00957   }
00958 
00959   //
00960   // Create an orthogonal form of B using a modified Gram Schmidt algorithm
00961   //
00962 
00963   RCP<MultiVectorBase<double> > R;
00964   sillyModifiedGramSchmidt( &*B, &R );
00965 
00966   // Above:
00967   // 1) On output, B will have orthonomal columns which makes it a good basis
00968   // 2) We just discard the "R" factor since we don't need it for anything 
00969 
00970   B_ = B;
00971 
00972 }
00973 
00974 
00975 template<class Scalar>
00976 void DefaultLumpedParameterModelEvaluator<Scalar>::updateNominalValuesAndBounds() const
00977 {
00978 
00979   typedef ScalarTraits<Scalar> ST;
00980   typedef ModelEvaluatorBase MEB;
00981 
00982   if (nominalValuesAndBoundsUpdated_)
00983     return;
00984 
00985   finishInitialization();
00986 
00987   const RCP<const ModelEvaluator<Scalar> >
00988     thyraModel = this->getUnderlyingModel();
00989 
00990   const MEB::InArgs<Scalar> origNominalValues = thyraModel->getNominalValues();
00991   const MEB::InArgs<Scalar> origLowerBounds = thyraModel->getLowerBounds();
00992   const MEB::InArgs<Scalar> origUpperBounds = thyraModel->getUpperBounds();
00993 
00994   // p_orig_base
00995 
00996   if (nominalValueIsParameterBase_) {
00997     const RCP<const VectorBase<Scalar> >
00998       p_orig_init = origNominalValues.get_p(p_idx_);
00999     TEST_FOR_EXCEPTION(
01000       is_null(p_orig_init), std::logic_error,
01001       "Error, if the user requested that the nominal values be used\n"
01002       "as the base vector p_orig_base then that vector has to exist!" );
01003     p_orig_base_ = p_orig_init->clone_v();
01004   }
01005   else {
01006     TEST_FOR_EXCEPTION(
01007       true, std::logic_error,
01008       "Error, we don't handle reading in the parameter base vector yet!" );
01009   }
01010 
01011   // Nominal values
01012 
01013   nominalValues_ = origNominalValues;
01014   
01015   if (nominalValueIsParameterBase_) {
01016     // A value of p==0 will give p_orig = p_orig_init!
01017     const RCP<VectorBase<Scalar> >
01018       p_init = createMember(B_->domain());
01019     assign( &*p_init, ST::zero() );
01020     nominalValues_.set_p(p_idx_,p_init);
01021   }
01022   else {
01023     TEST_FOR_EXCEPTION(
01024       true, std::logic_error,
01025       "Error, we don't handle creating p_init when p_orig_base != p_orig_init yet!" );
01026   }
01027 
01028   // Bounds
01029 
01030   lowerBounds_ = origLowerBounds;
01031   upperBounds_ = origUpperBounds;
01032 
01033   lowerBounds_.set_p(p_idx_,Teuchos::null);
01034   upperBounds_.set_p(p_idx_,Teuchos::null);
01035 
01036   if (!ignoreParameterBounds_) {
01037     const RCP<const VectorBase<Scalar> >
01038       p_orig_l = origLowerBounds.get_p(p_idx_),
01039       p_orig_u = origUpperBounds.get_p(p_idx_);
01040     if ( !is_null(p_orig_l) || !is_null(p_orig_u) ) {
01041       TEST_FOR_EXCEPTION(
01042         true, std::logic_error,
01043         "Error, we don't handle bounds on p_orig yet!" );
01044     }
01045   }
01046 
01047   nominalValuesAndBoundsUpdated_ = true;
01048 
01049 }
01050 
01051 
01052 template<class Scalar>
01053 RCP<VectorBase<Scalar> >
01054 DefaultLumpedParameterModelEvaluator<Scalar>::map_from_p_to_p_orig(
01055   const VectorBase<Scalar> &p
01056   ) const
01057 {
01058   // p_orig = B*p + p_orig_base
01059   RCP<VectorBase<Scalar> > p_orig = createMember(B_->range());
01060   apply( *B_, NOTRANS, p, &*p_orig );
01061   Vp_V( &*p_orig, *p_orig_base_ );
01062   return p_orig;
01063 }
01064 
01065 
01066 template<class Scalar>
01067 void DefaultLumpedParameterModelEvaluator<Scalar>::setupWrappedParamDerivOutArgs(
01068   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs, // in
01069   ModelEvaluatorBase::OutArgs<Scalar> *wrappedOutArgs_inout // in/out
01070   ) const
01071 {
01072 
01073   typedef ModelEvaluatorBase MEB;
01074   typedef MEB::Derivative<Scalar> Deriv;
01075 
01076   TEST_FOR_EXCEPT(wrappedOutArgs_inout==0);
01077   MEB::OutArgs<Scalar> &wrappedOutArgs = *wrappedOutArgs_inout;
01078     
01079   Deriv DfDp;
01080   if ( !(DfDp=outArgs.get_DfDp(p_idx_)).isEmpty() ) {
01081     wrappedOutArgs.set_DfDp(p_idx_,create_deriv_wrt_p_orig(DfDp,MEB::DERIV_MV_BY_COL));
01082   }
01083 
01084   const int Ng = outArgs.Ng();
01085   for ( int j = 0; j < Ng; ++j ) {
01086     Deriv DgDp;
01087     if ( !(DgDp=outArgs.get_DgDp(j,p_idx_)).isEmpty() ) {
01088       wrappedOutArgs.set_DgDp(
01089         j, p_idx_,
01090         create_deriv_wrt_p_orig(DgDp,DgDp.getMultiVectorOrientation())
01091         );
01092     }
01093   }
01094 
01095 }
01096 
01097 
01098 template<class Scalar>
01099 ModelEvaluatorBase::Derivative<Scalar>
01100 DefaultLumpedParameterModelEvaluator<Scalar>::create_deriv_wrt_p_orig(
01101   const ModelEvaluatorBase::Derivative<Scalar> &DhDp,
01102   const ModelEvaluatorBase::EDerivativeMultiVectorOrientation requiredOrientation
01103   ) const
01104 {
01105 
01106   typedef ModelEvaluatorBase MEB;
01107 
01108   const RCP<const MultiVectorBase<Scalar> >
01109     DhDp_mv = DhDp.getMultiVector();
01110   TEST_FOR_EXCEPTION(
01111     is_null(DhDp_mv) || (DhDp.getMultiVectorOrientation() != requiredOrientation),
01112     std::logic_error,
01113     "Error, we currently can't handle non-multi-vector derivatives!" );
01114 
01115   RCP<MultiVectorBase<Scalar> > DhDp_orig_mv;
01116   switch (requiredOrientation) {
01117     case MEB::DERIV_MV_BY_COL:
01118       // DhDp = DhDp_orig * B
01119       DhDp_orig_mv = createMembers(DhDp_mv->range(),B_->range()->dim());
01120       // Above, we could just request DhDp_orig as a LinearOpBase object since
01121       // we just need to apply it!
01122       break;
01123     case MEB::DERIV_TRANS_MV_BY_ROW:
01124       // (DhDp^T) = B^T * (DhDp_orig^T)  [DhDp_orig_mv is transposed!]
01125       DhDp_orig_mv = createMembers(B_->range(),DhDp_mv->domain()->dim());
01126       // Above, we really do need DhDp_orig as the gradient form multi-vector
01127       // since it must be the RHS for a linear operator apply!
01128       break;
01129     default:
01130       TEST_FOR_EXCEPT(true); // Should never get here!
01131   }
01132   
01133   return MEB::Derivative<Scalar>(DhDp_orig_mv,requiredOrientation);
01134   
01135 }
01136 
01137 
01138 template<class Scalar>
01139 void DefaultLumpedParameterModelEvaluator<Scalar>::assembleParamDerivOutArgs(
01140   const ModelEvaluatorBase::OutArgs<Scalar> &wrappedOutArgs, // in
01141   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs // in/out
01142   ) const
01143 {
01144 
01145   typedef ModelEvaluatorBase MEB;
01146   typedef MEB::Derivative<Scalar> Deriv;
01147     
01148   Deriv DfDp;
01149   if ( !(DfDp=outArgs.get_DfDp(p_idx_)).isEmpty() ) {
01150     assembleParamDeriv( wrappedOutArgs.get_DfDp(p_idx_), DfDp );
01151   }
01152 
01153   const int Ng = outArgs.Ng();
01154   for ( int j = 0; j < Ng; ++j ) {
01155     Deriv DgDp;
01156     if ( !(DgDp=outArgs.get_DgDp(j,p_idx_)).isEmpty() ) {
01157       assembleParamDeriv( wrappedOutArgs.get_DgDp(j,p_idx_), DgDp );
01158     }
01159   }
01160 
01161 }
01162 
01163 
01164 template<class Scalar>
01165 void DefaultLumpedParameterModelEvaluator<Scalar>::assembleParamDeriv(
01166   const ModelEvaluatorBase::Derivative<Scalar> &DhDp_orig, // in
01167   const ModelEvaluatorBase::Derivative<Scalar> &DhDp // in/out
01168   ) const
01169 {
01170 
01171   typedef ModelEvaluatorBase MEB;
01172 
01173   const RCP<const MultiVectorBase<Scalar> >
01174     DhDp_orig_mv = DhDp_orig.getMultiVector();
01175   TEST_FOR_EXCEPTION(
01176     is_null(DhDp_orig_mv), std::logic_error,
01177     "Error, we currently can't handle non-multi-vector derivatives!" );
01178 
01179   const RCP<MultiVectorBase<Scalar> >
01180     DhDp_mv = DhDp.getMultiVector();
01181   TEST_FOR_EXCEPTION(
01182     is_null(DhDp_mv), std::logic_error,
01183     "Error, we currently can't handle non-multi-vector derivatives!" );
01184 
01185   switch( DhDp_orig.getMultiVectorOrientation() ) {
01186     case MEB::DERIV_MV_BY_COL:
01187       // DhDp = DhDp_orig * B
01188 #ifdef TEUCHSO_DEBUG
01189       TEUCHOS_ASSERT(
01190         DhDp.getMultiVectorOrientation() == MEB::DERIV_MV_BY_COL );
01191 #endif
01192       apply( *DhDp_orig_mv, NOTRANS, *B_, &*DhDp_mv );
01193       // Above, we could generalize DhDp_oirg to just be a general linear
01194       // operator.
01195       break;
01196     case MEB::DERIV_TRANS_MV_BY_ROW:
01197       // (DhDp^T) = B^T * (DhDp_orig^T)  [DhDp_orig_mv is transposed!]
01198 #ifdef TEUCHSO_DEBUG
01199       TEUCHOS_ASSERT(
01200         DhDp.getMultiVectorOrientation() == MEB::DERIV_TRANS_MV_BY_ROW );
01201 #endif
01202       apply( *B_, CONJTRANS, *DhDp_orig_mv, &*DhDp_mv );
01203       break;
01204     default:
01205       TEST_FOR_EXCEPT(true); // Should never get here!
01206   }
01207 
01208 }
01209 
01210 
01211 } // namespace Thyra
01212 
01213 
01214 #endif // THYRA_DEFAUL_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines