Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultInverseModelEvaluator.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP
43#define THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP
44
45
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"
55#include "Teuchos_ParameterListAcceptor.hpp"
56#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
57#include "Teuchos_StandardParameterEntryValidators.hpp"
58#include "Teuchos_Time.hpp"
59
60
61namespace Thyra {
62
63
262template<class Scalar>
264 : virtual public ModelEvaluatorDelegatorBase<Scalar>
265 , virtual public Teuchos::ParameterListAcceptor
266{
267public:
268
271
274
277
279 STANDARD_CONST_COMPOSITION_MEMBERS( LinearOpBase<Scalar>, parameterRegularizationWeightingOp );
280
284
288
291
294
296 void initialize(
297 const RCP<ModelEvaluator<Scalar> > &thyraModel
298 );
299
301 void uninitialize(
302 RCP<ModelEvaluator<Scalar> > *thyraModel
303 );
304
306
309
317 void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
333
335
338
340 std::string description() const;
341
343
346
353
355
356private:
357
360
362 ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
364 void evalModelImpl(
367 ) const;
368
370
371private:
372
373 // ////////////////////////////////
374 // Private data members
375
376 mutable RCP<const Teuchos::ParameterList> validParamList_;
378
380
381 mutable ModelEvaluatorBase::InArgs<Scalar> prototypeInArgs_;
382 mutable ModelEvaluatorBase::OutArgs<Scalar> prototypeOutArgs_;
383 mutable bool usingObservationTargetAsParameter_;
384
385 int obs_idx_;
386 int p_idx_;
387
388
389 double observationMultiplier_;
390 double parameterMultiplier_;
391
392 bool observationTargetAsParameter_;
393
394 bool observationPassThrough_;
395
396 Teuchos::EVerbosityLevel localVerbLevel_;
397
398 mutable ParameterDrivenMultiVectorInput<Scalar> observationTargetReader_;
399 mutable ParameterDrivenMultiVectorInput<Scalar> parameterBaseReader_;
400
401 static const std::string ObservationIndex_name_;
402 static const int ObservationIndex_default_;
403
404 static const std::string ParameterSubvectorIndex_name_;
405 static const int ParameterSubvectorIndex_default_;
406
407 static const std::string ObservationMultiplier_name_;
408 static const double ObservationMultiplier_default_;
409
410 static const std::string ObservationTargetVector_name_;
411
412 static const std::string ObservationTargetAsParameter_name_;
413 static const bool ObservationTargetAsParameter_default_;
414
415 static const std::string ObservationPassThrough_name_;
416 static const bool ObservationPassThrough_default_;
417
418 static const std::string ParameterMultiplier_name_;
419 static const double ParameterMultiplier_default_;
420
421 static const std::string ParameterBaseVector_name_;
422
423 // ////////////////////////////////
424 // Private member functions
425
426 void initializeDefaults();
427
428 void initializeInArgsOutArgs() const;
429
430 RCP<const VectorSpaceBase<Scalar> > get_obs_space() const;
431
432};
433
434
439template<class Scalar>
442 const RCP<ModelEvaluator<Scalar> > &thyraModel
443 )
444{
447 inverseModel->initialize(thyraModel);
448 return inverseModel;
449}
450
451
452// /////////////////////////////////
453// Implementations
454
455
456// Static data members
457
458
459template<class Scalar>
460const std::string
462= "Observation Index";
463
464template<class Scalar>
465const int
467= -1;
468
469
470template<class Scalar>
471const std::string
473= "Parameter Subvector Ordinal";
474
475template<class Scalar>
476const int
478= 0;
479
480
481template<class Scalar>
482const std::string
484= "Observation Multiplier";
485
486template<class Scalar>
487const double
489= 1.0;
490
491
492template<class Scalar>
493const std::string
495= "Observation Target Vector";
496
497
498template<class Scalar>
499const std::string
501= "Observation Target as Parameter";
502
503template<class Scalar>
504const bool
506= false;
507
508
509template<class Scalar>
510const std::string
512= "Observation Pass Through";
513
514template<class Scalar>
515const bool
517= false;
518
519
520template<class Scalar>
521const std::string
523= "Parameter Multiplier";
524
525template<class Scalar>
526const double
528= 1e-6;
529
530
531template<class Scalar>
532const std::string
534= "Parameter Base Vector";
535
536
537// Constructors/initializers/accessors/utilities
538
539
540template<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)
548{}
549
550
551template<class Scalar>
553 const RCP<ModelEvaluator<Scalar> > &thyraModel
554 )
555{
557 inv_g_space_= thyraModel->get_x_space()->smallVecSpcFcty()->createVecSpc(1);
558 // Get ready for reinitalization
559 prototypeInArgs_ = ModelEvaluatorBase::InArgs<Scalar>();
560 prototypeOutArgs_ = ModelEvaluatorBase::OutArgs<Scalar>();
561}
562
563
564template<class Scalar>
566 RCP<ModelEvaluator<Scalar> > *thyraModel
567 )
568{
569 if(thyraModel) *thyraModel = this->getUnderlyingModel();
571}
572
573
574// Overridden from Teuchos::ParameterListAcceptor
575
576
577template<class Scalar>
579 RCP<Teuchos::ParameterList> const& paramList
580 )
581{
582
583 using Teuchos::Array;
584 using Teuchos::getParameterPtr;
585 using Teuchos::rcp;
586 using Teuchos::sublist;
587
588 // Validate and set the parameter list
589 TEUCHOS_TEST_FOR_EXCEPT(0==paramList.get());
590 paramList->validateParameters(*getValidParameters(),0);
591 paramList_ = paramList;
592
593 // Parameters for observation matching term
594 obs_idx_ = paramList_->get(
595 ObservationIndex_name_,ObservationIndex_default_);
596 observationPassThrough_ = paramList_->get(
597 ObservationPassThrough_name_, ObservationPassThrough_default_ );
598#ifdef TEUCHOS_DEBUG
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!"
603 );
604#endif
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());
613 vots_observationTargetReader(
614 rcp(&observationTargetReader_,false)
615 ,this->getOStream(),this->getVerbLevel()
616 );
617 observationTargetReader_.setParameterList(
618 sublist(paramList_,ObservationTargetVector_name_)
619 );
621 observationTarget;
622 observationTargetReader_.readVector(
623 "observation target vector",&observationTarget);
624 observationTarget_ = observationTarget;
625 }
626 }
627 else {
628 observationTargetAsParameter_ = false;
629 observationTarget_ = Teuchos::null;
630 }
631
632 // Parameters for parameter matching term
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_));
640 vots_parameterBaseReader(
641 rcp(&parameterBaseReader_,false)
642 ,this->getOStream(),this->getVerbLevel()
643 );
644 parameterBaseReader_.setParameterList(
645 sublist(paramList_,ParameterBaseVector_name_)
646 );
648 parameterBase;
649 parameterBaseReader_.readVector(
650 "parameter base vector",&parameterBase);
651 parameterBase_ = parameterBase;
652 }
653
654 // Verbosity settings
655 localVerbLevel_ = this->readLocalVerbosityLevelValidatedParameter(*paramList_);
656 Teuchos::readVerboseObjectSublist(&*paramList_,this);
657
658#ifdef TEUCHOS_DEBUG
659 paramList_->validateParameters(*getValidParameters(),0);
660#endif // TEUCHOS_DEBUG
661
662 // Get ready for reinitalization
663 prototypeInArgs_ = ModelEvaluatorBase::InArgs<Scalar>();
664 prototypeOutArgs_ = ModelEvaluatorBase::OutArgs<Scalar>();
665
666}
667
668
669template<class Scalar>
675
676
677template<class Scalar>
680{
681 RCP<Teuchos::ParameterList> _paramList = paramList_;
682 paramList_ = Teuchos::null;
683 return _paramList;
684}
685
686
687template<class Scalar>
690{
691 return paramList_;
692}
693
694
695template<class Scalar>
698{
699 if(validParamList_.get()==NULL) {
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)."
706 );
707 pl->set( ParameterSubvectorIndex_name_,ParameterSubvectorIndex_default_,
708 "The index of the parameter subvector that will be used in the\n"
709 "regularization term."
710 );
711 pl->set( ObservationMultiplier_name_,ObservationMultiplier_default_,
712 "observationMultiplier"
713 );
714 if(this->get_observationTargetIO().get())
715 observationTargetReader_.set_fileIO(this->get_observationTargetIO());
716 pl->sublist(ObservationTargetVector_name_).setParameters(
717 *observationTargetReader_.getValidParameters()
718 );
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"
722 "response function!"
723 );
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."
727 );
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()
734 );
735 this->setLocalVerbosityLevelValidatedParameter(&*pl);
736 Teuchos::setupVerboseObjectSublist(&*pl);
737 validParamList_ = pl;
738 }
739 return validParamList_;
740}
741
742
743// Overridden from ModelEvaulator.
744
745
746template<class Scalar>
749{
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);
755}
756
757
758template<class Scalar>
761{
762 if (prototypeOutArgs_.Np()==0)
763 initializeInArgsOutArgs();
764 if (j==prototypeOutArgs_.Ng()-1)
765 return inv_g_space_;
766 return this->getUnderlyingModel()->get_g_space(j);
767}
768
769
770template<class Scalar>
773{
774 if (prototypeInArgs_.Np()==0)
775 initializeInArgsOutArgs();
776 return prototypeInArgs_;
777}
778
779
780// Public functions overridden from Teuchos::Describable
781
782
783template<class Scalar>
785{
787 thyraModel = this->getUnderlyingModel();
788 std::ostringstream oss;
789 oss << "Thyra::DefaultInverseModelEvaluator{";
790 oss << "thyraModel=";
791 if(thyraModel.get())
792 oss << "\'"<<thyraModel->description()<<"\'";
793 else
794 oss << "NULL";
795 oss << "}";
796 return oss.str();
797}
798
799
800// Private functions overridden from ModelEvaulatorDefaultBase
801
802
803template<class Scalar>
806{
807 if (prototypeOutArgs_.Np()==0)
808 initializeInArgsOutArgs();
809 return prototypeOutArgs_;
810}
811
812
813template<class Scalar>
814void DefaultInverseModelEvaluator<Scalar>::evalModelImpl(
815 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
816 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
817 ) const
818{
819
820 using std::endl;
821 using Teuchos::rcp;
822 using Teuchos::rcp_const_cast;
823 using Teuchos::rcp_dynamic_cast;
824 using Teuchos::OSTab;
826 typedef typename ST::magnitudeType ScalarMag;
827 typedef ModelEvaluatorBase MEB;
828
829 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_LOCALVERBLEVEL_BEGIN(
830 "Thyra::DefaultInverseModelEvaluator",inArgs,outArgs,localVerbLevel_
831 );
832
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;
837
838 //
839 // A) See what needs to be computed
840 //
841
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
847 );
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
851 );
852 const bool computeInverseFunction = ( nonnull(g_inv_out)
853 || nonnull(DgDx_inv_trans_out) || nonnull(DgDp_inv_trans_out) );
854
855 //
856 // B) Compute all of the needed functions from the base model
857 //
858
859 if(trace)
860 *out << "\nComputing the base point and the observation(s) ...\n";
861
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 )
870 {
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_);
876 else
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);
880 }
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
884 );
885 wrappedOutArgs.set_DgDp(obs_idx_,p_idx_,wrapped_DoDp_trans);
886 }
887 // 2007/07/28: rabartl: Above, we really should check if these output
888 // arguments have already been set by the client. If they are, then we
889 // need to make sure that they are of the correct form or we need to throw
890 // an exception!
891 }
892
893 if (!wrappedOutArgs.isEmpty()) {
894 thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
895 }
896 else {
897 if(trace)
898 *out << "\nSkipping the evaluation of the underlying model since "
899 << "there is nothing to compute ...\n";
900 }
901
902 bool failed = wrappedOutArgs.isFailed();
903
904 //
905 // C) Assemble the final observation and paramter terms
906 //
907
908 if ( !failed && computeInverseFunction ) {
909
910 //
911 // Compute the inverse response function and its derivatives
912 //
913
914 RCP<const VectorBase<Scalar> >
915 x_in = inArgs.get_x(),
916 p_in = inArgs.get_p(p_idx_);
917
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() );
922
923 const RCP<const VectorSpaceBase<Scalar> >
924 o_space = get_obs_space(),
925 p_space = this->get_p_space(p_idx_);
926
927 const Ordinal
928 no = o_space->dim(),
929 np = p_space->dim();
930
931 if (trace)
932 *out << "\nno = " << no
933 << "\nnp = " << np
934 << endl;
935
936#ifdef TEUCHOS_DEBUG
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!"
941 );
942#endif
943
944 // Compute diff_o if needed
945 RCP<const VectorBase<Scalar> > o;
946 RCP<VectorBase<Scalar> > diff_o;
947 if( !observationPassThrough_ && ( nonnull(g_inv_out) || nonnull(DgDx_inv_trans_out) ) )
948 {
949 if (obs_idx_ < 0 ) o = x; else o = wrapped_o; // can't use ( test ? x : 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> >
954 observationTarget
955 = ( observationTargetAsParameter_
956 ? inArgs.get_p(inArgs.Np()-1)
958 );
959 if (is_null(observationTarget) ) {
960 observationTarget = observationTarget_;
961 if (trace)
962 *out << "\n||ot||inf = " << norm_inf(*observationTarget) << endl;
963 if (print_o)
964 *out << "\not = " << *observationTarget;
965 }
966 if (!is_null(observationTarget)) {
967 V_VmV( diff_o.ptr(), *o, *observationTarget );
968 }
969 else {
970 assign( diff_o.ptr(), *o );
971 }
972 if(trace) {
973 *out << "\n||diff_o||inf = " << norm_inf(*diff_o) << endl;
974 }
975 if (print_o) {
976 *out << "\ndiff_o = " << *diff_o;
977 }
978 }
979
980 // Compute diff_p if needed
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;
988 if(print_p) {
989 *out << "\npt = "
990 << Teuchos::describe(*parameterBase_,Teuchos::VERB_EXTREME);
991 }
992 V_VmV( diff_p.ptr(), *p, *parameterBase_ );
993 }
994 else {
995 assign( diff_p.ptr(), *p );
996 }
997 if(trace) {*out << "\n||diff_p|| = " << norm(*diff_p) << endl;}
998 if(print_p) {
999 *out << "\ndiff_p = "
1000 << Teuchos::describe(*diff_p, Teuchos::VERB_EXTREME);
1001 }
1002 }
1003
1004
1005 // Get and check Q_o and Q_p
1006
1007 RCP<const LinearOpBase<Scalar> >
1008 Q_o = this->get_observationMatchWeightingOp(),
1009 Q_p = this->get_parameterRegularizationWeightingOp();
1010
1011#ifdef TEUCHOS_DEBUG
1012 if (!is_null(Q_o)) {
1014 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1015 *Q_o->range(), *o_space
1016 );
1018 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1019 *Q_o->domain(), *o_space
1020 );
1021 }
1022 if (!is_null(Q_p)) {
1024 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1025 *Q_p->range(), *p_space
1026 );
1028 "Thyra::DefaultInverseModelEvaluator::evalModel(...)",
1029 *Q_p->domain(), *p_space
1030 );
1031 }
1032 // Note, we have not proved that Q_o and Q_p are s.p.d. but at least we
1033 // have established that that have the right range and domain spaces!
1034#endif
1035
1036 // Compute Q_o * diff_o
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()); // Should be same as domain!
1040 apply( *Q_o, NOTRANS, *diff_o, Q_o_diff_o.ptr() );
1041 }
1042
1043 // Compute Q_p * diff_p
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()); // Should be same as domain!
1047 apply( *Q_p, NOTRANS, *diff_p, Q_p_diff_p.ptr() );
1048 }
1049
1050 // Compute g_inv(x,p)
1051 if (nonnull(g_inv_out)) {
1052 if(trace)
1053 *out << "\nComputing inverse response function ginv = g(Np-1) ...\n";
1054 const Scalar observationTerm
1055 = ( observationPassThrough_
1056 ? get_ele(*wrapped_o,0) // ToDo; Verify that this is already a scalar
1057 : ( observationMultiplier_ != ST::zero()
1058 ? ( !is_null(Q_o)
1059 ? observationMultiplier_*0.5*dot(*diff_o,*Q_o_diff_o)
1060 : observationMultiplier_*(0.5/no)*dot(*diff_o,*diff_o)
1061 )
1062 : ST::zero()
1063 )
1064 );
1065 const Scalar parameterTerm
1066 = ( parameterMultiplier_ != ST::zero()
1067 ? ( !is_null(Q_p)
1068 ? parameterMultiplier_*0.5*dot(*diff_p,*Q_p_diff_p)
1069 : parameterMultiplier_*(0.5/np)*dot(*diff_p,*diff_p)
1070 )
1071 : ST::zero()
1072 );
1073 const Scalar g_inv_val = observationTerm+parameterTerm;
1074 if(trace)
1075 *out
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
1083 << "\n";
1084 set_ele(0, observationTerm+parameterTerm, g_inv_out.ptr());
1085 }
1086
1087 // Compute d(g_inv)/d(x)^T
1088 if (nonnull(DgDx_inv_trans_out)) {
1089 if(trace)
1090 *out << "\nComputing inverse response function derivative DginvDx^T:\n";
1091 if (!observationPassThrough_) {
1092 if( obs_idx_ < 0 ) {
1093 if (!is_null(Q_o)) {
1094 if (trace)
1095 *out << "\nDginvDx^T = observationMultiplier * Q_o * diff_o ...\n";
1096 V_StV(
1097 DgDx_inv_trans_out->col(0).ptr(),
1098 observationMultiplier_,
1099 *Q_o_diff_o
1100 );
1101 }
1102 else {
1103 if (trace)
1104 *out << "\nDginvDx^T = observationMultiplier * (1/no) * diff_o ...\n";
1105 V_StV(
1106 DgDx_inv_trans_out->col(0).ptr(),
1107 Scalar(observationMultiplier_*(1.0/no)),
1108 *diff_o
1109 );
1110 }
1111 }
1112 else {
1113 //if (trace)
1114 // *out << "\n||DoDx^T||inf = " << norms_inf(*wrapped_DoDx.getMultiVector()) << endl;
1115 if (print_o && print_x)
1116 *out << "\nDoDx = " << *wrapped_DoDx.getLinearOp();
1117 if (!is_null(Q_o)) {
1118 if (trace)
1119 *out << "\nDginvDx^T = observationMultiplier * DoDx^T * Q_o * diff_o ...\n";
1120 apply(
1121 *wrapped_DoDx.getLinearOp(), CONJTRANS,
1122 *Q_o_diff_o,
1123 DgDx_inv_trans_out->col(0).ptr(),
1124 observationMultiplier_
1125 );
1126 }
1127 else {
1128 if (trace)
1129 *out << "\nDginvDx^T = (observationMultiplier*(1/no)) * DoDx^T * diff_o ...\n";
1130 apply(
1131 *wrapped_DoDx.getLinearOp(), CONJTRANS,
1132 *diff_o,
1133 DgDx_inv_trans_out->col(0).ptr(),
1134 Scalar(observationMultiplier_*(1.0/no))
1135 );
1136 }
1137 }
1138 }
1139 else {
1140 if (trace)
1141 *out << "\nDginvDx^T = observationMultiplier * DoDx^T ...\n";
1142 V_StV(
1143 DgDx_inv_trans_out->col(0).ptr(), observationMultiplier_,
1144 *wrapped_DoDx.getMultiVector()->col(0)
1145 );
1146 }
1147 if(trace)
1148 *out << "\n||DginvDx^T||inf = " << norms_inf(*DgDx_inv_trans_out) << "\n";
1149 if (print_x)
1150 *out << "\nDginvDx^T = " << *DgDx_inv_trans_out;
1151 }
1152
1153 // Compute d(g_inv)/d(p)^T
1154 if (nonnull(DgDp_inv_trans_out)) {
1155 if(trace)
1156 *out << "\nComputing inverse response function derivative DginvDp^T ...\n";
1157 if (obs_idx_ >= 0) {
1158 if (trace)
1159 *out << "\n||DoDp^T|| = " << norms_inf(*wrapped_DoDp_trans.getMultiVector()) << endl;
1160 if (print_p)
1161 *out << "\nDoDp^T = " << Teuchos::describe(*wrapped_DoDp_trans.getMultiVector(),Teuchos::VERB_EXTREME);
1162 }
1163 if(trace)
1164 *out << "\nDginvDp^T = 0 ...\n";
1165 assign( DgDp_inv_trans_out->col(0).ptr(), ST::zero() );
1166 // DgDp^T += observationMultiplier * d(observationMatch)/d(p)^T
1167 if (!observationPassThrough_) {
1168 if ( obs_idx_ >= 0 ) {
1169 if ( !is_null(Q_o) ) {
1170 if(trace)
1171 *out << "\nDginvDp^T += observationMultiplier* * (DoDp^T) * Q_o * diff_o ...\n";
1172 apply(
1173 *wrapped_DoDp_trans.getMultiVector(), NOTRANS,
1174 *Q_o_diff_o,
1175 DgDp_inv_trans_out->col(0).ptr(),
1176 Scalar(observationMultiplier_*(1.0/no)),
1177 ST::one()
1178 );
1179 }
1180 else {
1181 if(trace)
1182 *out << "\nDgDp^T += observationMultiplier* * (DoDp^T) * Q_o * diff_o ...\n";
1183 apply(
1184 *wrapped_DoDp_trans.getMultiVector(), NOTRANS,
1185 *diff_o,
1186 DgDp_inv_trans_out->col(0).ptr(),
1187 Scalar(observationMultiplier_*(1.0/no)),
1188 ST::one()
1189 );
1190 }
1191 if(trace)
1192 *out << "\n||DginvDp^T||inf = " << norms_inf(*DgDp_inv_trans_out) << "\n";
1193 if (print_p)
1194 *out << "\nDginvDp^T = " << *DgDp_inv_trans_out;
1195 }
1196 else {
1197 // d(observationMatch)/d(p)^T = 0, nothing to do!
1198 }
1199 }
1200 else {
1201 if(trace)
1202 *out << "\nDginvDp^T += (observationMultiplier*(1/no)) * (DoDp^T) * diff_o ...\n";
1203 Vp_StV(
1204 DgDp_inv_trans_out->col(0).ptr(), observationMultiplier_,
1205 *wrapped_DoDp_trans.getMultiVector()->col(0)
1206 );
1207
1208 }
1209 // DgDp^T += parameterMultiplier * d(parameterRegularization)/d(p)^T
1210 if( parameterMultiplier_ != ST::zero() ) {
1211 if ( !is_null(Q_p) ) {
1212 if(trace)
1213 *out << "\nDginvDp^T += parameterMultiplier * Q_p * diff_p ...\n";
1214 Vp_StV(
1215 DgDp_inv_trans_out->col(0).ptr(),
1216 parameterMultiplier_,
1217 *Q_p_diff_p
1218 );
1219 }
1220 else {
1221 if(trace)
1222 *out << "\nDginvDp^T += (parameterMultiplier*(1.0/np)) * diff_p ...\n";
1223 Vp_StV(
1224 DgDp_inv_trans_out->col(0).ptr(),
1225 Scalar(parameterMultiplier_*(1.0/np)),
1226 *diff_p
1227 );
1228 }
1229 if(trace)
1230 *out << "\n||DginvDp^T||inf = " << norms_inf(*DgDp_inv_trans_out) << "\n";
1231 if (print_p)
1232 *out << "\nDginvDp^T = " << *DgDp_inv_trans_out;
1233 }
1234 else {
1235 // This term is zero so there is nothing to do!
1236 }
1237 }
1238
1239 }
1240
1241 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
1242
1243}
1244
1245
1246// private
1247
1248
1249template<class Scalar>
1250void DefaultInverseModelEvaluator<Scalar>::initializeDefaults()
1251{
1252 obs_idx_ = ObservationIndex_default_;
1253 p_idx_ = ParameterSubvectorIndex_default_;
1254 observationMultiplier_ = ObservationMultiplier_default_;
1255 parameterMultiplier_ = ParameterMultiplier_default_;
1256}
1257
1258
1259template<class Scalar>
1260void DefaultInverseModelEvaluator<Scalar>::initializeInArgsOutArgs() const
1261{
1262
1263 typedef ModelEvaluatorBase MEB;
1264
1265 const RCP<const ModelEvaluator<Scalar> >
1266 thyraModel = this->getUnderlyingModel();
1267
1268 const MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
1269 const int wrapped_Np = wrappedInArgs.Np();
1270
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_ );
1275 inArgs.setSupports(
1276 wrappedInArgs,
1277 wrapped_Np + ( usingObservationTargetAsParameter_ ? 1 : 0 )
1278 );
1279 prototypeInArgs_ = inArgs;
1280
1281 const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
1282 const int wrapped_Ng = wrappedOutArgs.Ng();
1283
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;
1291
1292}
1293
1294
1295template<class Scalar>
1296RCP<const VectorSpaceBase<Scalar> >
1297DefaultInverseModelEvaluator<Scalar>::get_obs_space() const
1298{
1299 return ( obs_idx_ < 0 ? this->get_x_space() : this->get_g_space(obs_idx_) );
1300}
1301
1302
1303} // namespace Thyra
1304
1305
1306#endif // THYRA_DEFAUL_INVERSE_MODEL_EVALUATOR_HPP
T * get() const
This class wraps any ModelEvaluator object and adds a simple, but fairly general, inverse response fu...
RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const
STANDARD_CONST_COMPOSITION_MEMBERS(VectorBase< Scalar >, observationTarget)
Observation target vector ot.
STANDARD_CONST_COMPOSITION_MEMBERS(LinearOpBase< Scalar >, observationMatchWeightingOp)
Observation match weighting operator Q_o.
STANDARD_NONCONST_COMPOSITION_MEMBERS(MultiVectorFileIOBase< Scalar >, observationTargetIO)
MultiVectorFileIOBase object used to read the observation target vector ot as directed by the paramet...
void initialize(const RCP< ModelEvaluator< Scalar > > &thyraModel)
RCP< const VectorSpaceBase< Scalar > > get_g_space(int j) const
STANDARD_CONST_COMPOSITION_MEMBERS(LinearOpBase< Scalar >, parameterRegularizationWeightingOp)
Parameter regulization weighting operator Q_p.
RCP< const Teuchos::ParameterList > getParameterList() const
RCP< const Teuchos::ParameterList > getValidParameters() const
STANDARD_NONCONST_COMPOSITION_MEMBERS(MultiVectorFileIOBase< Scalar >, parameterBaseIO)
MultiVectorFileIOBase object used to read the parameter base vector pt as directed by the parameter l...
RCP< DefaultInverseModelEvaluator< Scalar > > defaultInverseModelEvaluator(const RCP< ModelEvaluator< Scalar > > &thyraModel)
Non-member constructor.
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
STANDARD_CONST_COMPOSITION_MEMBERS(VectorBase< Scalar >, parameterBase)
Parameter base vector pt.
Base class for all linear operators.
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object.
This is a base class that delegetes almost all function to a wrapped model evaluator object.
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 strategy interface for reading and writing (multi)vector objects to and from files.
Concrete utility class that an ANA can use for reading in a (multi)vector as directed by a parameter ...
Abstract interface for finite-dimensional dense vectors.
#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)
bool nonnull(const std::shared_ptr< T > &p)
#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.
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).
T_To & dyn_cast(T_From &from)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)