Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
58 namespace Thyra {
59 
60 
239 template<class Scalar>
241  : virtual public ModelEvaluatorDelegatorBase<Scalar>
242  , virtual public Teuchos::ParameterListAcceptor
243 {
244 public:
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 
312 private:
313 
316 
318  ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
320  void evalModelImpl(
323  ) const;
324 
326 
327 private:
328 
329  // ////////////////////////////////
330  // Private data members
331 
332  mutable bool isInitialized_;
333  mutable bool nominalValuesAndBoundsUpdated_;
334 
335  mutable RCP<const Teuchos::ParameterList> validParamList_;
336  RCP<Teuchos::ParameterList> paramList_;
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 
431 template<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 
451 template<class Scalar>
452 const std::string
454 = "Parameter Subvector Index";
455 
456 template<class Scalar>
457 const int
459 = 0;
460 
461 
462 template<class Scalar>
463 const std::string
465 = "Auto-generate Basis Matrix";
466 
467 template<class Scalar>
468 const bool
470 = true;
471 
472 
473 template<class Scalar>
474 const std::string
476 = "Number of Basis Columns";
477 
478 template<class Scalar>
479 const int
481 = 1;
482 
483 
484 template<class Scalar>
485 const std::string
487 = "Nominal Value is Parameter Base";
488 
489 template<class Scalar>
490 const bool
492 = true;
493 
494 
495 template<class Scalar>
496 const std::string
498 = "Parameter Base Vector";
499 
500 
501 template<class Scalar>
502 const std::string
504 = "Ignore Parameter Bounds";
505 
506 template<class Scalar>
507 const bool
509 = false;
510 
511 
512 template<class Scalar>
513 const std::string
515 = "Dump Basis Matrix";
516 
517 template<class Scalar>
518 const bool
520 = false;
521 
522 
523 // Constructors/initializers/accessors/utilities
524 
525 
526 template<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 
540 template<class Scalar>
542  const RCP<ModelEvaluator<Scalar> > &thyraModel
543  )
544 {
545  isInitialized_ = false;
546  nominalValuesAndBoundsUpdated_ = false;
548 }
549 
550 
551 template<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 
565 template<class Scalar>
566 std::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 
586 template<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 
634 template<class Scalar>
637 {
638  return paramList_;
639 }
640 
641 
642 template<class Scalar>
645 {
646  RCP<Teuchos::ParameterList> _paramList = paramList_;
647  paramList_ = Teuchos::null;
648  return _paramList;
649 }
650 
651 
652 template<class Scalar>
655 {
656  return paramList_;
657 }
658 
659 
660 template<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 
706 template<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 
717 template<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 
728 template<class Scalar>
731 {
732  updateNominalValuesAndBounds();
733  return nominalValues_;
734 }
735 
736 
737 template<class Scalar>
740 {
741  updateNominalValuesAndBounds();
742  return lowerBounds_;
743 }
744 
745 
746 template<class Scalar>
749 {
750  updateNominalValuesAndBounds();
751  return upperBounds_;
752 }
753 
754 
755 template<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 
790 template<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 
804 template<class Scalar>
805 void 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 
895 template<class Scalar>
896 void 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 
934 template<class Scalar>
935 void 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 
989 template<class Scalar>
990 void 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 
1066 template<class Scalar>
1067 RCP<VectorBase<Scalar> >
1068 DefaultLumpedParameterModelEvaluator<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 
1080 template<class Scalar>
1081 void 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 
1112 template<class Scalar>
1113 ModelEvaluatorBase::Derivative<Scalar>
1114 DefaultLumpedParameterModelEvaluator<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 
1152 template<class Scalar>
1153 void 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 
1178 template<class Scalar>
1179 void 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
void reportFinalPoint(const ModelEvaluatorBase::InArgs< Scalar > &finalPoint, const bool wasSolved)
Pure abstract base interface for evaluating a stateless &quot;model&quot; that can be mapped into a number of d...
Decorator class that wraps any ModelEvaluator object and lumps parameters together using a linear bas...
bool is_null(const boost::shared_ptr< T > &p)
basic_OSTab< char > OSTab
RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...
RCP< const Array< std::string > > get_p_names(int l) const
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
This is a base class that delegetes almost all function to a wrapped model evaluator object...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Use the non-transposed operator.
void initialize(const RCP< ModelEvaluator< Scalar > > &model)
Initialize given a non-const model evaluator.
Simple aggregate class that stores a derivative object as a general linear operator or as a multi-vec...
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
T * get() const
void setParameterList(RCP< Teuchos::ParameterList > const &paramList)
ModelEvaluatorBase::InArgs< Scalar > getNominalValues() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Abstract interface for finite-dimensional dense vectors.
void initialize(const RCP< ModelEvaluator< Scalar > > &thyraModel)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Base subclass for ModelEvaluator that defines some basic types.
RCP< DefaultLumpedParameterModelEvaluator< Scalar > > defaultLumpedParameterModelEvaluator(const RCP< ModelEvaluator< Scalar > > &thyraModel)
Non-member constructor.
TypeTo as(const TypeFrom &t)
#define TEUCHOS_ASSERT(assertion_test)
void setModelEvalDescription(const std::string &modelEvalDescription)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 &lt;= l &amp;&amp; l &lt; this-&gt;Np().
Protected subclass of OutArgs that only ModelEvaluator subclasses can access to set up the selection ...
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object...