42 #ifndef THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP 43 #define THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP 46 #include "Thyra_ModelEvaluatorDelegatorBase.hpp" 47 #include "Thyra_ModelEvaluatorHelpers.hpp" 48 #include "Thyra_DetachedVectorView.hpp" 49 #include "Thyra_ParameterDrivenMultiVectorInput.hpp" 50 #include "Thyra_VectorSpaceFactoryBase.hpp" 51 #include "Thyra_MultiVectorStdOps.hpp" 52 #include "Thyra_AssertOp.hpp" 53 #include "Teuchos_StandardMemberCompositionMacros.hpp" 54 #include "Teuchos_StandardCompositionMacros.hpp" 55 #include "Teuchos_ParameterListAcceptor.hpp" 56 #include "Teuchos_VerboseObjectParameterListHelpers.hpp" 57 #include "Teuchos_StandardParameterEntryValidators.hpp" 58 #include "Teuchos_Time.hpp" 262 template<
class Scalar>
265 ,
virtual public Teuchos::ParameterListAcceptor
348 RCP<const VectorSpaceBase<Scalar> >
get_p_space(
int l)
const;
350 RCP<const VectorSpaceBase<Scalar> >
get_g_space(
int j)
const;
376 mutable RCP<const Teuchos::ParameterList> validParamList_;
377 RCP<Teuchos::ParameterList> paramList_;
379 RCP<const VectorSpaceBase<Scalar> > inv_g_space_;
383 mutable bool usingObservationTargetAsParameter_;
389 double observationMultiplier_;
390 double parameterMultiplier_;
392 bool observationTargetAsParameter_;
394 bool observationPassThrough_;
396 Teuchos::EVerbosityLevel localVerbLevel_;
401 static const std::string ObservationIndex_name_;
402 static const int ObservationIndex_default_;
404 static const std::string ParameterSubvectorIndex_name_;
405 static const int ParameterSubvectorIndex_default_;
407 static const std::string ObservationMultiplier_name_;
408 static const double ObservationMultiplier_default_;
410 static const std::string ObservationTargetVector_name_;
412 static const std::string ObservationTargetAsParameter_name_;
413 static const bool ObservationTargetAsParameter_default_;
415 static const std::string ObservationPassThrough_name_;
416 static const bool ObservationPassThrough_default_;
418 static const std::string ParameterMultiplier_name_;
419 static const double ParameterMultiplier_default_;
421 static const std::string ParameterBaseVector_name_;
426 void initializeDefaults();
428 void initializeInArgsOutArgs()
const;
430 RCP<const VectorSpaceBase<Scalar> > get_obs_space()
const;
439 template<
class Scalar>
440 RCP<DefaultInverseModelEvaluator<Scalar> >
445 RCP<DefaultInverseModelEvaluator<Scalar> >
447 inverseModel->initialize(thyraModel);
459 template<
class Scalar>
462 =
"Observation Index";
464 template<
class Scalar>
470 template<
class Scalar>
473 =
"Parameter Subvector Ordinal";
475 template<
class Scalar>
481 template<
class Scalar>
484 =
"Observation Multiplier";
486 template<
class Scalar>
492 template<
class Scalar>
495 =
"Observation Target Vector";
498 template<
class Scalar>
501 =
"Observation Target as Parameter";
503 template<
class Scalar>
509 template<
class Scalar>
512 =
"Observation Pass Through";
514 template<
class Scalar>
520 template<
class Scalar>
523 =
"Parameter Multiplier";
525 template<
class Scalar>
531 template<
class Scalar>
534 =
"Parameter Base Vector";
540 template<
class Scalar>
542 :usingObservationTargetAsParameter_(false), obs_idx_(-1),p_idx_(0),
543 observationTargetAsParameter_(false),
544 observationPassThrough_(ObservationPassThrough_default_),
545 localVerbLevel_(
Teuchos::VERB_DEFAULT),
546 observationMultiplier_(0.0),
547 parameterMultiplier_(0.0)
551 template<
class Scalar>
557 inv_g_space_= thyraModel->get_x_space()->smallVecSpcFcty()->createVecSpc(1);
564 template<
class Scalar>
569 if(thyraModel) *thyraModel = this->getUnderlyingModel();
577 template<
class Scalar>
579 RCP<Teuchos::ParameterList>
const& paramList
583 using Teuchos::Array;
584 using Teuchos::getParameterPtr;
586 using Teuchos::sublist;
589 TEUCHOS_TEST_FOR_EXCEPT(0==paramList.get());
590 paramList->validateParameters(*getValidParameters(),0);
591 paramList_ = paramList;
594 obs_idx_ = paramList_->get(
595 ObservationIndex_name_,ObservationIndex_default_);
596 observationPassThrough_ = paramList_->get(
597 ObservationPassThrough_name_, ObservationPassThrough_default_ );
599 TEUCHOS_TEST_FOR_EXCEPTION(
600 ( obs_idx_ < 0 && observationPassThrough_ ), std::logic_error,
601 "Error, the observation function index obs_idx = " << obs_idx_ <<
" is not\n" 602 "allowed when the observation is simply passed through!" 605 observationMultiplier_ = paramList_->get(
606 ObservationMultiplier_name_,ObservationMultiplier_default_);
607 if (!ObservationPassThrough_default_) {
608 observationTargetAsParameter_ = paramList_->get(
609 ObservationTargetAsParameter_name_, ObservationTargetAsParameter_default_ );
610 if(get_observationTargetIO().
get()) {
611 observationTargetReader_.set_vecSpc(get_obs_space());
612 Teuchos::VerboseObjectTempState<ParameterDrivenMultiVectorInput<Scalar> >
613 vots_observationTargetReader(
614 rcp(&observationTargetReader_,
false)
615 ,this->getOStream(),this->getVerbLevel()
617 observationTargetReader_.setParameterList(
618 sublist(paramList_,ObservationTargetVector_name_)
620 RCP<VectorBase<Scalar> >
622 observationTargetReader_.readVector(
623 "observation target vector",&observationTarget);
624 observationTarget_ = observationTarget;
628 observationTargetAsParameter_ =
false;
629 observationTarget_ = Teuchos::null;
633 p_idx_ = paramList_->get(
634 ParameterSubvectorIndex_name_,ParameterSubvectorIndex_default_);
635 parameterMultiplier_ = paramList_->get(
636 ParameterMultiplier_name_,ParameterMultiplier_default_);
637 if(get_parameterBaseIO().
get()) {
638 parameterBaseReader_.set_vecSpc(this->get_p_space(p_idx_));
639 Teuchos::VerboseObjectTempState<ParameterDrivenMultiVectorInput<Scalar> >
640 vots_parameterBaseReader(
641 rcp(¶meterBaseReader_,
false)
642 ,this->getOStream(),this->getVerbLevel()
644 parameterBaseReader_.setParameterList(
645 sublist(paramList_,ParameterBaseVector_name_)
647 RCP<VectorBase<Scalar> >
649 parameterBaseReader_.readVector(
650 "parameter base vector",¶meterBase);
651 parameterBase_ = parameterBase;
655 localVerbLevel_ = this->readLocalVerbosityLevelValidatedParameter(*paramList_);
656 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
659 paramList_->validateParameters(*getValidParameters(),0);
660 #endif // TEUCHOS_DEBUG 669 template<
class Scalar>
670 RCP<Teuchos::ParameterList>
677 template<
class Scalar>
678 RCP<Teuchos::ParameterList>
681 RCP<Teuchos::ParameterList> _paramList = paramList_;
682 paramList_ = Teuchos::null;
687 template<
class Scalar>
688 RCP<const Teuchos::ParameterList>
695 template<
class Scalar>
696 RCP<const Teuchos::ParameterList>
699 if(validParamList_.get()==NULL) {
700 RCP<Teuchos::ParameterList>
701 pl = Teuchos::rcp(
new Teuchos::ParameterList());
702 pl->set( ObservationIndex_name_,ObservationIndex_default_,
703 "The index of the observation function, obs_idx.\n" 704 "If obs_idx < 0, then the observation will be the state vector x.\n" 705 "If obs_idx >= 0, then the observation will be the response function g(obs_idx)." 707 pl->set( ParameterSubvectorIndex_name_,ParameterSubvectorIndex_default_,
708 "The index of the parameter subvector that will be used in the\n" 709 "regularization term." 711 pl->set( ObservationMultiplier_name_,ObservationMultiplier_default_,
712 "observationMultiplier" 714 if(this->get_observationTargetIO().
get())
715 observationTargetReader_.set_fileIO(this->get_observationTargetIO());
716 pl->sublist(ObservationTargetVector_name_).setParameters(
717 *observationTargetReader_.getValidParameters()
719 pl->set( ObservationPassThrough_name_, ObservationPassThrough_default_,
720 "If true, then the observation will just be used instead of the least-squares\n" 721 "function. This allows you to add a parameter regularization term to any existing\n" 724 pl->set( ObservationTargetAsParameter_name_, ObservationTargetAsParameter_default_,
725 "If true, then a parameter will be accepted for the state observation vector\n" 726 "to allow it to be set by an external client through the InArgs object." 728 pl->set( ParameterMultiplier_name_,ParameterMultiplier_default_,
729 "parameterMultiplier" );
730 if(this->get_parameterBaseIO().
get())
731 parameterBaseReader_.set_fileIO(this->get_parameterBaseIO());
732 pl->sublist(ParameterBaseVector_name_).setParameters(
733 *parameterBaseReader_.getValidParameters()
735 this->setLocalVerbosityLevelValidatedParameter(&*pl);
736 Teuchos::setupVerboseObjectSublist(&*pl);
737 validParamList_ = pl;
739 return validParamList_;
746 template<
class Scalar>
747 RCP<const VectorSpaceBase<Scalar> >
750 if (prototypeInArgs_.Np()==0)
751 initializeInArgsOutArgs();
752 if ( l == prototypeInArgs_.Np()-1 && usingObservationTargetAsParameter_ )
753 return get_obs_space();
754 return this->getUnderlyingModel()->get_p_space(l);
758 template<
class Scalar>
759 RCP<const VectorSpaceBase<Scalar> >
762 if (prototypeOutArgs_.Np()==0)
763 initializeInArgsOutArgs();
764 if (j==prototypeOutArgs_.Ng()-1)
766 return this->getUnderlyingModel()->get_g_space(j);
770 template<
class Scalar>
774 if (prototypeInArgs_.Np()==0)
775 initializeInArgsOutArgs();
776 return prototypeInArgs_;
783 template<
class Scalar>
786 const RCP<const ModelEvaluator<Scalar> >
787 thyraModel = this->getUnderlyingModel();
788 std::ostringstream oss;
789 oss <<
"Thyra::DefaultInverseModelEvaluator{";
790 oss <<
"thyraModel=";
792 oss <<
"\'"<<thyraModel->description()<<
"\'";
803 template<
class Scalar>
807 if (prototypeOutArgs_.Np()==0)
808 initializeInArgsOutArgs();
809 return prototypeOutArgs_;
813 template<
class Scalar>
814 void DefaultInverseModelEvaluator<Scalar>::evalModelImpl(
815 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
816 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
822 using Teuchos::rcp_const_cast;
823 using Teuchos::rcp_dynamic_cast;
824 using Teuchos::OSTab;
825 typedef Teuchos::ScalarTraits<Scalar> ST;
826 typedef typename ST::magnitudeType ScalarMag;
827 typedef ModelEvaluatorBase MEB;
829 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_LOCALVERBLEVEL_BEGIN(
830 "Thyra::DefaultInverseModelEvaluator",inArgs,outArgs,localVerbLevel_
833 const bool trace = out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_LOW);
834 const bool print_p = out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_MEDIUM);
835 const bool print_x = out.get() && includesVerbLevel(localVerbLevel,Teuchos::VERB_EXTREME);
836 const bool print_o = print_x;
842 const RCP<VectorBase<Scalar> >
843 g_inv_out = outArgs.get_g(outArgs.Ng()-1);
844 const RCP<MultiVectorBase<Scalar> >
845 DgDx_inv_trans_out = get_mv(
846 outArgs.get_DgDx(outArgs.Ng()-1),
"DgDx", MEB::DERIV_TRANS_MV_BY_ROW
848 const RCP<MultiVectorBase<Scalar> >
849 DgDp_inv_trans_out = get_mv(
850 outArgs.get_DgDp(outArgs.Ng()-1,p_idx_),
"DgDp", MEB::DERIV_TRANS_MV_BY_ROW
852 const bool computeInverseFunction = ( nonnull(g_inv_out)
853 || nonnull(DgDx_inv_trans_out) || nonnull(DgDp_inv_trans_out) );
860 *out <<
"\nComputing the base point and the observation(s) ...\n";
862 MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
863 wrappedInArgs.setArgs(inArgs,
true);
864 MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
865 wrappedOutArgs.setArgs(outArgs,
true);
866 RCP<VectorBase<Scalar> > wrapped_o;
867 MEB::Derivative<Scalar> wrapped_DoDx;
868 MEB::Derivative<Scalar> wrapped_DoDp_trans;
869 if( obs_idx_ >= 0 && computeInverseFunction )
871 wrapped_o = createMember(thyraModel->get_g_space(obs_idx_));
872 wrappedOutArgs.set_g(obs_idx_,wrapped_o);
873 if (nonnull(DgDx_inv_trans_out)) {
874 if (!observationPassThrough_)
875 wrapped_DoDx = thyraModel->create_DgDx_op(obs_idx_);
877 wrapped_DoDx = Thyra::create_DgDx_mv(
878 *thyraModel, obs_idx_, MEB::DERIV_TRANS_MV_BY_ROW );
879 wrappedOutArgs.set_DgDx(obs_idx_,wrapped_DoDx);
881 if (nonnull(DgDp_inv_trans_out)) {
882 wrapped_DoDp_trans = create_DgDp_mv(
883 *thyraModel, obs_idx_, p_idx_, MEB::DERIV_TRANS_MV_BY_ROW
885 wrappedOutArgs.set_DgDp(obs_idx_,p_idx_,wrapped_DoDp_trans);
893 if (!wrappedOutArgs.isEmpty()) {
894 thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
898 *out <<
"\nSkipping the evaluation of the underlying model since " 899 <<
"there is nothing to compute ...\n";
902 bool failed = wrappedOutArgs.isFailed();
908 if ( !failed && computeInverseFunction ) {
914 RCP<const VectorBase<Scalar> >
915 x_in = inArgs.get_x(),
916 p_in = inArgs.get_p(p_idx_);
918 const MEB::InArgs<Scalar> nominalValues = this->getNominalValues();
919 RCP<const VectorBase<Scalar> >
920 x = ( !is_null(x_in) ? x_in : nominalValues.get_x().assert_not_null() ),
921 p = ( !is_null(p_in) ? p_in : nominalValues.get_p(p_idx_).assert_not_null() );
923 const RCP<const VectorSpaceBase<Scalar> >
924 o_space = get_obs_space(),
925 p_space = this->get_p_space(p_idx_);
932 *out <<
"\nno = " << no
937 TEUCHOS_TEST_FOR_EXCEPTION(
938 observationPassThrough_ && no != 1, std::logic_error,
939 "Error, the observation function dimension no="<<no<<
" > 1 is not allowed" 940 " when the observation is passed through as the observation matching term!" 945 RCP<const VectorBase<Scalar> > o;
946 RCP<VectorBase<Scalar> > diff_o;
947 if( !observationPassThrough_ && ( nonnull(g_inv_out) || nonnull(DgDx_inv_trans_out) ) )
949 if (obs_idx_ < 0 ) o = x;
else o = wrapped_o;
950 if(trace) *out <<
"\n||o||inf = " << norm_inf(*o) << endl;
951 if (print_o) *out <<
"\no = " << *o;
952 diff_o = createMember(o_space);
953 RCP<const VectorBase<Scalar> >
955 = ( observationTargetAsParameter_
956 ? inArgs.get_p(inArgs.Np()-1)
959 if (is_null(observationTarget) ) {
960 observationTarget = observationTarget_;
962 *out <<
"\n||ot||inf = " << norm_inf(*observationTarget) << endl;
964 *out <<
"\not = " << *observationTarget;
966 if (!is_null(observationTarget)) {
967 V_VmV( diff_o.ptr(), *o, *observationTarget );
970 assign( diff_o.ptr(), *o );
973 *out <<
"\n||diff_o||inf = " << norm_inf(*diff_o) << endl;
976 *out <<
"\ndiff_o = " << *diff_o;
981 RCP<VectorBase<Scalar> > diff_p;
982 if ( nonnull(g_inv_out) || nonnull(DgDp_inv_trans_out)) {
983 if(trace) *out <<
"\n||p||inf = " << norm_inf(*p) << endl;
984 if(print_p) *out <<
"\np = " << Teuchos::describe(*p,Teuchos::VERB_EXTREME);
985 diff_p = createMember(p_space);
986 if (!is_null(parameterBase_) ) {
987 if(trace) *out <<
"\n||pt||inf = " << norm_inf(*parameterBase_) << endl;
990 << Teuchos::describe(*parameterBase_,Teuchos::VERB_EXTREME);
992 V_VmV( diff_p.ptr(), *p, *parameterBase_ );
995 assign( diff_p.ptr(), *p );
997 if(trace) {*out <<
"\n||diff_p|| = " << norm(*diff_p) << endl;}
999 *out <<
"\ndiff_p = " 1000 << Teuchos::describe(*diff_p, Teuchos::VERB_EXTREME);
1007 RCP<const LinearOpBase<Scalar> >
1008 Q_o = this->get_observationMatchWeightingOp(),
1009 Q_p = this->get_parameterRegularizationWeightingOp();
1011 #ifdef TEUCHOS_DEBUG 1012 if (!is_null(Q_o)) {
1014 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1015 *Q_o->range(), *o_space
1018 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1019 *Q_o->domain(), *o_space
1022 if (!is_null(Q_p)) {
1024 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1025 *Q_p->range(), *p_space
1028 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1029 *Q_p->domain(), *p_space
1037 RCP<VectorBase<Scalar> > Q_o_diff_o;
1038 if ( !is_null(Q_o) && !is_null(diff_o) ) {
1039 Q_o_diff_o = createMember(Q_o->range());
1040 apply( *Q_o,
NOTRANS, *diff_o, Q_o_diff_o.ptr() );
1044 RCP<VectorBase<Scalar> > Q_p_diff_p;
1045 if ( !is_null(Q_p) && !is_null(diff_p) ) {
1046 Q_p_diff_p = createMember(Q_p->range());
1047 apply( *Q_p,
NOTRANS, *diff_p, Q_p_diff_p.ptr() );
1051 if (nonnull(g_inv_out)) {
1053 *out <<
"\nComputing inverse response function ginv = g(Np-1) ...\n";
1054 const Scalar observationTerm
1055 = ( observationPassThrough_
1056 ? get_ele(*wrapped_o,0)
1057 : ( observationMultiplier_ != ST::zero()
1059 ? observationMultiplier_*0.5*dot(*diff_o,*Q_o_diff_o)
1060 : observationMultiplier_*(0.5/no)*dot(*diff_o,*diff_o)
1065 const Scalar parameterTerm
1066 = ( parameterMultiplier_ != ST::zero()
1068 ? parameterMultiplier_*0.5*dot(*diff_p,*Q_p_diff_p)
1069 : parameterMultiplier_*(0.5/np)*dot(*diff_p,*diff_p)
1073 const Scalar g_inv_val = observationTerm+parameterTerm;
1076 <<
"\nObservation matching term of ginv = g(Np-1):" 1077 <<
"\n observationMultiplier = " << observationMultiplier_
1078 <<
"\n observationMultiplier*observationMatch(x,p) = " << observationTerm
1079 <<
"\nParameter regularization term of ginv = g(Np-1):" 1080 <<
"\n parameterMultiplier = " << parameterMultiplier_
1081 <<
"\n parameterMultiplier*parameterRegularization(p) = " << parameterTerm
1082 <<
"\nginv = " << g_inv_val
1084 set_ele(0, observationTerm+parameterTerm, g_inv_out.ptr());
1088 if (nonnull(DgDx_inv_trans_out)) {
1090 *out <<
"\nComputing inverse response function derivative DginvDx^T:\n";
1091 if (!observationPassThrough_) {
1092 if( obs_idx_ < 0 ) {
1093 if (!is_null(Q_o)) {
1095 *out <<
"\nDginvDx^T = observationMultiplier * Q_o * diff_o ...\n";
1097 DgDx_inv_trans_out->col(0).ptr(),
1098 observationMultiplier_,
1104 *out <<
"\nDginvDx^T = observationMultiplier * (1/no) * diff_o ...\n";
1106 DgDx_inv_trans_out->col(0).ptr(),
1107 Scalar(observationMultiplier_*(1.0/no)),
1115 if (print_o && print_x)
1116 *out <<
"\nDoDx = " << *wrapped_DoDx.getLinearOp();
1117 if (!is_null(Q_o)) {
1119 *out <<
"\nDginvDx^T = observationMultiplier * DoDx^T * Q_o * diff_o ...\n";
1123 DgDx_inv_trans_out->col(0).ptr(),
1124 observationMultiplier_
1129 *out <<
"\nDginvDx^T = (observationMultiplier*(1/no)) * DoDx^T * diff_o ...\n";
1133 DgDx_inv_trans_out->col(0).ptr(),
1134 Scalar(observationMultiplier_*(1.0/no))
1141 *out <<
"\nDginvDx^T = observationMultiplier * DoDx^T ...\n";
1143 DgDx_inv_trans_out->col(0).ptr(), observationMultiplier_,
1144 *wrapped_DoDx.getMultiVector()->col(0)
1148 *out <<
"\n||DginvDx^T||inf = " << norms_inf(*DgDx_inv_trans_out) <<
"\n";
1150 *out <<
"\nDginvDx^T = " << *DgDx_inv_trans_out;
1154 if (nonnull(DgDp_inv_trans_out)) {
1156 *out <<
"\nComputing inverse response function derivative DginvDp^T ...\n";
1157 if (obs_idx_ >= 0) {
1159 *out <<
"\n||DoDp^T|| = " << norms_inf(*wrapped_DoDp_trans.getMultiVector()) << endl;
1161 *out <<
"\nDoDp^T = " << Teuchos::describe(*wrapped_DoDp_trans.getMultiVector(),Teuchos::VERB_EXTREME);
1164 *out <<
"\nDginvDp^T = 0 ...\n";
1165 assign( DgDp_inv_trans_out->col(0).ptr(), ST::zero() );
1167 if (!observationPassThrough_) {
1168 if ( obs_idx_ >= 0 ) {
1169 if ( !is_null(Q_o) ) {
1171 *out <<
"\nDginvDp^T += observationMultiplier* * (DoDp^T) * Q_o * diff_o ...\n";
1173 *wrapped_DoDp_trans.getMultiVector(),
NOTRANS,
1175 DgDp_inv_trans_out->col(0).ptr(),
1176 Scalar(observationMultiplier_*(1.0/no)),
1182 *out <<
"\nDgDp^T += observationMultiplier* * (DoDp^T) * Q_o * diff_o ...\n";
1184 *wrapped_DoDp_trans.getMultiVector(),
NOTRANS,
1186 DgDp_inv_trans_out->col(0).ptr(),
1187 Scalar(observationMultiplier_*(1.0/no)),
1192 *out <<
"\n||DginvDp^T||inf = " << norms_inf(*DgDp_inv_trans_out) <<
"\n";
1194 *out <<
"\nDginvDp^T = " << *DgDp_inv_trans_out;
1202 *out <<
"\nDginvDp^T += (observationMultiplier*(1/no)) * (DoDp^T) * diff_o ...\n";
1204 DgDp_inv_trans_out->col(0).ptr(), observationMultiplier_,
1205 *wrapped_DoDp_trans.getMultiVector()->col(0)
1210 if( parameterMultiplier_ != ST::zero() ) {
1211 if ( !is_null(Q_p) ) {
1213 *out <<
"\nDginvDp^T += parameterMultiplier * Q_p * diff_p ...\n";
1215 DgDp_inv_trans_out->col(0).ptr(),
1216 parameterMultiplier_,
1222 *out <<
"\nDginvDp^T += (parameterMultiplier*(1.0/np)) * diff_p ...\n";
1224 DgDp_inv_trans_out->col(0).ptr(),
1225 Scalar(parameterMultiplier_*(1.0/np)),
1230 *out <<
"\n||DginvDp^T||inf = " << norms_inf(*DgDp_inv_trans_out) <<
"\n";
1232 *out <<
"\nDginvDp^T = " << *DgDp_inv_trans_out;
1241 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
1249 template<
class Scalar>
1250 void DefaultInverseModelEvaluator<Scalar>::initializeDefaults()
1252 obs_idx_ = ObservationIndex_default_;
1253 p_idx_ = ParameterSubvectorIndex_default_;
1254 observationMultiplier_ = ObservationMultiplier_default_;
1255 parameterMultiplier_ = ParameterMultiplier_default_;
1259 template<
class Scalar>
1260 void DefaultInverseModelEvaluator<Scalar>::initializeInArgsOutArgs()
const 1263 typedef ModelEvaluatorBase MEB;
1265 const RCP<const ModelEvaluator<Scalar> >
1266 thyraModel = this->getUnderlyingModel();
1268 const MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
1269 const int wrapped_Np = wrappedInArgs.Np();
1271 MEB::InArgsSetup<Scalar> inArgs;
1272 inArgs.setModelEvalDescription(this->description());
1273 const bool supports_x = wrappedInArgs.supports(MEB::IN_ARG_x);
1274 usingObservationTargetAsParameter_ = ( supports_x && observationTargetAsParameter_ );
1277 wrapped_Np + ( usingObservationTargetAsParameter_ ? 1 : 0 )
1279 prototypeInArgs_ = inArgs;
1281 const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
1282 const int wrapped_Ng = wrappedOutArgs.Ng();
1284 MEB::OutArgsSetup<Scalar> outArgs;
1285 outArgs.setModelEvalDescription(inArgs.modelEvalDescription());
1286 outArgs.set_Np_Ng( prototypeInArgs_.Np(), wrapped_Ng+1 );
1287 outArgs.setSupports(wrappedOutArgs);
1288 outArgs.setSupports(MEB::OUT_ARG_DgDx,wrapped_Ng,MEB::DERIV_TRANS_MV_BY_ROW);
1289 outArgs.setSupports(MEB::OUT_ARG_DgDp,wrapped_Ng,p_idx_,MEB::DERIV_TRANS_MV_BY_ROW);
1290 prototypeOutArgs_ = outArgs;
1295 template<
class Scalar>
1296 RCP<const VectorSpaceBase<Scalar> >
1297 DefaultInverseModelEvaluator<Scalar>::get_obs_space()
const 1299 return ( obs_idx_ < 0 ? this->get_x_space() : this->get_g_space(obs_idx_) );
1306 #endif // THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible...
RCP< const Teuchos::ParameterList > getParameterList() const
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...
void uninitialize()
Uninitialize.
RCP< const Teuchos::ParameterList > getValidParameters() const
This is a base class that delegetes almost all function to a wrapped model evaluator object...
Use the non-transposed operator.
void initialize(const RCP< ModelEvaluator< Scalar > > &model)
Initialize given a non-const model evaluator.
RCP< const VectorSpaceBase< Scalar > > get_g_space(int j) const
void initialize(const RCP< ModelEvaluator< Scalar > > &thyraModel)
ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
DefaultInverseModelEvaluator()
RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const
RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Abstract strategy interface for reading and writing (multi)vector objects to and from files...
Abstract interface for finite-dimensional dense vectors.
RCP< DefaultInverseModelEvaluator< Scalar > > defaultInverseModelEvaluator(const RCP< ModelEvaluator< Scalar > > &thyraModel)
Non-member constructor.
Base class for all linear operators.
STANDARD_CONST_COMPOSITION_MEMBERS(VectorBase< Scalar >, observationTarget)
Observation target vector ot.
This class wraps any ModelEvaluator object and adds a simple, but fairly general, inverse response fu...
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
RCP< Teuchos::ParameterList > getNonconstParameterList()
std::string description() const
STANDARD_NONCONST_COMPOSITION_MEMBERS(MultiVectorFileIOBase< Scalar >, observationTargetIO)
MultiVectorFileIOBase object used to read the observation target vector ot as directed by the paramet...
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...