42#ifndef THYRA_DEFAUL_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
43#define THYRA_DEFAUL_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
46#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
47#include "Thyra_ModelEvaluatorHelpers.hpp"
48#include "Thyra_DetachedVectorView.hpp"
49#include "Teuchos_ParameterListAcceptor.hpp"
50#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
51#include "Teuchos_Time.hpp"
52#include "Teuchos_Assert.hpp"
53#include "Teuchos_as.hpp"
55#include "sillyModifiedGramSchmidt.hpp"
239template<
class Scalar>
332 mutable bool isInitialized_;
333 mutable bool nominalValuesAndBoundsUpdated_;
340 bool autogenerateBasisMatrix_;
341 int numberOfBasisColumns_;
342 bool nominalValueIsParameterBase_;
343 bool ignoreParameterBounds_;
345 bool dumpBasisMatrix_;
358 static const std::string ParameterSubvectorIndex_name_;
359 static const int ParameterSubvectorIndex_default_;
361 static const std::string AutogenerateBasisMatrix_name_;
362 static const bool AutogenerateBasisMatrix_default_;
364 static const std::string NumberOfBasisColumns_name_;
365 static const int NumberOfBasisColumns_default_;
367 static const std::string NominalValueIsParameterBase_name_;
368 static const bool NominalValueIsParameterBase_default_;
370 static const std::string ParameterBaseVector_name_;
372 static const std::string IgnoreParameterBounds_name_;
373 static const bool IgnoreParameterBounds_default_;
375 static const std::string DumpBasisMatrix_name_;
376 static const bool DumpBasisMatrix_default_;
385 void finishInitialization()
const;
388 void generateParameterBasisMatrix()
const;
392 void updateNominalValuesAndBounds()
const;
399 void setupWrappedParamDerivOutArgs(
406 create_deriv_wrt_p_orig(
413 void assembleParamDerivOutArgs(
419 void assembleParamDeriv(
431template<
class Scalar>
439 paramLumpedModel->initialize(thyraModel);
440 return paramLumpedModel;
451template<
class Scalar>
454=
"Parameter Subvector Index";
456template<
class Scalar>
462template<
class Scalar>
465=
"Auto-generate Basis Matrix";
467template<
class Scalar>
473template<
class Scalar>
476=
"Number of Basis Columns";
478template<
class Scalar>
484template<
class Scalar>
487=
"Nominal Value is Parameter Base";
489template<
class Scalar>
495template<
class Scalar>
498=
"Parameter Base Vector";
501template<
class Scalar>
504=
"Ignore Parameter Bounds";
506template<
class Scalar>
512template<
class Scalar>
515=
"Dump Basis Matrix";
517template<
class Scalar>
526template<
class Scalar>
528 :isInitialized_(false),
529 nominalValuesAndBoundsUpdated_(false),
530 p_idx_(ParameterSubvectorIndex_default_),
531 autogenerateBasisMatrix_(AutogenerateBasisMatrix_default_),
532 numberOfBasisColumns_(NumberOfBasisColumns_default_),
533 nominalValueIsParameterBase_(NominalValueIsParameterBase_default_),
534 ignoreParameterBounds_(IgnoreParameterBounds_default_),
536 dumpBasisMatrix_(DumpBasisMatrix_default_)
540template<
class Scalar>
545 isInitialized_ =
false;
546 nominalValuesAndBoundsUpdated_ =
false;
551template<
class Scalar>
556 isInitialized_ =
false;
557 if(thyraModel) *thyraModel = this->getUnderlyingModel();
565template<
class Scalar>
570 thyraModel = this->getUnderlyingModel();
571 std::ostringstream oss;
572 oss <<
"Thyra::DefaultLumpedParameterModelEvaluator{";
573 oss <<
"thyraModel=";
575 oss <<
"\'"<<thyraModel->description()<<
"\'";
586template<
class Scalar>
592 using Teuchos::getParameterPtr;
594 using Teuchos::sublist;
596 isInitialized_ =
false;
597 nominalValuesAndBoundsUpdated_ =
false;
601 paramList->validateParameters(*getValidParameters(),0);
602 paramList_ = paramList;
605 p_idx_ = paramList_->
get(
606 ParameterSubvectorIndex_name_, ParameterSubvectorIndex_default_ );
607 autogenerateBasisMatrix_ = paramList_->get(
608 AutogenerateBasisMatrix_name_, AutogenerateBasisMatrix_default_ );
609 if (autogenerateBasisMatrix_) {
610 numberOfBasisColumns_ = paramList_->get(
611 NumberOfBasisColumns_name_, NumberOfBasisColumns_default_ );
613 nominalValueIsParameterBase_ = paramList_->get(
614 NominalValueIsParameterBase_name_, NominalValueIsParameterBase_default_ );
615 if (!nominalValueIsParameterBase_) {
618 ignoreParameterBounds_ = paramList_->get(
619 IgnoreParameterBounds_name_, IgnoreParameterBounds_default_ );
620 dumpBasisMatrix_ = paramList_->get(
621 DumpBasisMatrix_name_, DumpBasisMatrix_default_ );
624 localVerbLevel_ = this->readLocalVerbosityLevelValidatedParameter(*paramList_);
625 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
628 paramList_->validateParameters(*getValidParameters(),0);
634template<
class Scalar>
642template<
class Scalar>
647 paramList_ = Teuchos::null;
652template<
class Scalar>
660template<
class Scalar>
664 if(validParamList_.get()==NULL) {
667 pl->set( ParameterSubvectorIndex_name_, ParameterSubvectorIndex_default_,
668 "Determines the index of the parameter subvector in the underlying model\n"
669 "for which the reduced basis representation will be determined." );
670 pl->set( AutogenerateBasisMatrix_name_, AutogenerateBasisMatrix_default_,
671 "If true, then a basis matrix will be auto-generated for a given number\n"
672 " of basis vectors." );
673 pl->set( NumberOfBasisColumns_name_, NumberOfBasisColumns_default_,
674 "If a basis is auto-generated, then this parameter gives the number\n"
675 "of columns in the basis matrix that will be created. Warning! This\n"
676 "number must be less than or equal to the number of original parameters\n"
677 "or an exception will be thrown!" );
678 pl->set( NominalValueIsParameterBase_name_, NominalValueIsParameterBase_default_,
679 "If true, then the nominal values for the full parameter subvector from the\n"
680 "underlying model will be used for p_orig_base. This allows p==0 to give\n"
681 "the nominal values for the parameters." );
689 pl->set( IgnoreParameterBounds_name_, IgnoreParameterBounds_default_,
690 "If true, then any bounds on the parameter subvector will be ignored." );
691 pl->set( DumpBasisMatrix_name_, DumpBasisMatrix_default_,
692 "If true, then the basis matrix will be printed the first time it is created\n"
693 "as part of the verbose output and as part of the Describable::describe(...)\n"
694 "output for any verbositiy level >= \"low\"." );
695 this->setLocalVerbosityLevelValidatedParameter(&*pl);
696 Teuchos::setupVerboseObjectSublist(&*pl);
697 validParamList_ = pl;
699 return validParamList_;
706template<
class Scalar>
710 finishInitialization();
713 return this->getUnderlyingModel()->get_p_space(l);
717template<
class Scalar>
721 finishInitialization();
723 return Teuchos::null;
724 return this->getUnderlyingModel()->get_p_names(l);
728template<
class Scalar>
732 updateNominalValuesAndBounds();
733 return nominalValues_;
737template<
class Scalar>
741 updateNominalValuesAndBounds();
746template<
class Scalar>
750 updateNominalValuesAndBounds();
755template<
class Scalar>
765 updateNominalValuesAndBounds();
768 thyraModel = this->getNonconstUnderlyingModel();
773 MEB::InArgs<Scalar> wrappedFinalPoint = thyraModel->createInArgs();
774 wrappedFinalPoint.setArgs(finalPoint);
778 if (!is_null(p=finalPoint.
get_p(p_idx_))) {
779 wrappedFinalPoint.set_p(p_idx_, map_from_p_to_p_orig(*p));
782 thyraModel->reportFinalPoint(wrappedFinalPoint,wasSolved);
790template<
class Scalar>
795 outArgs = this->getUnderlyingModel()->createOutArgs();
804template<
class Scalar>
805void DefaultLumpedParameterModelEvaluator<Scalar>::evalModelImpl(
806 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
807 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
817 using Teuchos::rcp_const_cast;
818 using Teuchos::rcp_dynamic_cast;
821 typedef typename ST::magnitudeType ScalarMag;
822 typedef ModelEvaluatorBase MEB;
825 updateNominalValuesAndBounds();
827 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_LOCALVERBLEVEL_BEGIN(
828 "Thyra::DefaultLumpedParameterModelEvaluator",inArgs,outArgs,localVerbLevel_
838 MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
839 wrappedInArgs.setArgs(inArgs);
842 RCP<const VectorBase<Scalar> > p;
843 if (!is_null(p=wrappedInArgs.get_p(p_idx_))) {
851 wrappedInArgs.set_p(p_idx_,map_from_p_to_p_orig(*p));
862 MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
863 wrappedOutArgs.setArgs(outArgs);
867 setupWrappedParamDerivOutArgs(outArgs,&wrappedOutArgs);
874 *out <<
"\nEvaluating the fully parameterized underlying model ...\n";
877 thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
885 assembleParamDerivOutArgs(wrappedOutArgs,outArgs);
887 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
895template<
class Scalar>
896void DefaultLumpedParameterModelEvaluator<Scalar>::finishInitialization()
const
899 typedef ScalarTraits<Scalar> ST;
900 typedef ModelEvaluatorBase MEB;
909 const RCP<const ModelEvaluator<Scalar> >
910 thyraModel = this->getUnderlyingModel();
913 is_null(thyraModel), std::logic_error,
914 "Error, the underlying model evaluator must be set!" );
920 if (autogenerateBasisMatrix_) {
921 generateParameterBasisMatrix();
925 true, std::logic_error,
926 "Error, we don't handle a client-set parameter basis matrix yet!" );
929 isInitialized_ =
true;
934template<
class Scalar>
935void DefaultLumpedParameterModelEvaluator<Scalar>::generateParameterBasisMatrix()
const
939 typedef ScalarTraits<Scalar> ST;
941 const RCP<const ModelEvaluator<Scalar> >
942 thyraModel = this->getUnderlyingModel();
944 const RCP<const VectorSpaceBase<Scalar> >
945 p_orig_space = thyraModel->get_p_space(p_idx_);
947 const Ordinal p_orig_dim = p_orig_space->dim();
950 !( 1 <= numberOfBasisColumns_ && numberOfBasisColumns_ <= p_orig_dim ),
952 "Error, the number of basis columns = " << numberOfBasisColumns_ <<
" does not\n"
953 "fall in the range [1,"<<p_orig_dim<<
"]!" );
964 const RCP<MultiVectorBase<Scalar> >
965 B = createMembers(p_orig_space,numberOfBasisColumns_);
966 assign( B->col(0).ptr(), ST::one() );
967 if (numberOfBasisColumns_ > 1) {
968 seed_randomize<double>(0);
969 Thyra::randomize( as<Scalar>(0.5*ST::one()), as<Scalar>(1.5*ST::one()),
977 RCP<MultiVectorBase<double> > R;
978 sillyModifiedGramSchmidt( B.ptr(), Teuchos::outArg(R) );
989template<
class Scalar>
990void DefaultLumpedParameterModelEvaluator<Scalar>::updateNominalValuesAndBounds()
const
993 typedef ScalarTraits<Scalar> ST;
994 typedef ModelEvaluatorBase MEB;
996 if (nominalValuesAndBoundsUpdated_)
999 finishInitialization();
1001 const RCP<const ModelEvaluator<Scalar> >
1002 thyraModel = this->getUnderlyingModel();
1004 const MEB::InArgs<Scalar> origNominalValues = thyraModel->getNominalValues();
1005 const MEB::InArgs<Scalar> origLowerBounds = thyraModel->getLowerBounds();
1006 const MEB::InArgs<Scalar> origUpperBounds = thyraModel->getUpperBounds();
1010 if (nominalValueIsParameterBase_) {
1011 const RCP<const VectorBase<Scalar> >
1012 p_orig_init = origNominalValues.get_p(p_idx_);
1014 is_null(p_orig_init), std::logic_error,
1015 "Error, if the user requested that the nominal values be used\n"
1016 "as the base vector p_orig_base then that vector has to exist!" );
1017 p_orig_base_ = p_orig_init->clone_v();
1021 true, std::logic_error,
1022 "Error, we don't handle reading in the parameter base vector yet!" );
1027 nominalValues_ = origNominalValues;
1029 if (nominalValueIsParameterBase_) {
1031 const RCP<VectorBase<Scalar> >
1032 p_init = createMember(B_->domain());
1033 assign( p_init.ptr(), ST::zero() );
1034 nominalValues_.set_p(p_idx_, p_init);
1038 true, std::logic_error,
1039 "Error, we don't handle creating p_init when p_orig_base != p_orig_init yet!" );
1044 lowerBounds_ = origLowerBounds;
1045 upperBounds_ = origUpperBounds;
1047 lowerBounds_.set_p(p_idx_,Teuchos::null);
1048 upperBounds_.set_p(p_idx_,Teuchos::null);
1050 if (!ignoreParameterBounds_) {
1051 const RCP<const VectorBase<Scalar> >
1052 p_orig_l = origLowerBounds.get_p(p_idx_),
1053 p_orig_u = origUpperBounds.get_p(p_idx_);
1056 true, std::logic_error,
1057 "Error, we don't handle bounds on p_orig yet!" );
1061 nominalValuesAndBoundsUpdated_ =
true;
1066template<
class Scalar>
1067RCP<VectorBase<Scalar> >
1068DefaultLumpedParameterModelEvaluator<Scalar>::map_from_p_to_p_orig(
1069 const VectorBase<Scalar> &p
1073 const RCP<VectorBase<Scalar> > p_orig = createMember(B_->range());
1074 apply( *B_,
NOTRANS, p, p_orig.ptr() );
1075 Vp_V( p_orig.ptr(), *p_orig_base_ );
1080template<
class Scalar>
1081void DefaultLumpedParameterModelEvaluator<Scalar>::setupWrappedParamDerivOutArgs(
1082 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs,
1083 ModelEvaluatorBase::OutArgs<Scalar> *wrappedOutArgs_inout
1087 typedef ModelEvaluatorBase MEB;
1088 typedef MEB::Derivative<Scalar> Deriv;
1091 MEB::OutArgs<Scalar> &wrappedOutArgs = *wrappedOutArgs_inout;
1094 if ( !(DfDp=outArgs.get_DfDp(p_idx_)).isEmpty() ) {
1095 wrappedOutArgs.set_DfDp(p_idx_,create_deriv_wrt_p_orig(DfDp,MEB::DERIV_MV_BY_COL));
1098 const int Ng = outArgs.Ng();
1099 for (
int j = 0; j < Ng; ++j ) {
1101 if ( !(DgDp=outArgs.get_DgDp(j,p_idx_)).isEmpty() ) {
1102 wrappedOutArgs.set_DgDp(
1104 create_deriv_wrt_p_orig(DgDp,DgDp.getMultiVectorOrientation())
1112template<
class Scalar>
1113ModelEvaluatorBase::Derivative<Scalar>
1114DefaultLumpedParameterModelEvaluator<Scalar>::create_deriv_wrt_p_orig(
1115 const ModelEvaluatorBase::Derivative<Scalar> &DhDp,
1120 typedef ModelEvaluatorBase MEB;
1122 const RCP<const MultiVectorBase<Scalar> >
1123 DhDp_mv = DhDp.getMultiVector();
1125 is_null(DhDp_mv) || (DhDp.getMultiVectorOrientation() != requiredOrientation),
1127 "Error, we currently can't handle non-multi-vector derivatives!" );
1129 RCP<MultiVectorBase<Scalar> > DhDp_orig_mv;
1130 switch (requiredOrientation) {
1131 case MEB::DERIV_MV_BY_COL:
1133 DhDp_orig_mv = createMembers(DhDp_mv->range(),B_->range()->dim());
1137 case MEB::DERIV_TRANS_MV_BY_ROW:
1139 DhDp_orig_mv = createMembers(B_->range(),DhDp_mv->domain()->dim());
1147 return MEB::Derivative<Scalar>(DhDp_orig_mv,requiredOrientation);
1152template<
class Scalar>
1153void DefaultLumpedParameterModelEvaluator<Scalar>::assembleParamDerivOutArgs(
1154 const ModelEvaluatorBase::OutArgs<Scalar> &wrappedOutArgs,
1155 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
1159 typedef ModelEvaluatorBase MEB;
1160 typedef MEB::Derivative<Scalar> Deriv;
1163 if ( !(DfDp=outArgs.get_DfDp(p_idx_)).isEmpty() ) {
1164 assembleParamDeriv( wrappedOutArgs.get_DfDp(p_idx_), DfDp );
1167 const int Ng = outArgs.Ng();
1168 for (
int j = 0; j < Ng; ++j ) {
1170 if ( !(DgDp=outArgs.get_DgDp(j,p_idx_)).isEmpty() ) {
1171 assembleParamDeriv( wrappedOutArgs.get_DgDp(j,p_idx_), DgDp );
1178template<
class Scalar>
1179void DefaultLumpedParameterModelEvaluator<Scalar>::assembleParamDeriv(
1180 const ModelEvaluatorBase::Derivative<Scalar> &DhDp_orig,
1181 const ModelEvaluatorBase::Derivative<Scalar> &DhDp
1185 typedef ModelEvaluatorBase MEB;
1187 const RCP<const MultiVectorBase<Scalar> >
1188 DhDp_orig_mv = DhDp_orig.getMultiVector();
1190 is_null(DhDp_orig_mv), std::logic_error,
1191 "Error, we currently can't handle non-multi-vector derivatives!" );
1193 const RCP<MultiVectorBase<Scalar> >
1194 DhDp_mv = DhDp.getMultiVector();
1196 is_null(DhDp_mv), std::logic_error,
1197 "Error, we currently can't handle non-multi-vector derivatives!" );
1199 switch( DhDp_orig.getMultiVectorOrientation() ) {
1200 case MEB::DERIV_MV_BY_COL:
1204 DhDp.getMultiVectorOrientation() == MEB::DERIV_MV_BY_COL );
1206 apply( *DhDp_orig_mv,
NOTRANS, *B_, DhDp_mv.ptr() );
1210 case MEB::DERIV_TRANS_MV_BY_ROW:
1214 DhDp.getMultiVectorOrientation() == MEB::DERIV_TRANS_MV_BY_ROW );
1216 apply( *B_,
CONJTRANS, *DhDp_orig_mv, DhDp_mv.ptr() );
Decorator class that wraps any ModelEvaluator object and lumps parameters together using a linear bas...
RCP< const Array< std::string > > get_p_names(int l) const
RCP< Teuchos::ParameterList > getNonconstParameterList()
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
.
RCP< DefaultLumpedParameterModelEvaluator< Scalar > > defaultLumpedParameterModelEvaluator(const RCP< ModelEvaluator< Scalar > > &thyraModel)
Non-member constructor.
RCP< const Teuchos::ParameterList > getParameterList() const
RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const
ModelEvaluatorBase::InArgs< Scalar > getLowerBounds() const
ModelEvaluatorBase::InArgs< Scalar > getUpperBounds() const
RCP< Teuchos::ParameterList > unsetParameterList()
DefaultLumpedParameterModelEvaluator()
RCP< const Teuchos::ParameterList > getValidParameters() const
ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
std::string description() const
void initialize(const RCP< ModelEvaluator< Scalar > > &thyraModel)
void reportFinalPoint(const ModelEvaluatorBase::InArgs< Scalar > &finalPoint, const bool wasSolved)
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 <= l && l < this->Np().
Protected subclass of OutArgs that only ModelEvaluator subclasses can access to set up the selection ...
void setModelEvalDescription(const std::string &modelEvalDescription)
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
Base subclass for ModelEvaluator that defines some basic types.
EDerivativeMultiVectorOrientation
This is a base class that delegetes almost all function to a wrapped model evaluator object.
void uninitialize()
Uninitialize.
void initialize(const RCP< ModelEvaluator< Scalar > > &model)
Initialize given a non-const model evaluator.
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Abstract interface for finite-dimensional dense vectors.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool is_null(const std::shared_ptr< T > &p)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
@ NOTRANS
Use the non-transposed operator.
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
TypeTo as(const TypeFrom &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)