Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultLumpedParameterModelEvaluator.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_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
43#define THYRA_DEFAUL_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
44
45
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"
54
55#include "sillyModifiedGramSchmidt.hpp" // This is just an example!
56
57
58namespace Thyra {
59
60
239template<class Scalar>
241 : virtual public ModelEvaluatorDelegatorBase<Scalar>
242 , virtual public Teuchos::ParameterListAcceptor
243{
244public:
245
248
251
253 void initialize(
254 const RCP<ModelEvaluator<Scalar> > &thyraModel
255 );
256
258 void uninitialize(
259 RCP<ModelEvaluator<Scalar> > *thyraModel
260 );
261
262 // 2007/07/30: rabartl: ToDo: Add functions to set and get the underlying
263 // basis matrix!
264
266
269
271 std::string description() const;
272
274
277
279 void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
288
290
293
305 void reportFinalPoint(
306 const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
307 const bool wasSolved
308 );
309
311
312private:
313
316
318 ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
320 void evalModelImpl(
323 ) const;
324
326
327private:
328
329 // ////////////////////////////////
330 // Private data members
331
332 mutable bool isInitialized_;
333 mutable bool nominalValuesAndBoundsUpdated_;
334
335 mutable RCP<const Teuchos::ParameterList> validParamList_;
337
338 // Parameters read from the parameter list
339 int p_idx_;
340 bool autogenerateBasisMatrix_;
341 int numberOfBasisColumns_;
342 bool nominalValueIsParameterBase_;
343 bool ignoreParameterBounds_;
344 Teuchos::EVerbosityLevel localVerbLevel_;
345 bool dumpBasisMatrix_;
346
347 // Reduced affine parameter model
349 mutable RCP<const VectorBase<Scalar> > p_orig_base_;
350
351 // Nominal values and bounds
352 mutable ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
353 mutable ModelEvaluatorBase::InArgs<Scalar> lowerBounds_;
354 mutable ModelEvaluatorBase::InArgs<Scalar> upperBounds_;
355
356 // Static
357
358 static const std::string ParameterSubvectorIndex_name_;
359 static const int ParameterSubvectorIndex_default_;
360
361 static const std::string AutogenerateBasisMatrix_name_;
362 static const bool AutogenerateBasisMatrix_default_;
363
364 static const std::string NumberOfBasisColumns_name_;
365 static const int NumberOfBasisColumns_default_;
366
367 static const std::string NominalValueIsParameterBase_name_;
368 static const bool NominalValueIsParameterBase_default_;
369
370 static const std::string ParameterBaseVector_name_;
371
372 static const std::string IgnoreParameterBounds_name_;
373 static const bool IgnoreParameterBounds_default_;
374
375 static const std::string DumpBasisMatrix_name_;
376 static const bool DumpBasisMatrix_default_;
377
378 // ////////////////////////////////
379 // Private member functions
380
381 // These functions are used to implement late initialization so that the
382 // need for clients to order function calls is reduced.
383
384 // Finish enough initialization to defined spaces etc.
385 void finishInitialization() const;
386
387 // Generate the parameter basis matrix B.
388 void generateParameterBasisMatrix() const;
389
390 // Finish all of initialization needed to define nominal values, bounds, and
391 // p_orig_base. This calls finishInitialization().
392 void updateNominalValuesAndBounds() const;
393
394 // Map from p -> p_orig.
396 map_from_p_to_p_orig( const VectorBase<Scalar> &p ) const;
397
398 // Set up the arguments for DhDp_orig to be computed by the underlying model.
399 void setupWrappedParamDerivOutArgs(
400 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs, // in
401 ModelEvaluatorBase::OutArgs<Scalar> *wrappedOutArgs // in/out
402 ) const;
403
404 // Create DhDp_orig needed to assembled DhDp
406 create_deriv_wrt_p_orig(
409 ) const;
410
411 // After DhDp_orig has been computed, assemble DhDp or DhDp^T for all deriv
412 // output arguments.
413 void assembleParamDerivOutArgs(
414 const ModelEvaluatorBase::OutArgs<Scalar> &wrappedOutArgs, // in
415 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs // in/out
416 ) const;
417
418 // Given a single DhDp_orig, assemble DhDp
419 void assembleParamDeriv(
420 const ModelEvaluatorBase::Derivative<Scalar> &DhDp_orig, // in
421 const ModelEvaluatorBase::Derivative<Scalar> &DhDp // in/out
422 ) const;
423
424};
425
426
431template<class Scalar>
434 const RCP<ModelEvaluator<Scalar> > &thyraModel
435 )
436{
439 paramLumpedModel->initialize(thyraModel);
440 return paramLumpedModel;
441}
442
443
444// /////////////////////////////////
445// Implementations
446
447
448// Static data members
449
450
451template<class Scalar>
452const std::string
454= "Parameter Subvector Index";
455
456template<class Scalar>
457const int
459= 0;
460
461
462template<class Scalar>
463const std::string
465= "Auto-generate Basis Matrix";
466
467template<class Scalar>
468const bool
470= true;
471
472
473template<class Scalar>
474const std::string
476= "Number of Basis Columns";
477
478template<class Scalar>
479const int
481= 1;
482
483
484template<class Scalar>
485const std::string
487= "Nominal Value is Parameter Base";
488
489template<class Scalar>
490const bool
492= true;
493
494
495template<class Scalar>
496const std::string
498= "Parameter Base Vector";
499
500
501template<class Scalar>
502const std::string
504= "Ignore Parameter Bounds";
505
506template<class Scalar>
507const bool
509= false;
510
511
512template<class Scalar>
513const std::string
515= "Dump Basis Matrix";
516
517template<class Scalar>
518const bool
520= false;
521
522
523// Constructors/initializers/accessors/utilities
524
525
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_),
535 localVerbLevel_(Teuchos::VERB_DEFAULT),
536 dumpBasisMatrix_(DumpBasisMatrix_default_)
537{}
538
539
540template<class Scalar>
542 const RCP<ModelEvaluator<Scalar> > &thyraModel
543 )
544{
545 isInitialized_ = false;
546 nominalValuesAndBoundsUpdated_ = false;
548}
549
550
551template<class Scalar>
553 RCP<ModelEvaluator<Scalar> > *thyraModel
554 )
555{
556 isInitialized_ = false;
557 if(thyraModel) *thyraModel = this->getUnderlyingModel();
559}
560
561
562// Public functions overridden from Teuchos::Describable
563
564
565template<class Scalar>
566std::string
568{
570 thyraModel = this->getUnderlyingModel();
571 std::ostringstream oss;
572 oss << "Thyra::DefaultLumpedParameterModelEvaluator{";
573 oss << "thyraModel=";
574 if(thyraModel.get())
575 oss << "\'"<<thyraModel->description()<<"\'";
576 else
577 oss << "NULL";
578 oss << "}";
579 return oss.str();
580}
581
582
583// Overridden from Teuchos::ParameterListAcceptor
584
585
586template<class Scalar>
588 RCP<Teuchos::ParameterList> const& paramList
589 )
590{
591
592 using Teuchos::getParameterPtr;
593 using Teuchos::rcp;
594 using Teuchos::sublist;
595
596 isInitialized_ = false;
597 nominalValuesAndBoundsUpdated_ = false;
598
599 // Validate and set the parameter list
600 TEUCHOS_TEST_FOR_EXCEPT(is_null(paramList));
601 paramList->validateParameters(*getValidParameters(),0);
602 paramList_ = paramList;
603
604 // Read in parameters
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_ );
612 }
613 nominalValueIsParameterBase_ = paramList_->get(
614 NominalValueIsParameterBase_name_, NominalValueIsParameterBase_default_ );
615 if (!nominalValueIsParameterBase_) {
616 TEUCHOS_TEST_FOR_EXCEPT("ToDo: Implement reading parameter base vector from file!");
617 }
618 ignoreParameterBounds_ = paramList_->get(
619 IgnoreParameterBounds_name_, IgnoreParameterBounds_default_ );
620 dumpBasisMatrix_ = paramList_->get(
621 DumpBasisMatrix_name_, DumpBasisMatrix_default_ );
622
623 // Verbosity settings
624 localVerbLevel_ = this->readLocalVerbosityLevelValidatedParameter(*paramList_);
625 Teuchos::readVerboseObjectSublist(&*paramList_,this);
626
627#ifdef TEUCHOS_DEBUG
628 paramList_->validateParameters(*getValidParameters(),0);
629#endif
630
631}
632
633
634template<class Scalar>
640
641
642template<class Scalar>
645{
646 RCP<Teuchos::ParameterList> _paramList = paramList_;
647 paramList_ = Teuchos::null;
648 return _paramList;
649}
650
651
652template<class Scalar>
658
659
660template<class Scalar>
663{
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." );
682 /*
683 if(this->get_parameterBaseIO().get())
684 parameterBaseReader_.set_fileIO(this->get_parameterBaseIO());
685 pl->sublist(ParameterBaseVector_name_).setParameters(
686 *parameterBaseReader_.getValidParameters()
687 );
688 */
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;
698 }
699 return validParamList_;
700}
701
702
703// Overridden from ModelEvaulator.
704
705
706template<class Scalar>
709{
710 finishInitialization();
711 if (l == p_idx_)
712 return B_->domain();
713 return this->getUnderlyingModel()->get_p_space(l);
714}
715
716
717template<class Scalar>
720{
721 finishInitialization();
722 if (l == p_idx_)
723 return Teuchos::null; // Names for these parameters would be meaningless!
724 return this->getUnderlyingModel()->get_p_names(l);
725}
726
727
728template<class Scalar>
731{
732 updateNominalValuesAndBounds();
733 return nominalValues_;
734}
735
736
737template<class Scalar>
740{
741 updateNominalValuesAndBounds();
742 return lowerBounds_;
743}
744
745
746template<class Scalar>
749{
750 updateNominalValuesAndBounds();
751 return upperBounds_;
752}
753
754
755template<class Scalar>
757 const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
758 const bool wasSolved
759 )
760{
761
762 typedef ModelEvaluatorBase MEB;
763
764 // Make sure that everything has been initialized
765 updateNominalValuesAndBounds();
766
768 thyraModel = this->getNonconstUnderlyingModel();
769
770 // By default, copy all input arguments since they will all be the same
771 // except for the given reduced p. We will then replace the reduced p with
772 // p_orig below.
773 MEB::InArgs<Scalar> wrappedFinalPoint = thyraModel->createInArgs();
774 wrappedFinalPoint.setArgs(finalPoint);
775
776 // Replace p with p_orig.
778 if (!is_null(p=finalPoint.get_p(p_idx_))) {
779 wrappedFinalPoint.set_p(p_idx_, map_from_p_to_p_orig(*p));
780 }
781
782 thyraModel->reportFinalPoint(wrappedFinalPoint,wasSolved);
783
784}
785
786
787// Private functions overridden from ModelEvaulatorDefaultBase
788
789
790template<class Scalar>
793{
795 outArgs = this->getUnderlyingModel()->createOutArgs();
796 outArgs.setModelEvalDescription(this->description());
797 return outArgs;
798 // 2007/07/31: rabartl: ToDo: We need to manually set the forms of the
799 // derivatives that this class object will support! This needs to be based
800 // on tests of what the forms of derivatives the underlying model supports.
801}
802
803
804template<class Scalar>
805void DefaultLumpedParameterModelEvaluator<Scalar>::evalModelImpl(
806 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
807 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
808 ) const
809{
810
811 // This routine is pretty simple for the most part. By default, we just
812 // pass everything through to the underlying model evaluator except for
813 // arguments reated to the parameter subvector with index
814 // p_idx_.
815
816 using Teuchos::rcp;
817 using Teuchos::rcp_const_cast;
818 using Teuchos::rcp_dynamic_cast;
819 using Teuchos::OSTab;
821 typedef typename ST::magnitudeType ScalarMag;
822 typedef ModelEvaluatorBase MEB;
823
824 // Make sure that everything has been initialized
825 updateNominalValuesAndBounds();
826
827 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_LOCALVERBLEVEL_BEGIN(
828 "Thyra::DefaultLumpedParameterModelEvaluator",inArgs,outArgs,localVerbLevel_
829 );
830
831 //
832 // A) Setup InArgs
833 //
834
835 // By default, copy all input arguments since they will all be the same
836 // except for the given reduced p. We will then replace the reduced p with
837 // p_orig below.
838 MEB::InArgs<Scalar> wrappedInArgs = thyraModel->createInArgs();
839 wrappedInArgs.setArgs(inArgs);
840
841 // Replace p with p_orig.
842 RCP<const VectorBase<Scalar> > p;
843 if (!is_null(p=wrappedInArgs.get_p(p_idx_))) {
844 if (
845 dumpBasisMatrix_
846 && includesVerbLevel(localVerbLevel,Teuchos::VERB_MEDIUM)
847 )
848 {
849 *out << "\nB = " << Teuchos::describe(*B_,Teuchos::VERB_EXTREME);
850 }
851 wrappedInArgs.set_p(p_idx_,map_from_p_to_p_orig(*p));
852 }
853
854 //
855 // B) Setup OutArgs
856 //
857
858 // By default, copy all output arguments since they will all be the same
859 // except for those derivatives w.r.t. p(p_idx). We will then replace the
860 // derivative objects w.r.t. given reduced p with the derivarive objects
861 // w.r.t. p_orig below.
862 MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
863 wrappedOutArgs.setArgs(outArgs);
864
865 // Set derivative output arguments for p_orig if derivatives for p are
866 // reqeusted in outArgs
867 setupWrappedParamDerivOutArgs(outArgs,&wrappedOutArgs);
868
869 //
870 // C) Evaluate the underlying model functions
871 //
872
873 if (includesVerbLevel(localVerbLevel,Teuchos::VERB_LOW))
874 *out << "\nEvaluating the fully parameterized underlying model ...\n";
875 // Compute the underlying functions in terms of p_orig, including
876 // derivatives w.r.t. p_orig.
877 thyraModel->evalModel(wrappedInArgs,wrappedOutArgs);
878
879 //
880 // D) Postprocess the output arguments
881 //
882
883 // Assemble the derivatives for p given derivatives for p_orig computed
884 // above.
885 assembleParamDerivOutArgs(wrappedOutArgs,outArgs);
886
887 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
888
889}
890
891
892// private
893
894
895template<class Scalar>
896void DefaultLumpedParameterModelEvaluator<Scalar>::finishInitialization() const
897{
898
899 typedef ScalarTraits<Scalar> ST;
900 typedef ModelEvaluatorBase MEB;
901
902 if (isInitialized_)
903 return;
904
905 //
906 // A) Get the underlying model
907 //
908
909 const RCP<const ModelEvaluator<Scalar> >
910 thyraModel = this->getUnderlyingModel();
911
913 is_null(thyraModel), std::logic_error,
914 "Error, the underlying model evaluator must be set!" );
915
916 //
917 // B) Create B for the reduced affine model for the given parameter subvector
918 //
919
920 if (autogenerateBasisMatrix_) {
921 generateParameterBasisMatrix();
922 }
923 else {
925 true, std::logic_error,
926 "Error, we don't handle a client-set parameter basis matrix yet!" );
927 }
928
929 isInitialized_ = true;
930
931}
932
933
934template<class Scalar>
935void DefaultLumpedParameterModelEvaluator<Scalar>::generateParameterBasisMatrix() const
936{
937
938 using Teuchos::as;
939 typedef ScalarTraits<Scalar> ST;
940
941 const RCP<const ModelEvaluator<Scalar> >
942 thyraModel = this->getUnderlyingModel();
943
944 const RCP<const VectorSpaceBase<Scalar> >
945 p_orig_space = thyraModel->get_p_space(p_idx_);
946
947 const Ordinal p_orig_dim = p_orig_space->dim();
948
950 !( 1 <= numberOfBasisColumns_ && numberOfBasisColumns_ <= p_orig_dim ),
951 std::logic_error,
952 "Error, the number of basis columns = " << numberOfBasisColumns_ << " does not\n"
953 "fall in the range [1,"<<p_orig_dim<<"]!" );
954
955 //
956 // Create and randomize B
957 //
958 // Here we make the first column all ones and then randomize columns 1
959 // through numberOfBasisColumns_-1 so that the average entry is 1.0 with a
960 // spread of 1.0. This is just to give as well a scaled starting matrix as
961 // possible that will hopefully yeild a well scaled orthonomal B after we
962 // are finished.
963
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()),
970 B.ptr() );
971 }
972
973 //
974 // Create an orthogonal form of B using a modified Gram Schmidt algorithm
975 //
976
977 RCP<MultiVectorBase<double> > R;
978 sillyModifiedGramSchmidt( B.ptr(), Teuchos::outArg(R) );
979
980 // Above:
981 // 1) On output, B will have orthonomal columns which makes it a good basis
982 // 2) We just discard the "R" factor since we don't need it for anything
983
984 B_ = B;
985
986}
987
988
989template<class Scalar>
990void DefaultLumpedParameterModelEvaluator<Scalar>::updateNominalValuesAndBounds() const
991{
992
993 typedef ScalarTraits<Scalar> ST;
994 typedef ModelEvaluatorBase MEB;
995
996 if (nominalValuesAndBoundsUpdated_)
997 return;
998
999 finishInitialization();
1000
1001 const RCP<const ModelEvaluator<Scalar> >
1002 thyraModel = this->getUnderlyingModel();
1003
1004 const MEB::InArgs<Scalar> origNominalValues = thyraModel->getNominalValues();
1005 const MEB::InArgs<Scalar> origLowerBounds = thyraModel->getLowerBounds();
1006 const MEB::InArgs<Scalar> origUpperBounds = thyraModel->getUpperBounds();
1007
1008 // p_orig_base
1009
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();
1018 }
1019 else {
1021 true, std::logic_error,
1022 "Error, we don't handle reading in the parameter base vector yet!" );
1023 }
1024
1025 // Nominal values
1026
1027 nominalValues_ = origNominalValues;
1028
1029 if (nominalValueIsParameterBase_) {
1030 // A value of p==0 will give p_orig = p_orig_init!
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);
1035 }
1036 else {
1038 true, std::logic_error,
1039 "Error, we don't handle creating p_init when p_orig_base != p_orig_init yet!" );
1040 }
1041
1042 // Bounds
1043
1044 lowerBounds_ = origLowerBounds;
1045 upperBounds_ = origUpperBounds;
1046
1047 lowerBounds_.set_p(p_idx_,Teuchos::null);
1048 upperBounds_.set_p(p_idx_,Teuchos::null);
1049
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_);
1054 if ( !is_null(p_orig_l) || !is_null(p_orig_u) ) {
1056 true, std::logic_error,
1057 "Error, we don't handle bounds on p_orig yet!" );
1058 }
1059 }
1060
1061 nominalValuesAndBoundsUpdated_ = true;
1062
1063}
1064
1065
1066template<class Scalar>
1067RCP<VectorBase<Scalar> >
1068DefaultLumpedParameterModelEvaluator<Scalar>::map_from_p_to_p_orig(
1069 const VectorBase<Scalar> &p
1070 ) const
1071{
1072 // p_orig = B*p + p_orig_base
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_ );
1076 return p_orig;
1077}
1078
1079
1080template<class Scalar>
1081void DefaultLumpedParameterModelEvaluator<Scalar>::setupWrappedParamDerivOutArgs(
1082 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs, // in
1083 ModelEvaluatorBase::OutArgs<Scalar> *wrappedOutArgs_inout // in/out
1084 ) const
1085{
1086
1087 typedef ModelEvaluatorBase MEB;
1088 typedef MEB::Derivative<Scalar> Deriv;
1089
1090 TEUCHOS_TEST_FOR_EXCEPT(wrappedOutArgs_inout==0);
1091 MEB::OutArgs<Scalar> &wrappedOutArgs = *wrappedOutArgs_inout;
1092
1093 Deriv DfDp;
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));
1096 }
1097
1098 const int Ng = outArgs.Ng();
1099 for ( int j = 0; j < Ng; ++j ) {
1100 Deriv DgDp;
1101 if ( !(DgDp=outArgs.get_DgDp(j,p_idx_)).isEmpty() ) {
1102 wrappedOutArgs.set_DgDp(
1103 j, p_idx_,
1104 create_deriv_wrt_p_orig(DgDp,DgDp.getMultiVectorOrientation())
1105 );
1106 }
1107 }
1108
1109}
1110
1111
1112template<class Scalar>
1113ModelEvaluatorBase::Derivative<Scalar>
1114DefaultLumpedParameterModelEvaluator<Scalar>::create_deriv_wrt_p_orig(
1115 const ModelEvaluatorBase::Derivative<Scalar> &DhDp,
1117 ) const
1118{
1119
1120 typedef ModelEvaluatorBase MEB;
1121
1122 const RCP<const MultiVectorBase<Scalar> >
1123 DhDp_mv = DhDp.getMultiVector();
1125 is_null(DhDp_mv) || (DhDp.getMultiVectorOrientation() != requiredOrientation),
1126 std::logic_error,
1127 "Error, we currently can't handle non-multi-vector derivatives!" );
1128
1129 RCP<MultiVectorBase<Scalar> > DhDp_orig_mv;
1130 switch (requiredOrientation) {
1131 case MEB::DERIV_MV_BY_COL:
1132 // DhDp = DhDp_orig * B
1133 DhDp_orig_mv = createMembers(DhDp_mv->range(),B_->range()->dim());
1134 // Above, we could just request DhDp_orig as a LinearOpBase object since
1135 // we just need to apply it!
1136 break;
1137 case MEB::DERIV_TRANS_MV_BY_ROW:
1138 // (DhDp^T) = B^T * (DhDp_orig^T) [DhDp_orig_mv is transposed!]
1139 DhDp_orig_mv = createMembers(B_->range(),DhDp_mv->domain()->dim());
1140 // Above, we really do need DhDp_orig as the gradient form multi-vector
1141 // since it must be the RHS for a linear operator apply!
1142 break;
1143 default:
1144 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
1145 }
1146
1147 return MEB::Derivative<Scalar>(DhDp_orig_mv,requiredOrientation);
1148
1149}
1150
1151
1152template<class Scalar>
1153void DefaultLumpedParameterModelEvaluator<Scalar>::assembleParamDerivOutArgs(
1154 const ModelEvaluatorBase::OutArgs<Scalar> &wrappedOutArgs, // in
1155 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs // in/out
1156 ) const
1157{
1158
1159 typedef ModelEvaluatorBase MEB;
1160 typedef MEB::Derivative<Scalar> Deriv;
1161
1162 Deriv DfDp;
1163 if ( !(DfDp=outArgs.get_DfDp(p_idx_)).isEmpty() ) {
1164 assembleParamDeriv( wrappedOutArgs.get_DfDp(p_idx_), DfDp );
1165 }
1166
1167 const int Ng = outArgs.Ng();
1168 for ( int j = 0; j < Ng; ++j ) {
1169 Deriv DgDp;
1170 if ( !(DgDp=outArgs.get_DgDp(j,p_idx_)).isEmpty() ) {
1171 assembleParamDeriv( wrappedOutArgs.get_DgDp(j,p_idx_), DgDp );
1172 }
1173 }
1174
1175}
1176
1177
1178template<class Scalar>
1179void DefaultLumpedParameterModelEvaluator<Scalar>::assembleParamDeriv(
1180 const ModelEvaluatorBase::Derivative<Scalar> &DhDp_orig, // in
1181 const ModelEvaluatorBase::Derivative<Scalar> &DhDp // in/out
1182 ) const
1183{
1184
1185 typedef ModelEvaluatorBase MEB;
1186
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!" );
1192
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!" );
1198
1199 switch( DhDp_orig.getMultiVectorOrientation() ) {
1200 case MEB::DERIV_MV_BY_COL:
1201 // DhDp = DhDp_orig * B
1202#ifdef TEUCHSO_DEBUG
1204 DhDp.getMultiVectorOrientation() == MEB::DERIV_MV_BY_COL );
1205#endif
1206 apply( *DhDp_orig_mv, NOTRANS, *B_, DhDp_mv.ptr() );
1207 // Above, we could generalize DhDp_oirg to just be a general linear
1208 // operator.
1209 break;
1210 case MEB::DERIV_TRANS_MV_BY_ROW:
1211 // (DhDp^T) = B^T * (DhDp_orig^T) [DhDp_orig_mv is transposed!]
1212#ifdef TEUCHSO_DEBUG
1214 DhDp.getMultiVectorOrientation() == MEB::DERIV_TRANS_MV_BY_ROW );
1215#endif
1216 apply( *B_, CONJTRANS, *DhDp_orig_mv, DhDp_mv.ptr() );
1217 break;
1218 default:
1219 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
1220 }
1221
1222}
1223
1224
1225} // namespace Thyra
1226
1227
1228#endif // THYRA_DEFAUL_LUMPED_PARAMETER_LUMPED_MODEL_EVALUATOR_HPP
T * get() const
Decorator class that wraps any ModelEvaluator object and lumps parameters together using a linear bas...
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
.
RCP< DefaultLumpedParameterModelEvaluator< Scalar > > defaultLumpedParameterModelEvaluator(const RCP< ModelEvaluator< Scalar > > &thyraModel)
Non-member constructor.
RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) 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.
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 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)
T_To & dyn_cast(T_From &from)
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)