MoochoPack_MoochoThyraSolver.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "MoochoPack_MoochoThyraSolver.hpp"
00030 #include "NLPInterfacePack_NLPDirectThyraModelEvaluator.hpp"
00031 #include "NLPInterfacePack_NLPFirstOrderThyraModelEvaluator.hpp"
00032 #include "Thyra_DefaultFiniteDifferenceModelEvaluator.hpp"
00033 #include "Thyra_DefaultStateEliminationModelEvaluator.hpp"
00034 #include "Thyra_DefaultEvaluationLoggerModelEvaluator.hpp"
00035 #include "Thyra_DefaultInverseModelEvaluator.hpp"
00036 #include "Thyra_DefaultLumpedParameterModelEvaluator.hpp"
00037 #include "Thyra_DefaultSpmdMultiVectorFileIO.hpp"
00038 #include "Thyra_DampenedNewtonNonlinearSolver.hpp"
00039 #include "Thyra_ModelEvaluatorHelpers.hpp"
00040 #include "Thyra_VectorStdOps.hpp"
00041 #include "Teuchos_StandardParameterEntryValidators.hpp"
00042 #include "Teuchos_XMLParameterListHelpers.hpp"
00043 
00044 namespace {
00045 
00046 //
00047 // ParameterList parameters and sublists
00048 //
00049 
00050 const std::string SolveMode_name = "Solve Mode";
00051 const Teuchos::RCP<
00052   Teuchos::StringToIntegralParameterEntryValidator<
00053     MoochoPack::MoochoThyraSolver::ESolveMode
00054   >
00055 >
00056 solveModeValidator = Teuchos::rcp(
00057   new Teuchos::StringToIntegralParameterEntryValidator<MoochoPack::MoochoThyraSolver::ESolveMode>(
00058     Teuchos::tuple<std::string>(
00059       "Forward Solve"
00060       ,"Optimize"
00061       )
00062     ,Teuchos::tuple<std::string>(
00063       "Only solve state equaitons f(x,p)=0 for states x\n"
00064       "given fixed parameters values p."
00065       ,"Solve the simulation constrained optimization problem\n"
00066       "  min  g(x,p)\n"
00067       "  s.t. f(x,p)=0\n"
00068       "for the state varaibles x and parameters p."
00069       )
00070     ,Teuchos::tuple<MoochoPack::MoochoThyraSolver::ESolveMode>(
00071       MoochoPack::MoochoThyraSolver::SOLVE_MODE_FORWARD
00072       ,MoochoPack::MoochoThyraSolver::SOLVE_MODE_OPTIMIZE
00073       )
00074     ,""
00075     )
00076   );
00077 const std::string SolveMode_default = "Optimize";
00078 
00079 const std::string NLPType_name = "NLP Type";
00080 const Teuchos::RCP<
00081   Teuchos::StringToIntegralParameterEntryValidator<
00082     MoochoPack::MoochoThyraSolver::ENLPType
00083     >
00084   >
00085 nlpTypeValidator = Teuchos::rcp(
00086   new Teuchos::StringToIntegralParameterEntryValidator<MoochoPack::MoochoThyraSolver::ENLPType>(
00087     Teuchos::tuple<std::string>(
00088       "First Order"
00089       ,"Direct"
00090       )
00091     ,Teuchos::tuple<std::string>(
00092       "Support the NLPInterfacePack::NLPFirstOrder interface which assumes\n"
00093       "that full adjoints for the objective and constraint derivatives are\n"
00094       "available."
00095       ,"Support the NLPInterfacePack::NLPDirect interface which only assumes\n"
00096       "that forward or direct sensitivities and state solves are supported."
00097       )
00098     ,Teuchos::tuple<MoochoPack::MoochoThyraSolver::ENLPType>(
00099       MoochoPack::MoochoThyraSolver::NLP_TYPE_FIRST_ORDER
00100       ,MoochoPack::MoochoThyraSolver::NLP_TYPE_DIRECT
00101       )
00102     ,""
00103     )
00104   );
00105 const std::string NLPType_default = "First Order";
00106 
00107 const std::string NonlinearlyEliminateStates_name = "Nonlinearly Eliminate States";
00108 const bool NonlinearlyEliminateStates_default = false;
00109 
00110 const std::string UseFiniteDifferencesForObjective_name = "Use Finite Differences For Objective";
00111 const bool UseFiniteDifferencesForObjective_default = false;
00112 
00113 const std::string ObjectiveFiniteDifferenceSettings_name = "Objective Finite Difference Settings";
00114 
00115 const std::string UseFiniteDifferencesForConstraints_name = "Use Finite Differences For Constraints";
00116 const bool UseFiniteDifferencesForConstraints_default = false;
00117 
00118 const std::string ConstraintsFiniteDifferenceSettings_name = "Constraints Finite Difference Settings";
00119 
00120 const std::string FwdNewtonTol_name = "Forward Newton Tolerance";
00121 const double FwdNewtonTol_default = -1.0;
00122 
00123 const std::string FwdNewtonMaxIters_name = "Forward Newton Max Iters";
00124 const int FwdNewtonMaxIters_default = 20;
00125 
00126 const std::string ForwardNewtonDampening_name = "Forward Newton Dampening";
00127 const bool ForwardNewtonDampening_default = true;
00128 
00129 const std::string FwdNewtonMaxLineSearchIters_name = "Forward Newton Max Line Search Iters";
00130 const int FwdNewtonMaxLineSearchIters_default = 20;
00131 
00132 const std::string UseBuiltInInverseObjectiveFunction_name = "Use Built-in Inverse Objective Function";
00133 const bool UseBuiltInInverseObjectiveFunction_default = false;
00134 
00135 const std::string InverseObjectiveFunctionSettings_name = "Inverse Objective Function Settings";
00136 
00137 const std::string UseParameterLumping_name = "Use Parameter Lumping";
00138 const bool UseParameterLumping_default = false;
00139 
00140 const std::string LumpedParametersSettings_name = "Lumped Parameters Settings";
00141 
00142 const std::string OutputFileTag_name = "Output File Tag";
00143 const std::string OutputFileTag_default = "";
00144 
00145 const std::string ShowModelEvaluatorTrace_name = "Show Model Evaluator Trace";
00146 const bool ShowModelEvaluatorTrace_default = "false";
00147 
00148 const std::string StateGuess_name = "State Guess";
00149 
00150 const std::string ParamGuess_name = "Parameter Guess";
00151 
00152 const std::string ParamLowerBounds_name = "Parameter Lower Bounds";
00153 
00154 const std::string ParamUpperBounds_name = "Parameter Upper Bounds";
00155 
00156 const std::string StateSoluFileBaseName_name = "State Solution File Base Name";
00157 const std::string StateSoluFileBaseName_default = "";
00158 
00159 const std::string ParamSoluFileBaseName_name = "Parameters Solution File Base Name";
00160 const std::string ParamSoluFileBaseName_default = "";
00161 
00162 } // namespace
00163 
00164 namespace MoochoPack {
00165 
00166 // Constructors/initialization
00167 
00168 MoochoThyraSolver::MoochoThyraSolver(
00169   const std::string    &paramsXmlFileName
00170   ,const std::string   &extraParamsXmlString
00171   ,const std::string   &paramsUsedXmlOutFileName
00172   ,const std::string   &paramsXmlFileNameOption
00173   ,const std::string   &extraParamsXmlStringOption
00174   ,const std::string   &paramsUsedXmlOutFileNameOption
00175   )
00176   :paramsXmlFileName_(paramsXmlFileName)
00177   ,extraParamsXmlString_(extraParamsXmlString)
00178   ,paramsUsedXmlOutFileName_(paramsUsedXmlOutFileName)
00179   ,paramsXmlFileNameOption_(paramsXmlFileNameOption)
00180   ,extraParamsXmlStringOption_(extraParamsXmlStringOption)
00181   ,paramsUsedXmlOutFileNameOption_(paramsUsedXmlOutFileNameOption)
00182   ,stateVectorIO_(Teuchos::rcp(new Thyra::DefaultSpmdMultiVectorFileIO<value_type>))
00183   ,parameterVectorIO_(Teuchos::rcp(new Thyra::DefaultSpmdMultiVectorFileIO<value_type>))
00184   ,solveMode_(SOLVE_MODE_OPTIMIZE)
00185   ,nlpType_(NLP_TYPE_FIRST_ORDER)
00186   ,nonlinearlyElimiateStates_(false)
00187   ,use_finite_diff_for_obj_(false)
00188   ,use_finite_diff_for_con_(false)
00189   ,fwd_newton_tol_(-1.0)
00190   ,fwd_newton_max_iters_(20)
00191   ,fwd_newton_dampening_(false)
00192   ,fwd_newton_max_ls_iters_(20)
00193   ,useInvObjFunc_(false)
00194   ,useParameterLumping_(false)
00195   ,outputFileTag_("")
00196   ,showModelEvaluatorTrace_(false)
00197   ,stateSoluFileBase_("")
00198   ,paramSoluFileBase_("")
00199 {}
00200 
00201 MoochoThyraSolver::~MoochoThyraSolver()
00202 {}
00203 
00204 void MoochoThyraSolver::setupCLP(
00205   Teuchos::CommandLineProcessor *clp
00206   )
00207 {
00208   TEST_FOR_EXCEPT(0==clp);
00209   solver_.setup_commandline_processor(clp);
00210   clp->setOption(
00211     paramsXmlFileNameOption().c_str(),&paramsXmlFileName_
00212     ,"Name of an XML file containing parameters for linear solver options to be appended first."
00213     );
00214   clp->setOption(
00215     extraParamsXmlStringOption().c_str(),&extraParamsXmlString_
00216     ,"An XML string containing linear solver parameters to be appended second."
00217     );
00218   clp->setOption(
00219     paramsUsedXmlOutFileNameOption().c_str(),&paramsUsedXmlOutFileName_
00220     ,"Name of an XML file that can be written with the parameter list after it has been used on completion of this program."
00221     );
00222 }
00223 
00224 void MoochoThyraSolver::readParameters( std::ostream *out_arg )
00225 {
00226   Teuchos::RCP<Teuchos::FancyOStream>
00227     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00228   Teuchos::OSTab tab(out);
00229   if(out.get()) *out << "\nMoochoThyraSolver::readParameters(...):\n";
00230   Teuchos::OSTab tab2(out);
00231   Teuchos::RCP<Teuchos::ParameterList>
00232     paramList = this->getParameterList();
00233   if(!paramList.get()) {
00234     if(out.get()) *out << "\nCreating a new Teuchos::ParameterList ...\n";
00235     paramList = Teuchos::rcp(new Teuchos::ParameterList("MoochoThyraSolver"));
00236   }
00237   if(paramsXmlFileName().length()) {
00238     if(out.get()) *out << "\nReading parameters from XML file \""<<paramsXmlFileName()<<"\" ...\n";
00239     Teuchos::updateParametersFromXmlFile(paramsXmlFileName(),&*paramList);
00240   }
00241   if(extraParamsXmlString().length()) {
00242     if(out.get())
00243       *out << "\nAppending extra parameters from the XML string \""<<extraParamsXmlString()<<"\" ...\n";
00244     Teuchos::updateParametersFromXmlString(extraParamsXmlString(),&*paramList);
00245   }
00246   if( paramsXmlFileName().length() || extraParamsXmlString().length() ) {
00247     typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
00248     if(out.get()) {
00249       *out  << "\nUpdated parameter list:\n";
00250       paramList->print(
00251         *out,PLPrintOptions().indent(2).showTypes(true)
00252         );
00253     }
00254     this->setParameterList(paramList);
00255   }
00256 }
00257 
00258 // Overridden from ParameterListAcceptor
00259 
00260 void MoochoThyraSolver::setParameterList(
00261   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00262   )
00263 {
00264   TEST_FOR_EXCEPT(!paramList.get());
00265   paramList->validateParameters(*getValidParameters(),0); // Just validate my level!
00266   paramList_ = paramList;
00267   solveMode_ = solveModeValidator->getIntegralValue(
00268     *paramList_,SolveMode_name,SolveMode_default);
00269   nlpType_ = nlpTypeValidator->getIntegralValue(
00270     *paramList_,NLPType_name,NLPType_default);
00271   nonlinearlyElimiateStates_ = paramList_->get(
00272     NonlinearlyEliminateStates_name,NonlinearlyEliminateStates_default);
00273   use_finite_diff_for_obj_ = paramList_->get(
00274     UseFiniteDifferencesForObjective_name,UseFiniteDifferencesForObjective_default);
00275   use_finite_diff_for_con_ = paramList_->get(
00276     UseFiniteDifferencesForConstraints_name,UseFiniteDifferencesForConstraints_default);
00277   fwd_newton_tol_ = paramList_->get(
00278     FwdNewtonTol_name,FwdNewtonTol_default);
00279   fwd_newton_max_iters_ = paramList_->get(
00280     FwdNewtonMaxIters_name,FwdNewtonMaxIters_default);
00281   fwd_newton_dampening_ = paramList_->get(
00282     ForwardNewtonDampening_name,ForwardNewtonDampening_default);
00283   fwd_newton_max_ls_iters_ = paramList_->get(
00284     FwdNewtonMaxLineSearchIters_name,FwdNewtonMaxLineSearchIters_default);
00285   useInvObjFunc_ = paramList_->get(
00286     UseBuiltInInverseObjectiveFunction_name,UseBuiltInInverseObjectiveFunction_default);
00287   useParameterLumping_ = paramList_->get(
00288     UseParameterLumping_name, UseParameterLumping_default);
00289   outputFileTag_ = paramList->get(
00290     OutputFileTag_name,OutputFileTag_default);
00291   solver_.set_output_file_tag(outputFileTag_);
00292   showModelEvaluatorTrace_ = paramList->get(
00293     ShowModelEvaluatorTrace_name,ShowModelEvaluatorTrace_default);
00294   x_reader_.setParameterList(
00295     sublist(paramList_,StateGuess_name)
00296     );
00297   p_reader_.setParameterList(
00298     sublist(paramList_,ParamGuess_name)
00299     );
00300   p_l_reader_.setParameterList(
00301     sublist(paramList_,ParamLowerBounds_name)
00302     );
00303   p_u_reader_.setParameterList(
00304     sublist(paramList_,ParamUpperBounds_name)
00305     );
00306   stateSoluFileBase_ = paramList_->get(
00307     StateSoluFileBaseName_name,StateSoluFileBaseName_default);
00308   paramSoluFileBase_ = paramList_->get(
00309     ParamSoluFileBaseName_name,ParamSoluFileBaseName_default);
00310 #ifdef TEUCHOS_DEBUG
00311   paramList->validateParameters(*getValidParameters(),0); // Just validate my level!
00312 #endif
00313 }
00314 
00315 Teuchos::RCP<Teuchos::ParameterList>
00316 MoochoThyraSolver::getParameterList()
00317 {
00318   return paramList_;
00319 }
00320 
00321 Teuchos::RCP<Teuchos::ParameterList>
00322 MoochoThyraSolver::unsetParameterList()
00323 {
00324   Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
00325   paramList_ = Teuchos::null;
00326   return _paramList;
00327 }
00328 
00329 Teuchos::RCP<const Teuchos::ParameterList>
00330 MoochoThyraSolver::getParameterList() const
00331 {
00332   return paramList_;
00333 }
00334 
00335 Teuchos::RCP<const Teuchos::ParameterList>
00336 MoochoThyraSolver::getValidParameters() const
00337 {
00338   static Teuchos::RCP<Teuchos::ParameterList> pl;
00339   if(pl.get()==NULL) {
00340     pl = Teuchos::rcp(new Teuchos::ParameterList());
00341     pl->set(
00342       SolveMode_name,SolveMode_default
00343       ,"The type of solve to perform."
00344       ,solveModeValidator
00345       );
00346     pl->set(
00347       NLPType_name,NLPType_default
00348       ,"The type of MOOCHO NLP subclass to use."
00349       ,nlpTypeValidator
00350       );
00351     pl->set(
00352       NonlinearlyEliminateStates_name,NonlinearlyEliminateStates_default
00353       ,"If true, then the model's state equations and state variables\n"
00354       "are nonlinearlly eliminated using a forward solver."
00355       );
00356     pl->set(
00357       UseFiniteDifferencesForObjective_name,UseFiniteDifferencesForObjective_default
00358       ,"Use finite differences for missing objective function derivatives (Direct NLP only).\n"
00359       "See the options in the sublist \"" + ObjectiveFiniteDifferenceSettings_name + "\"."
00360       );
00361     {
00362       Thyra::DirectionalFiniteDiffCalculator<Scalar> dfdcalc;
00363       {
00364         Teuchos::ParameterList
00365           &fdSublist = pl->sublist(ObjectiveFiniteDifferenceSettings_name);
00366         fdSublist.setParameters(*dfdcalc.getValidParameters());
00367       }
00368       pl->set(
00369         UseFiniteDifferencesForConstraints_name,UseFiniteDifferencesForConstraints_default
00370         ,"Use  finite differences for missing constraint derivatives (Direct NLP only).\n"
00371         "See the   options in the sublist \"" + ConstraintsFiniteDifferenceSettings_name + "\"."
00372         );
00373       {
00374         Teuchos::ParameterList
00375           &fdSublist = pl->sublist(ConstraintsFiniteDifferenceSettings_name);
00376         fdSublist.setParameters(*dfdcalc.getValidParameters());
00377       }
00378     }
00379     pl->set(
00380       FwdNewtonTol_name,FwdNewtonTol_default
00381       ,"Tolarance used for the forward state solver in eliminating\n"
00382       "the state equations/variables."
00383       );
00384     pl->set(
00385       FwdNewtonMaxIters_name,FwdNewtonMaxIters_default
00386       ,"Maximum number of iterations allowed for the forward state\n"
00387       "solver in eliminating the state equations/variables."
00388       );
00389     pl->set(
00390       ForwardNewtonDampening_name, ForwardNewtonDampening_default,
00391       "If true, then the state elimination nonlinear solver will\n"
00392       "use a dampened line search.  Otherwise, it will just take fulls steps."
00393       );
00394     pl->set(
00395       FwdNewtonMaxLineSearchIters_name, FwdNewtonMaxLineSearchIters_default,
00396       "Maximum number of linea search iterations per newton iteration\n"
00397       "allowed for the forward state solver in eliminating the state equations/variables."
00398       );
00399     pl->set(
00400       UseBuiltInInverseObjectiveFunction_name,UseBuiltInInverseObjectiveFunction_default
00401       ,"Use a built-in form of a simple inverse objection function instead\n"
00402       "of a a response function contained in the underlying model evaluator\n"
00403       "object itself.  The settings are contained in the sublist\n"
00404       "\""+InverseObjectiveFunctionSettings_name+"\".\n"
00405       "Note that this feature allows the client to form a useful type\n"
00406       "of optimization problem just with a model that supports only the\n"
00407       "parameterized state function f(x,p)=0."
00408       );
00409     {
00410       Teuchos::RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00411         inverseModel = rcp(new Thyra::DefaultInverseModelEvaluator<Scalar>());
00412       pl->sublist(
00413         InverseObjectiveFunctionSettings_name,false
00414         ,"Settings for the built-in inverse objective function.\n"
00415         "See the outer parameter \""+UseBuiltInInverseObjectiveFunction_name+"\"."
00416         ).setParameters(*inverseModel->getValidParameters());
00417     }
00418     pl->set(
00419       UseParameterLumping_name,UseParameterLumping_default,
00420       "Use a reduced basis to lump optimization parameters as\n"
00421       "p_orig = P_basis * p.  If set to true, then the settings\n"
00422       "in \""+LumpedParametersSettings_name+"\" determine how the\n"
00423       "parameter basis is set.  This feature can be used to safely\n"
00424       "regularize a problem if there are linearly dependent parameters\n"
00425       "and will generally speed up the optimiztation algorithms."
00426       );
00427     {
00428       Teuchos::RCP<Thyra::DefaultLumpedParameterModelEvaluator<Scalar> >
00429         lumpedParamModel = rcp(new Thyra::DefaultLumpedParameterModelEvaluator<Scalar>);
00430       pl->sublist(
00431         LumpedParametersSettings_name,false
00432         ,"Settings for parameter lumping.\n"
00433         "See the outer parameter \""+UseParameterLumping_name+"\"."
00434         ).setParameters(*lumpedParamModel->getValidParameters());
00435     }
00436     pl->set(OutputFileTag_name,OutputFileTag_default,
00437       "A tag that is attached to every output file that is created by the\n"
00438       "solver.  If empty \"\", then no tag is used.  This option simply is\n"
00439       "passed into the set_output_context(...) function on the underlying\n"
00440       "MoochoPack::MoochoSolver object.  Therefore, this same parameter\n"
00441       "can be set in code as well without going through the parameter list.");
00442     pl->set(ShowModelEvaluatorTrace_name,ShowModelEvaluatorTrace_default
00443       ,"Determine if a trace of the objective function will be shown or not\n"
00444       "when the NLP is evaluated."
00445       );
00446     if(this->get_stateVectorIO().get())
00447       x_reader_.set_fileIO(this->get_stateVectorIO());
00448     pl->sublist(StateGuess_name).setParameters(*x_reader_.getValidParameters());
00449     if(this->get_parameterVectorIO().get()) {
00450       p_reader_.set_fileIO(this->get_parameterVectorIO());
00451       p_l_reader_.set_fileIO(this->get_parameterVectorIO());
00452       p_u_reader_.set_fileIO(this->get_parameterVectorIO());
00453       pl->sublist(ParamGuess_name).setParameters(*p_reader_.getValidParameters());
00454       pl->sublist(ParamLowerBounds_name).setParameters(*p_l_reader_.getValidParameters());
00455       pl->sublist(ParamUpperBounds_name).setParameters(*p_u_reader_.getValidParameters());
00456     }
00457     pl->set(
00458       StateSoluFileBaseName_name,StateSoluFileBaseName_default
00459       ,"If specified, a file with this basename will be written to with\n"
00460       "the final value of the state variables.  A different file for each\n"
00461       "process will be created.  Note that these files can be used for the\n"
00462       "initial guess for the state variables."
00463       );
00464     pl->set(
00465       ParamSoluFileBaseName_name,ParamSoluFileBaseName_default
00466       ,"If specified, a file with this basename will be written to with\n"
00467       "the final value of the parameters.  A different file for each\n"
00468       "process will be created.  Note that these files can be used for the\n"
00469       "initial guess for the parameters."
00470       );
00471   }
00472   return pl;
00473 }
00474 
00475 // Misc Access/Setup
00476 
00477 void MoochoThyraSolver::setSolveMode( const ESolveMode solveMode )
00478 {
00479   solveMode_ = solveMode;
00480 }
00481 
00482 MoochoThyraSolver::ESolveMode
00483 MoochoThyraSolver::getSolveMode() const
00484 {
00485   return solveMode_;
00486 }
00487 
00488 MoochoSolver& MoochoThyraSolver::getSolver()
00489 {
00490   return solver_;
00491 }
00492 
00493 const MoochoSolver& MoochoThyraSolver::getSolver() const
00494 {
00495   return solver_;
00496 }
00497 
00498 // Model specification, setup, solve, and solution extraction.
00499 
00500 void MoochoThyraSolver::setModel(
00501   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> > &origModel,
00502   const int p_idx,
00503   const int g_idx
00504   )
00505 {
00506 
00507   using Teuchos::rcp;
00508   using Teuchos::sublist;
00509   using NLPInterfacePack::NLP;
00510   using NLPInterfacePack::NLPDirectThyraModelEvaluator;
00511   using NLPInterfacePack::NLPFirstOrderThyraModelEvaluator;
00512 
00513   origModel_ = origModel;
00514   p_idx_ = p_idx;
00515   g_idx_ = g_idx;
00516 
00517   const int procRank = Teuchos::GlobalMPISession::getRank();
00518   //const int numProcs = Teuchos::GlobalMPISession::getNProc();
00519 
00520   //
00521   // Wrap the orginal model in different decorators
00522   //
00523 
00524   outerModel_ = origModel_;
00525 
00526   // Inverse objective (and parameter regularization)?
00527   if (useInvObjFunc_) {
00528     Teuchos::RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00529       inverseModel = Thyra::defaultInverseModelEvaluator<Scalar>(outerModel_);
00530     inverseModel->setVerbLevel(Teuchos::VERB_LOW);
00531     inverseModel->set_observationTargetIO(get_stateVectorIO());
00532     inverseModel->set_parameterBaseIO(get_parameterVectorIO());
00533     inverseModel->setParameterList(
00534       Teuchos::sublist(paramList_,InverseObjectiveFunctionSettings_name) );
00535     outerModel_ = inverseModel; 
00536     g_idx_ = inverseModel->Ng()-1;
00537   }
00538   
00539   // Evaluation loging
00540   Teuchos::RCP<std::ostream>
00541     modelEvalLogOut = Teuchos::fancyOStream(
00542       solver_.generate_output_file("ModelEvaluationLog")
00543       );
00544   Teuchos::RCP<Thyra::DefaultEvaluationLoggerModelEvaluator<Scalar> >
00545     loggerThyraModel
00546     = rcp(
00547       new Thyra::DefaultEvaluationLoggerModelEvaluator<Scalar>(
00548         outerModel_,modelEvalLogOut
00549         )
00550       );
00551   outerModel_ = loggerThyraModel; 
00552   
00553   // Manipulating the nominal values
00554   nominalModel_
00555     = rcp(
00556       new Thyra::DefaultNominalBoundsOverrideModelEvaluator<Scalar>(outerModel_,Teuchos::null)
00557       );
00558   outerModel_ = nominalModel_; 
00559 
00560   // Capturing the final point
00561   finalPointModel_
00562     = rcp(
00563       new Thyra::DefaultFinalPointCaptureModelEvaluator<value_type>(outerModel_)
00564       );
00565   outerModel_ = finalPointModel_;
00566 
00567   // Parameter lumping?
00568   if (useParameterLumping_) {
00569     Teuchos::RCP<Thyra::DefaultLumpedParameterModelEvaluator<Scalar> >
00570       lumpedParamModel = Thyra::defaultLumpedParameterModelEvaluator<Scalar>(outerModel_);
00571     lumpedParamModel->setVerbLevel(Teuchos::VERB_LOW);
00572     lumpedParamModel->setParameterList(
00573       sublist(paramList_,LumpedParametersSettings_name));
00574     outerModel_ = lumpedParamModel; 
00575   }
00576   // Note, above we put parameter lumping on the very top so that everything
00577   // like the final point capture and the nominal bounds overrider deal with
00578   // the original parameters, not the lumped parameters.
00579 
00580   //
00581   // Create the NLP
00582   //
00583     
00584   Teuchos::RCP<NLP> nlp;
00585 
00586   switch(solveMode_) {
00587     case SOLVE_MODE_FORWARD: {
00588       RCP<NLPFirstOrderThyraModelEvaluator>
00589         nlpFirstOrder = rcp(
00590           new NLPFirstOrderThyraModelEvaluator(outerModel_,-1,-1)
00591           );
00592       nlpFirstOrder->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00593       nlp = nlpFirstOrder;
00594       break;
00595     }
00596     case SOLVE_MODE_OPTIMIZE: {
00597       // Setup finite difference object
00598       RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > objDirecFiniteDiffCalculator;
00599       if(use_finite_diff_for_obj_) {
00600         objDirecFiniteDiffCalculator = rcp(new Thyra::DirectionalFiniteDiffCalculator<Scalar>());
00601         if(paramList_.get())
00602           objDirecFiniteDiffCalculator->setParameterList(
00603             Teuchos::sublist(paramList_,ObjectiveFiniteDifferenceSettings_name)
00604             );
00605       }
00606       RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > conDirecFiniteDiffCalculator;
00607       if(use_finite_diff_for_con_) {
00608         conDirecFiniteDiffCalculator = rcp(new Thyra::DirectionalFiniteDiffCalculator<Scalar>());
00609         if(paramList_.get())
00610           conDirecFiniteDiffCalculator->setParameterList(
00611             Teuchos::sublist(paramList_,ConstraintsFiniteDifferenceSettings_name)
00612             );
00613       }
00614       if( nonlinearlyElimiateStates_ ) {
00615         // Create a Thyra::NonlinearSolverBase object to solve and eliminate the
00616         // state variables and the state equations
00617         Teuchos::RCP<Thyra::DampenedNewtonNonlinearSolver<Scalar> >
00618           stateSolver = rcp(new Thyra::DampenedNewtonNonlinearSolver<Scalar>()); // ToDo: Replace with MOOCHO!
00619         stateSolver->defaultTol(fwd_newton_tol_);
00620         stateSolver->defaultMaxNewtonIterations(fwd_newton_max_iters_);
00621         stateSolver->useDampenedLineSearch(fwd_newton_dampening_);
00622         stateSolver->maxLineSearchIterations(fwd_newton_max_ls_iters_);
00623         // Create the reduced Thyra::ModelEvaluator object for p -> g_hat(p)
00624         Teuchos::RCP<Thyra::DefaultStateEliminationModelEvaluator<Scalar> >
00625           reducedThyraModel = rcp(new Thyra::DefaultStateEliminationModelEvaluator<Scalar>(outerModel_,stateSolver));
00626         Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
00627           finalReducedThyraModel;
00628         if(use_finite_diff_for_obj_) {
00629           // Create the finite-difference wrapped Thyra::ModelEvaluator object
00630           Teuchos::RCP<Thyra::DefaultFiniteDifferenceModelEvaluator<Scalar> >
00631             fdReducedThyraModel = Thyra::defaultFiniteDifferenceModelEvaluator<Scalar>(
00632               reducedThyraModel, objDirecFiniteDiffCalculator
00633               );
00634           finalReducedThyraModel = fdReducedThyraModel;
00635         }
00636         else {
00637           finalReducedThyraModel = reducedThyraModel;
00638         }
00639         // Wrap the reduced NAND Thyra::ModelEvaluator object in an NLP object
00640         RCP<NLPFirstOrderThyraModelEvaluator>
00641           nlpFirstOrder = rcp(
00642             new NLPFirstOrderThyraModelEvaluator(finalReducedThyraModel,p_idx_,g_idx_)
00643             );
00644         nlpFirstOrder->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00645         nlp = nlpFirstOrder;
00646       }
00647       else {
00648         switch(nlpType_) {
00649           case NLP_TYPE_DIRECT: {
00650             Teuchos::RCP<NLPDirectThyraModelEvaluator>
00651               nlpDirect = rcp(
00652                 new NLPDirectThyraModelEvaluator(
00653                   outerModel_,p_idx_,g_idx_
00654                   ,objDirecFiniteDiffCalculator
00655                   ,conDirecFiniteDiffCalculator
00656                   )
00657                 );
00658             nlpDirect->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00659             nlp = nlpDirect;
00660             break;
00661           }
00662           case NLP_TYPE_FIRST_ORDER: {
00663             RCP<NLPFirstOrderThyraModelEvaluator>
00664               nlpFirstOrder = rcp(
00665                 new NLPFirstOrderThyraModelEvaluator(outerModel_,p_idx_,g_idx_)
00666                 );
00667             nlpFirstOrder->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00668             nlp = nlpFirstOrder;
00669             break;
00670           }
00671           default:
00672             TEST_FOR_EXCEPT(true);
00673         }
00674       }
00675       break;
00676     }
00677     default:
00678       TEST_FOR_EXCEPT(true);
00679   }
00680     
00681   // Set the NLP
00682   solver_.set_nlp(nlp);
00683 
00684 }
00685 
00686 const Teuchos::RCP<Thyra::ModelEvaluator<value_type> >
00687 MoochoThyraSolver::getOrigModel() const
00688 {
00689   return origModel_;
00690 }
00691 
00692 const Teuchos::RCP<Thyra::ModelEvaluator<value_type> >
00693 MoochoThyraSolver::getOuterModel() const
00694 {
00695   return outerModel_;
00696 }
00697 
00698 void MoochoThyraSolver::readInitialGuess(
00699   std::ostream *out_arg
00700   )
00701 {
00702   using Teuchos::OSTab;
00703   using Teuchos::RCP;
00704   using Thyra::clone;
00705   typedef Thyra::ModelEvaluatorBase MEB;
00706 
00707   RCP<Teuchos::FancyOStream>
00708     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00709   RCP<MEB::InArgs<value_type> >
00710     initialGuess = clone(origModel_->getNominalValues()),
00711     lowerBounds = clone(origModel_->getLowerBounds()),
00712     upperBounds = clone(origModel_->getUpperBounds());
00713 
00714   RCP<const Thyra::VectorSpaceBase<value_type> >
00715     x_space = origModel_->get_x_space();
00716   if( 0 != x_space.get() ) {
00717     x_reader_.set_vecSpc(origModel_->get_x_space());
00718     if(this->get_stateVectorIO().get())
00719       x_reader_.set_fileIO(this->get_stateVectorIO());
00720     Teuchos::VerboseObjectTempState<Thyra::ParameterDrivenMultiVectorInput<value_type> >
00721       vots_x_reader(RCP<Thyra::ParameterDrivenMultiVectorInput<value_type> >(&x_reader_,false),out,Teuchos::VERB_LOW);
00722     initialGuess->set_x(
00723       readVectorOverride(
00724         x_reader_,"initial guess for the state \'x\'",initialGuess->get_x()
00725         )
00726       );
00727   }
00728   if( origModel_->Np() > 0 ) {
00729     p_reader_.set_vecSpc(origModel_->get_p_space(p_idx_));
00730     p_l_reader_.set_vecSpc(p_reader_.get_vecSpc());
00731     p_u_reader_.set_vecSpc(p_reader_.get_vecSpc());
00732     if(this->get_parameterVectorIO().get()) {
00733       p_reader_.set_fileIO(this->get_parameterVectorIO());
00734       p_l_reader_.set_fileIO(p_reader_.get_fileIO());
00735       p_u_reader_.set_fileIO(p_reader_.get_fileIO());
00736     }
00737     Teuchos::VerboseObjectTempState<Thyra::ParameterDrivenMultiVectorInput<value_type> >
00738       vots_p_reader(RCP<Thyra::ParameterDrivenMultiVectorInput<value_type> >(&p_reader_,false),out,Teuchos::VERB_LOW);
00739     initialGuess->set_p(
00740       p_idx_,
00741       readVectorOverride(
00742         p_reader_,"initial guess for the parameters \'p\'",
00743         initialGuess->get_p(p_idx_)
00744         )
00745       );
00746     lowerBounds->set_p(
00747       p_idx_,
00748       readVectorOverride(
00749         p_l_reader_,"lower bounds for the parameters \'p\'",lowerBounds->get_p(p_idx_)
00750         )
00751       );
00752     upperBounds->set_p(
00753       p_idx_,
00754       readVectorOverride(
00755         p_u_reader_,"upper bounds for the parameters \'p\'",upperBounds->get_p(p_idx_)
00756         )
00757       );
00758   }
00759   nominalModel_->setNominalValues(initialGuess);
00760   nominalModel_->setLowerBounds(lowerBounds);
00761   nominalModel_->setUpperBounds(upperBounds);
00762 }
00763 
00764 void MoochoThyraSolver::setInitialGuess(
00765   const Teuchos::RCP<const Thyra::ModelEvaluatorBase::InArgs<value_type> > &initialGuess
00766   )
00767 {
00768   nominalModel_->setNominalValues(initialGuess);
00769 }
00770 
00771 void MoochoThyraSolver::setInitialGuess(
00772   const Thyra::ModelEvaluatorBase::InArgs<value_type> &initialGuess
00773   )
00774 {
00775   nominalModel_->setNominalValues(
00776     Teuchos::rcp(new Thyra::ModelEvaluatorBase::InArgs<value_type>(initialGuess))
00777     );
00778 }
00779   
00780 MoochoSolver::ESolutionStatus MoochoThyraSolver::solve()
00781 {
00782   using Teuchos::RCP; using Teuchos::null; using Teuchos::describe;
00783   solver_.update_solver();
00784   std::ostringstream os;
00785   os
00786     << "\n**********************************"
00787     << "\n*** MoochoThyraSolver::solve() ***"
00788     << "\n**********************************\n";
00789   const RCP<const Thyra::VectorSpaceBase<value_type> >
00790     x_space = outerModel_->get_x_space(),
00791     p_space = (
00792       ( p_idx_ >= 0 && outerModel_->Np() > 0 )
00793       ? outerModel_->get_p_space(p_idx_)
00794       : null
00795       );
00796   if( !is_null(x_space) )
00797     os << "\nx_space: " << describe(*x_space,Teuchos::VERB_MEDIUM);
00798   if( !is_null(p_space) )
00799     os << "\np_space: " << describe(*p_space,Teuchos::VERB_MEDIUM);
00800   *solver_.get_console_out() << os.str();
00801   *solver_.get_summary_out() << os.str();
00802   *solver_.get_journal_out() << os.str();
00803   outerModel_->setOStream(Teuchos::getFancyOStream(solver_.get_journal_out()));
00804   return solver_.solve_nlp();
00805 }
00806 
00807 const Thyra::ModelEvaluatorBase::InArgs<value_type>&
00808 MoochoThyraSolver::getFinalPoint() const
00809 {
00810   return finalPointModel_->getFinalPoint();
00811 }
00812 
00813 void MoochoThyraSolver::writeFinalSolution(
00814   std::ostream *out_arg
00815   ) const
00816 {
00817   using Teuchos::OSTab;
00818   Teuchos::RCP<Teuchos::FancyOStream>
00819     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00820   if( stateSoluFileBase_ != "" && finalPointModel_->getFinalPoint().get_x().get() ) {
00821     if(out.get())
00822       *out << "\nWriting the state solution \'x\' to the file(s) with base name \""<<stateSoluFileBase_<<"\" ...\n";
00823     stateVectorIO().writeMultiVectorToFile(
00824       *finalPointModel_->getFinalPoint().get_x(),stateSoluFileBase_
00825       );
00826   }
00827   if(
00828     ( "" != paramSoluFileBase_ )
00829     && ( origModel_->Np() > 0 )
00830     && ( 0 != finalPointModel_->getFinalPoint().get_p(p_idx_).get() )
00831     )
00832   {
00833     if(out.get())
00834       *out << "\nWriting the parameter solution \'p\' to the file(s) with base name \""<<paramSoluFileBase_<<"\" ...\n";
00835     parameterVectorIO().writeMultiVectorToFile(
00836       *finalPointModel_->getFinalPoint().get_p(p_idx_),paramSoluFileBase_
00837       );
00838   }
00839 }
00840 
00841 void MoochoThyraSolver::writeParamsFile(
00842   const std::string &outputXmlFileName
00843   ) const
00844 {
00845   std::string xmlOutputFile
00846     = ( outputXmlFileName.length() ? outputXmlFileName : paramsUsedXmlOutFileName() );
00847   if( paramList_.get() && xmlOutputFile.length() ) {
00848     Teuchos::writeParameterListToXmlFile(*paramList_,xmlOutputFile);
00849   }
00850 }
00851 
00852 } // namespace MoochoPack

Generated on Tue Oct 20 12:51:15 2009 for MOOCHO/Thyra Adapter Software by doxygen 1.4.7