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 OutputOnAllProcesses_name = "Output On All Processes";
00146 const bool OutputOnAllProcesses_default = false;
00147 
00148 const std::string ShowModelEvaluatorTrace_name = "Show Model Evaluator Trace";
00149 const bool ShowModelEvaluatorTrace_default = "false";
00150 
00151 const std::string StateGuess_name = "State Guess";
00152 
00153 const std::string ParamGuess_name = "Parameter Guess";
00154 
00155 const std::string ParamLowerBounds_name = "Parameter Lower Bounds";
00156 
00157 const std::string ParamUpperBounds_name = "Parameter Upper Bounds";
00158 
00159 const std::string StateSoluFileBaseName_name = "State Solution File Base Name";
00160 const std::string StateSoluFileBaseName_default = "";
00161 
00162 const std::string ParamSoluFileBaseName_name = "Parameters Solution File Base Name";
00163 const std::string ParamSoluFileBaseName_default = "";
00164 
00165 } // namespace
00166 
00167 namespace MoochoPack {
00168 
00169 // Constructors/initialization
00170 
00171 MoochoThyraSolver::MoochoThyraSolver(
00172   const std::string    &paramsXmlFileName
00173   ,const std::string   &extraParamsXmlString
00174   ,const std::string   &paramsUsedXmlOutFileName
00175   ,const std::string   &paramsXmlFileNameOption
00176   ,const std::string   &extraParamsXmlStringOption
00177   ,const std::string   &paramsUsedXmlOutFileNameOption
00178   )
00179   :paramsXmlFileName_(paramsXmlFileName)
00180   ,extraParamsXmlString_(extraParamsXmlString)
00181   ,paramsUsedXmlOutFileName_(paramsUsedXmlOutFileName)
00182   ,paramsXmlFileNameOption_(paramsXmlFileNameOption)
00183   ,extraParamsXmlStringOption_(extraParamsXmlStringOption)
00184   ,paramsUsedXmlOutFileNameOption_(paramsUsedXmlOutFileNameOption)
00185   ,stateVectorIO_(Teuchos::rcp(new Thyra::DefaultSpmdMultiVectorFileIO<value_type>))
00186   ,parameterVectorIO_(Teuchos::rcp(new Thyra::DefaultSpmdMultiVectorFileIO<value_type>))
00187   ,solveMode_(SOLVE_MODE_OPTIMIZE)
00188   ,nlpType_(NLP_TYPE_FIRST_ORDER)
00189   ,nonlinearlyElimiateStates_(false)
00190   ,use_finite_diff_for_obj_(false)
00191   ,use_finite_diff_for_con_(false)
00192   ,fwd_newton_tol_(-1.0)
00193   ,fwd_newton_max_iters_(20)
00194   ,fwd_newton_dampening_(false)
00195   ,fwd_newton_max_ls_iters_(20)
00196   ,useInvObjFunc_(false)
00197   ,useParameterLumping_(false)
00198   ,outputFileTag_("")
00199   ,showModelEvaluatorTrace_(false)
00200   ,stateSoluFileBase_("")
00201   ,paramSoluFileBase_("")
00202 {}
00203 
00204 MoochoThyraSolver::~MoochoThyraSolver()
00205 {}
00206 
00207 void MoochoThyraSolver::setupCLP(
00208   Teuchos::CommandLineProcessor *clp
00209   )
00210 {
00211   TEST_FOR_EXCEPT(0==clp);
00212   solver_.setup_commandline_processor(clp);
00213   clp->setOption(
00214     paramsXmlFileNameOption().c_str(),&paramsXmlFileName_
00215     ,"Name of an XML file containing parameters for linear solver options to be appended first."
00216     );
00217   clp->setOption(
00218     extraParamsXmlStringOption().c_str(),&extraParamsXmlString_
00219     ,"An XML string containing linear solver parameters to be appended second."
00220     );
00221   clp->setOption(
00222     paramsUsedXmlOutFileNameOption().c_str(),&paramsUsedXmlOutFileName_
00223     ,"Name of an XML file that can be written with the parameter list after it has been used on completion of this program."
00224     );
00225 }
00226 
00227 void MoochoThyraSolver::readParameters( std::ostream *out_arg )
00228 {
00229   Teuchos::RCP<Teuchos::FancyOStream>
00230     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00231   Teuchos::OSTab tab(out);
00232   if(out.get()) *out << "\nMoochoThyraSolver::readParameters(...):\n";
00233   Teuchos::OSTab tab2(out);
00234   Teuchos::RCP<Teuchos::ParameterList>
00235     paramList = this->getNonconstParameterList();
00236   if(!paramList.get()) {
00237     if(out.get()) *out << "\nCreating a new Teuchos::ParameterList ...\n";
00238     paramList = Teuchos::rcp(new Teuchos::ParameterList("MoochoThyraSolver"));
00239   }
00240   if(paramsXmlFileName().length()) {
00241     if(out.get()) *out << "\nReading parameters from XML file \""<<paramsXmlFileName()<<"\" ...\n";
00242     Teuchos::updateParametersFromXmlFile(paramsXmlFileName(),&*paramList);
00243   }
00244   if(extraParamsXmlString().length()) {
00245     if(out.get())
00246       *out << "\nAppending extra parameters from the XML string \""<<extraParamsXmlString()<<"\" ...\n";
00247     Teuchos::updateParametersFromXmlString(extraParamsXmlString(),&*paramList);
00248   }
00249   if( paramsXmlFileName().length() || extraParamsXmlString().length() ) {
00250     typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
00251     if(out.get()) {
00252       *out  << "\nUpdated parameter list:\n";
00253       paramList->print(
00254         *out,PLPrintOptions().indent(2).showTypes(true)
00255         );
00256     }
00257     this->setParameterList(paramList);
00258   }
00259 }
00260 
00261 // Overridden from ParameterListAcceptor
00262 
00263 void MoochoThyraSolver::setParameterList(
00264   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00265   )
00266 {
00267   typedef MoochoPack::MoochoSolver MS;
00268   TEST_FOR_EXCEPT(!paramList.get());
00269   paramList->validateParameters(*getValidParameters(),0); // Just validate my level!
00270   paramList_ = paramList;
00271   solveMode_ = solveModeValidator->getIntegralValue(
00272     *paramList_,SolveMode_name,SolveMode_default);
00273   nlpType_ = nlpTypeValidator->getIntegralValue(
00274     *paramList_,NLPType_name,NLPType_default);
00275   nonlinearlyElimiateStates_ = paramList_->get(
00276     NonlinearlyEliminateStates_name,NonlinearlyEliminateStates_default);
00277   use_finite_diff_for_obj_ = paramList_->get(
00278     UseFiniteDifferencesForObjective_name,UseFiniteDifferencesForObjective_default);
00279   use_finite_diff_for_con_ = paramList_->get(
00280     UseFiniteDifferencesForConstraints_name,UseFiniteDifferencesForConstraints_default);
00281   fwd_newton_tol_ = paramList_->get(
00282     FwdNewtonTol_name,FwdNewtonTol_default);
00283   fwd_newton_max_iters_ = paramList_->get(
00284     FwdNewtonMaxIters_name,FwdNewtonMaxIters_default);
00285   fwd_newton_dampening_ = paramList_->get(
00286     ForwardNewtonDampening_name,ForwardNewtonDampening_default);
00287   fwd_newton_max_ls_iters_ = paramList_->get(
00288     FwdNewtonMaxLineSearchIters_name,FwdNewtonMaxLineSearchIters_default);
00289   useInvObjFunc_ = paramList_->get(
00290     UseBuiltInInverseObjectiveFunction_name,UseBuiltInInverseObjectiveFunction_default);
00291   useParameterLumping_ = paramList_->get(
00292     UseParameterLumping_name, UseParameterLumping_default);
00293   outputFileTag_ = paramList->get(
00294     OutputFileTag_name,OutputFileTag_default);
00295   solver_.set_output_file_tag(outputFileTag_);
00296   solver_.set_output_context("",
00297     paramList_->get(OutputOnAllProcesses_name, OutputOnAllProcesses_default)
00298     ? MS::OUTPUT_TO_BLACK_HOLE_FALSE
00299     : MS::OUTPUT_TO_BLACK_HOLE_DEFAULT
00300     );
00301   showModelEvaluatorTrace_ = paramList->get(
00302     ShowModelEvaluatorTrace_name,ShowModelEvaluatorTrace_default);
00303   x_reader_.setParameterList(
00304     sublist(paramList_,StateGuess_name)
00305     );
00306   p_reader_.setParameterList(
00307     sublist(paramList_,ParamGuess_name)
00308     );
00309   p_l_reader_.setParameterList(
00310     sublist(paramList_,ParamLowerBounds_name)
00311     );
00312   p_u_reader_.setParameterList(
00313     sublist(paramList_,ParamUpperBounds_name)
00314     );
00315   stateSoluFileBase_ = paramList_->get(
00316     StateSoluFileBaseName_name,StateSoluFileBaseName_default);
00317   paramSoluFileBase_ = paramList_->get(
00318     ParamSoluFileBaseName_name,ParamSoluFileBaseName_default);
00319 #ifdef TEUCHOS_DEBUG
00320   paramList->validateParameters(*getValidParameters(),0); // Just validate my level!
00321 #endif
00322 }
00323 
00324 Teuchos::RCP<Teuchos::ParameterList>
00325 MoochoThyraSolver::getNonconstParameterList()
00326 {
00327   return paramList_;
00328 }
00329 
00330 Teuchos::RCP<Teuchos::ParameterList>
00331 MoochoThyraSolver::unsetParameterList()
00332 {
00333   Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
00334   paramList_ = Teuchos::null;
00335   return _paramList;
00336 }
00337 
00338 Teuchos::RCP<const Teuchos::ParameterList>
00339 MoochoThyraSolver::getParameterList() const
00340 {
00341   return paramList_;
00342 }
00343 
00344 Teuchos::RCP<const Teuchos::ParameterList>
00345 MoochoThyraSolver::getValidParameters() const
00346 {
00347   static Teuchos::RCP<Teuchos::ParameterList> pl;
00348   if(pl.get()==NULL) {
00349     pl = Teuchos::rcp(new Teuchos::ParameterList());
00350     pl->set(
00351       SolveMode_name,SolveMode_default
00352       ,"The type of solve to perform."
00353       ,solveModeValidator
00354       );
00355     pl->set(
00356       NLPType_name,NLPType_default
00357       ,"The type of MOOCHO NLP subclass to use."
00358       ,nlpTypeValidator
00359       );
00360     pl->set(
00361       NonlinearlyEliminateStates_name,NonlinearlyEliminateStates_default
00362       ,"If true, then the model's state equations and state variables\n"
00363       "are nonlinearlly eliminated using a forward solver."
00364       );
00365     pl->set(
00366       UseFiniteDifferencesForObjective_name,UseFiniteDifferencesForObjective_default
00367       ,"Use finite differences for missing objective function derivatives (Direct NLP only).\n"
00368       "See the options in the sublist \"" + ObjectiveFiniteDifferenceSettings_name + "\"."
00369       );
00370     {
00371       Thyra::DirectionalFiniteDiffCalculator<Scalar> dfdcalc;
00372       {
00373         Teuchos::ParameterList
00374           &fdSublist = pl->sublist(ObjectiveFiniteDifferenceSettings_name);
00375         fdSublist.setParameters(*dfdcalc.getValidParameters());
00376       }
00377       pl->set(
00378         UseFiniteDifferencesForConstraints_name,UseFiniteDifferencesForConstraints_default
00379         ,"Use  finite differences for missing constraint derivatives (Direct NLP only).\n"
00380         "See the   options in the sublist \"" + ConstraintsFiniteDifferenceSettings_name + "\"."
00381         );
00382       {
00383         Teuchos::ParameterList
00384           &fdSublist = pl->sublist(ConstraintsFiniteDifferenceSettings_name);
00385         fdSublist.setParameters(*dfdcalc.getValidParameters());
00386       }
00387     }
00388     pl->set(
00389       FwdNewtonTol_name,FwdNewtonTol_default
00390       ,"Tolarance used for the forward state solver in eliminating\n"
00391       "the state equations/variables."
00392       );
00393     pl->set(
00394       FwdNewtonMaxIters_name,FwdNewtonMaxIters_default
00395       ,"Maximum number of iterations allowed for the forward state\n"
00396       "solver in eliminating the state equations/variables."
00397       );
00398     pl->set(
00399       ForwardNewtonDampening_name, ForwardNewtonDampening_default,
00400       "If true, then the state elimination nonlinear solver will\n"
00401       "use a dampened line search.  Otherwise, it will just take fulls steps."
00402       );
00403     pl->set(
00404       FwdNewtonMaxLineSearchIters_name, FwdNewtonMaxLineSearchIters_default,
00405       "Maximum number of linea search iterations per newton iteration\n"
00406       "allowed for the forward state solver in eliminating the state equations/variables."
00407       );
00408     pl->set(
00409       UseBuiltInInverseObjectiveFunction_name,UseBuiltInInverseObjectiveFunction_default
00410       ,"Use a built-in form of a simple inverse objection function instead\n"
00411       "of a a response function contained in the underlying model evaluator\n"
00412       "object itself.  The settings are contained in the sublist\n"
00413       "\""+InverseObjectiveFunctionSettings_name+"\".\n"
00414       "Note that this feature allows the client to form a useful type\n"
00415       "of optimization problem just with a model that supports only the\n"
00416       "parameterized state function f(x,p)=0."
00417       );
00418     {
00419       Teuchos::RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00420         inverseModel = rcp(new Thyra::DefaultInverseModelEvaluator<Scalar>());
00421       pl->sublist(
00422         InverseObjectiveFunctionSettings_name,false
00423         ,"Settings for the built-in inverse objective function.\n"
00424         "See the outer parameter \""+UseBuiltInInverseObjectiveFunction_name+"\"."
00425         ).setParameters(*inverseModel->getValidParameters());
00426     }
00427     pl->set(
00428       UseParameterLumping_name,UseParameterLumping_default,
00429       "Use a reduced basis to lump optimization parameters as\n"
00430       "p_orig = P_basis * p.  If set to true, then the settings\n"
00431       "in \""+LumpedParametersSettings_name+"\" determine how the\n"
00432       "parameter basis is set.  This feature can be used to safely\n"
00433       "regularize a problem if there are linearly dependent parameters\n"
00434       "and will generally speed up the optimiztation algorithms."
00435       );
00436     {
00437       Teuchos::RCP<Thyra::DefaultLumpedParameterModelEvaluator<Scalar> >
00438         lumpedParamModel = rcp(new Thyra::DefaultLumpedParameterModelEvaluator<Scalar>);
00439       pl->sublist(
00440         LumpedParametersSettings_name,false
00441         ,"Settings for parameter lumping.\n"
00442         "See the outer parameter \""+UseParameterLumping_name+"\"."
00443         ).setParameters(*lumpedParamModel->getValidParameters());
00444     }
00445     pl->set(OutputFileTag_name,OutputFileTag_default,
00446       "A tag that is attached to every output file that is created by the\n"
00447       "solver.  If empty \"\", then no tag is used.  This option simply is\n"
00448       "passed into the set_output_context(...) function on the underlying\n"
00449       "MoochoPack::MoochoSolver object.  Therefore, this same parameter\n"
00450       "can be set in code as well without going through the parameter list.");
00451     pl->set(OutputOnAllProcesses_name, OutputOnAllProcesses_default,
00452       "If set to true, then the console output and all MOOCHO files will be\n"
00453       "written for every process.  This helps in debugging very hard\n"
00454       "parallel problems.");
00455     pl->set(ShowModelEvaluatorTrace_name,ShowModelEvaluatorTrace_default
00456       ,"Determine if a trace of the objective function will be shown or not\n"
00457       "when the NLP is evaluated."
00458       );
00459     if(this->get_stateVectorIO().get())
00460       x_reader_.set_fileIO(this->get_stateVectorIO());
00461     pl->sublist(StateGuess_name).setParameters(*x_reader_.getValidParameters());
00462     if(this->get_parameterVectorIO().get()) {
00463       p_reader_.set_fileIO(this->get_parameterVectorIO());
00464       p_l_reader_.set_fileIO(this->get_parameterVectorIO());
00465       p_u_reader_.set_fileIO(this->get_parameterVectorIO());
00466       pl->sublist(ParamGuess_name).setParameters(*p_reader_.getValidParameters());
00467       pl->sublist(ParamLowerBounds_name).setParameters(*p_l_reader_.getValidParameters());
00468       pl->sublist(ParamUpperBounds_name).setParameters(*p_u_reader_.getValidParameters());
00469     }
00470     pl->set(
00471       StateSoluFileBaseName_name,StateSoluFileBaseName_default
00472       ,"If specified, a file with this basename will be written to with\n"
00473       "the final value of the state variables.  A different file for each\n"
00474       "process will be created.  Note that these files can be used for the\n"
00475       "initial guess for the state variables."
00476       );
00477     pl->set(
00478       ParamSoluFileBaseName_name,ParamSoluFileBaseName_default
00479       ,"If specified, a file with this basename will be written to with\n"
00480       "the final value of the parameters.  A different file for each\n"
00481       "process will be created.  Note that these files can be used for the\n"
00482       "initial guess for the parameters."
00483       );
00484   }
00485   return pl;
00486 }
00487 
00488 // Misc Access/Setup
00489 
00490 void MoochoThyraSolver::setSolveMode( const ESolveMode solveMode )
00491 {
00492   solveMode_ = solveMode;
00493 }
00494 
00495 MoochoThyraSolver::ESolveMode
00496 MoochoThyraSolver::getSolveMode() const
00497 {
00498   return solveMode_;
00499 }
00500 
00501 MoochoSolver& MoochoThyraSolver::getSolver()
00502 {
00503   return solver_;
00504 }
00505 
00506 const MoochoSolver& MoochoThyraSolver::getSolver() const
00507 {
00508   return solver_;
00509 }
00510 
00511 // Model specification, setup, solve, and solution extraction.
00512 
00513 void MoochoThyraSolver::setModel(
00514   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> > &origModel,
00515   const int p_idx,
00516   const int g_idx
00517   )
00518 {
00519 
00520   using Teuchos::rcp;
00521   using Teuchos::sublist;
00522   using NLPInterfacePack::NLP;
00523   using NLPInterfacePack::NLPDirectThyraModelEvaluator;
00524   using NLPInterfacePack::NLPFirstOrderThyraModelEvaluator;
00525 
00526   origModel_ = origModel;
00527   p_idx_ = p_idx;
00528   g_idx_ = g_idx;
00529 
00530   const int procRank = Teuchos::GlobalMPISession::getRank();
00531   //const int numProcs = Teuchos::GlobalMPISession::getNProc();
00532 
00533   //
00534   // Wrap the orginal model in different decorators
00535   //
00536 
00537   outerModel_ = origModel_;
00538 
00539   // Inverse objective (and parameter regularization)?
00540   if (useInvObjFunc_) {
00541     Teuchos::RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00542       inverseModel = Thyra::defaultInverseModelEvaluator<Scalar>(outerModel_);
00543     inverseModel->setVerbLevel(Teuchos::VERB_LOW);
00544     inverseModel->set_observationTargetIO(get_stateVectorIO());
00545     inverseModel->set_parameterBaseIO(get_parameterVectorIO());
00546     inverseModel->setParameterList(
00547       Teuchos::sublist(paramList_,InverseObjectiveFunctionSettings_name) );
00548     outerModel_ = inverseModel; 
00549     g_idx_ = inverseModel->Ng()-1;
00550   }
00551   
00552   // Evaluation loging
00553   Teuchos::RCP<std::ostream>
00554     modelEvalLogOut = Teuchos::fancyOStream(
00555       solver_.generate_output_file("ModelEvaluationLog")
00556       );
00557   Teuchos::RCP<Thyra::DefaultEvaluationLoggerModelEvaluator<Scalar> >
00558     loggerThyraModel
00559     = rcp(
00560       new Thyra::DefaultEvaluationLoggerModelEvaluator<Scalar>(
00561         outerModel_,modelEvalLogOut
00562         )
00563       );
00564   outerModel_ = loggerThyraModel; 
00565   
00566   // Manipulating the nominal values
00567   nominalModel_
00568     = rcp(
00569       new Thyra::DefaultNominalBoundsOverrideModelEvaluator<Scalar>(outerModel_,Teuchos::null)
00570       );
00571   outerModel_ = nominalModel_; 
00572 
00573   // Capturing the final point
00574   finalPointModel_
00575     = rcp(
00576       new Thyra::DefaultFinalPointCaptureModelEvaluator<value_type>(outerModel_)
00577       );
00578   outerModel_ = finalPointModel_;
00579 
00580   // Parameter lumping?
00581   if (useParameterLumping_) {
00582     Teuchos::RCP<Thyra::DefaultLumpedParameterModelEvaluator<Scalar> >
00583       lumpedParamModel = Thyra::defaultLumpedParameterModelEvaluator<Scalar>(outerModel_);
00584     lumpedParamModel->setVerbLevel(Teuchos::VERB_LOW);
00585     lumpedParamModel->setParameterList(
00586       sublist(paramList_,LumpedParametersSettings_name));
00587     outerModel_ = lumpedParamModel; 
00588   }
00589   // Note, above we put parameter lumping on the very top so that everything
00590   // like the final point capture and the nominal bounds overrider deal with
00591   // the original parameters, not the lumped parameters.
00592 
00593   //
00594   // Create the NLP
00595   //
00596     
00597   Teuchos::RCP<NLP> nlp;
00598 
00599   switch(solveMode_) {
00600     case SOLVE_MODE_FORWARD: {
00601       RCP<NLPFirstOrderThyraModelEvaluator>
00602         nlpFirstOrder = rcp(
00603           new NLPFirstOrderThyraModelEvaluator(outerModel_,-1,-1)
00604           );
00605       nlpFirstOrder->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00606       nlp = nlpFirstOrder;
00607       break;
00608     }
00609     case SOLVE_MODE_OPTIMIZE: {
00610       // Setup finite difference object
00611       RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > objDirecFiniteDiffCalculator;
00612       if(use_finite_diff_for_obj_) {
00613         objDirecFiniteDiffCalculator = rcp(new Thyra::DirectionalFiniteDiffCalculator<Scalar>());
00614         if(paramList_.get())
00615           objDirecFiniteDiffCalculator->setParameterList(
00616             Teuchos::sublist(paramList_,ObjectiveFiniteDifferenceSettings_name)
00617             );
00618       }
00619       RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > conDirecFiniteDiffCalculator;
00620       if(use_finite_diff_for_con_) {
00621         conDirecFiniteDiffCalculator = rcp(new Thyra::DirectionalFiniteDiffCalculator<Scalar>());
00622         if(paramList_.get())
00623           conDirecFiniteDiffCalculator->setParameterList(
00624             Teuchos::sublist(paramList_,ConstraintsFiniteDifferenceSettings_name)
00625             );
00626       }
00627       if( nonlinearlyElimiateStates_ ) {
00628         // Create a Thyra::NonlinearSolverBase object to solve and eliminate the
00629         // state variables and the state equations
00630         Teuchos::RCP<Thyra::DampenedNewtonNonlinearSolver<Scalar> >
00631           stateSolver = rcp(new Thyra::DampenedNewtonNonlinearSolver<Scalar>()); // ToDo: Replace with MOOCHO!
00632         stateSolver->defaultTol(fwd_newton_tol_);
00633         stateSolver->defaultMaxNewtonIterations(fwd_newton_max_iters_);
00634         stateSolver->useDampenedLineSearch(fwd_newton_dampening_);
00635         stateSolver->maxLineSearchIterations(fwd_newton_max_ls_iters_);
00636         // Create the reduced Thyra::ModelEvaluator object for p -> g_hat(p)
00637         Teuchos::RCP<Thyra::DefaultStateEliminationModelEvaluator<Scalar> >
00638           reducedThyraModel = rcp(new Thyra::DefaultStateEliminationModelEvaluator<Scalar>(outerModel_,stateSolver));
00639         Teuchos::RCP<Thyra::ModelEvaluator<Scalar> >
00640           finalReducedThyraModel;
00641         if(use_finite_diff_for_obj_) {
00642           // Create the finite-difference wrapped Thyra::ModelEvaluator object
00643           Teuchos::RCP<Thyra::DefaultFiniteDifferenceModelEvaluator<Scalar> >
00644             fdReducedThyraModel = Thyra::defaultFiniteDifferenceModelEvaluator<Scalar>(
00645               reducedThyraModel, objDirecFiniteDiffCalculator
00646               );
00647           finalReducedThyraModel = fdReducedThyraModel;
00648         }
00649         else {
00650           finalReducedThyraModel = reducedThyraModel;
00651         }
00652         // Wrap the reduced NAND Thyra::ModelEvaluator object in an NLP object
00653         RCP<NLPFirstOrderThyraModelEvaluator>
00654           nlpFirstOrder = rcp(
00655             new NLPFirstOrderThyraModelEvaluator(finalReducedThyraModel,p_idx_,g_idx_)
00656             );
00657         nlpFirstOrder->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00658         nlp = nlpFirstOrder;
00659       }
00660       else {
00661         switch(nlpType_) {
00662           case NLP_TYPE_DIRECT: {
00663             Teuchos::RCP<NLPDirectThyraModelEvaluator>
00664               nlpDirect = rcp(
00665                 new NLPDirectThyraModelEvaluator(
00666                   outerModel_,p_idx_,g_idx_
00667                   ,objDirecFiniteDiffCalculator
00668                   ,conDirecFiniteDiffCalculator
00669                   )
00670                 );
00671             nlpDirect->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00672             nlp = nlpDirect;
00673             break;
00674           }
00675           case NLP_TYPE_FIRST_ORDER: {
00676             RCP<NLPFirstOrderThyraModelEvaluator>
00677               nlpFirstOrder = rcp(
00678                 new NLPFirstOrderThyraModelEvaluator(outerModel_,p_idx_,g_idx_)
00679                 );
00680             nlpFirstOrder->showModelEvaluatorTrace(showModelEvaluatorTrace_);
00681             nlp = nlpFirstOrder;
00682             break;
00683           }
00684           default:
00685             TEST_FOR_EXCEPT(true);
00686         }
00687       }
00688       break;
00689     }
00690     default:
00691       TEST_FOR_EXCEPT(true);
00692   }
00693     
00694   // Set the NLP
00695   solver_.set_nlp(nlp);
00696 
00697 }
00698 
00699 const Teuchos::RCP<Thyra::ModelEvaluator<value_type> >
00700 MoochoThyraSolver::getOrigModel() const
00701 {
00702   return origModel_;
00703 }
00704 
00705 const Teuchos::RCP<Thyra::ModelEvaluator<value_type> >
00706 MoochoThyraSolver::getOuterModel() const
00707 {
00708   return outerModel_;
00709 }
00710 
00711 void MoochoThyraSolver::readInitialGuess(
00712   std::ostream *out_arg
00713   )
00714 {
00715   using Teuchos::OSTab;
00716   using Teuchos::RCP;
00717   using Thyra::clone;
00718   typedef Thyra::ModelEvaluatorBase MEB;
00719 
00720   RCP<Teuchos::FancyOStream>
00721     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00722   RCP<MEB::InArgs<value_type> >
00723     initialGuess = clone(origModel_->getNominalValues()),
00724     lowerBounds = clone(origModel_->getLowerBounds()),
00725     upperBounds = clone(origModel_->getUpperBounds());
00726 
00727   RCP<const Thyra::VectorSpaceBase<value_type> >
00728     x_space = origModel_->get_x_space();
00729   if( 0 != x_space.get() ) {
00730     x_reader_.set_vecSpc(origModel_->get_x_space());
00731     if(this->get_stateVectorIO().get())
00732       x_reader_.set_fileIO(this->get_stateVectorIO());
00733     Teuchos::VerboseObjectTempState<Thyra::ParameterDrivenMultiVectorInput<value_type> >
00734       vots_x_reader(RCP<Thyra::ParameterDrivenMultiVectorInput<value_type> >(&x_reader_,false),out,Teuchos::VERB_LOW);
00735     initialGuess->set_x(
00736       readVectorOverride(
00737         x_reader_,"initial guess for the state \'x\'",initialGuess->get_x()
00738         )
00739       );
00740   }
00741   if( origModel_->Np() > 0 ) {
00742     p_reader_.set_vecSpc(origModel_->get_p_space(p_idx_));
00743     p_l_reader_.set_vecSpc(p_reader_.get_vecSpc());
00744     p_u_reader_.set_vecSpc(p_reader_.get_vecSpc());
00745     if(this->get_parameterVectorIO().get()) {
00746       p_reader_.set_fileIO(this->get_parameterVectorIO());
00747       p_l_reader_.set_fileIO(p_reader_.get_fileIO());
00748       p_u_reader_.set_fileIO(p_reader_.get_fileIO());
00749     }
00750     Teuchos::VerboseObjectTempState<Thyra::ParameterDrivenMultiVectorInput<value_type> >
00751       vots_p_reader(RCP<Thyra::ParameterDrivenMultiVectorInput<value_type> >(&p_reader_,false),out,Teuchos::VERB_LOW);
00752     initialGuess->set_p(
00753       p_idx_,
00754       readVectorOverride(
00755         p_reader_,"initial guess for the parameters \'p\'",
00756         initialGuess->get_p(p_idx_)
00757         )
00758       );
00759     lowerBounds->set_p(
00760       p_idx_,
00761       readVectorOverride(
00762         p_l_reader_,"lower bounds for the parameters \'p\'",lowerBounds->get_p(p_idx_)
00763         )
00764       );
00765     upperBounds->set_p(
00766       p_idx_,
00767       readVectorOverride(
00768         p_u_reader_,"upper bounds for the parameters \'p\'",upperBounds->get_p(p_idx_)
00769         )
00770       );
00771   }
00772   nominalModel_->setNominalValues(initialGuess);
00773   nominalModel_->setLowerBounds(lowerBounds);
00774   nominalModel_->setUpperBounds(upperBounds);
00775 }
00776 
00777 void MoochoThyraSolver::setInitialGuess(
00778   const Teuchos::RCP<const Thyra::ModelEvaluatorBase::InArgs<value_type> > &initialGuess
00779   )
00780 {
00781   nominalModel_->setNominalValues(initialGuess);
00782 }
00783 
00784 void MoochoThyraSolver::setInitialGuess(
00785   const Thyra::ModelEvaluatorBase::InArgs<value_type> &initialGuess
00786   )
00787 {
00788   nominalModel_->setNominalValues(
00789     Teuchos::rcp(new Thyra::ModelEvaluatorBase::InArgs<value_type>(initialGuess))
00790     );
00791 }
00792   
00793 MoochoSolver::ESolutionStatus MoochoThyraSolver::solve()
00794 {
00795   using Teuchos::RCP; using Teuchos::null; using Teuchos::describe;
00796   solver_.update_solver();
00797   std::ostringstream os;
00798   os
00799     << "\n**********************************"
00800     << "\n*** MoochoThyraSolver::solve() ***"
00801     << "\n**********************************\n";
00802   const RCP<const Thyra::VectorSpaceBase<value_type> >
00803     x_space = outerModel_->get_x_space(),
00804     p_space = (
00805       ( p_idx_ >= 0 && outerModel_->Np() > 0 )
00806       ? outerModel_->get_p_space(p_idx_)
00807       : null
00808       );
00809   if( !is_null(x_space) )
00810     os << "\nx_space: " << describe(*x_space,Teuchos::VERB_MEDIUM);
00811   if( !is_null(p_space) )
00812     os << "\np_space: " << describe(*p_space,Teuchos::VERB_MEDIUM);
00813   *solver_.get_console_out() << os.str();
00814   *solver_.get_summary_out() << os.str();
00815   *solver_.get_journal_out() << os.str();
00816   outerModel_->setOStream(Teuchos::getFancyOStream(solver_.get_journal_out()));
00817   return solver_.solve_nlp();
00818 }
00819 
00820 const Thyra::ModelEvaluatorBase::InArgs<value_type>&
00821 MoochoThyraSolver::getFinalPoint() const
00822 {
00823   return finalPointModel_->getFinalPoint();
00824 }
00825 
00826 void MoochoThyraSolver::writeFinalSolution(
00827   std::ostream *out_arg
00828   ) const
00829 {
00830   using Teuchos::OSTab;
00831   Teuchos::RCP<Teuchos::FancyOStream>
00832     out = Teuchos::getFancyOStream(Teuchos::rcp(out_arg,false));
00833   if( stateSoluFileBase_ != "" && finalPointModel_->getFinalPoint().get_x().get() ) {
00834     if(out.get())
00835       *out << "\nWriting the state solution \'x\' to the file(s) with base name \""<<stateSoluFileBase_<<"\" ...\n";
00836     stateVectorIO().writeMultiVectorToFile(
00837       *finalPointModel_->getFinalPoint().get_x(),stateSoluFileBase_
00838       );
00839   }
00840   if(
00841     ( "" != paramSoluFileBase_ )
00842     && ( origModel_->Np() > 0 )
00843     && ( 0 != finalPointModel_->getFinalPoint().get_p(p_idx_).get() )
00844     )
00845   {
00846     if(out.get())
00847       *out << "\nWriting the parameter solution \'p\' to the file(s) with base name \""<<paramSoluFileBase_<<"\" ...\n";
00848     parameterVectorIO().writeMultiVectorToFile(
00849       *finalPointModel_->getFinalPoint().get_p(p_idx_),paramSoluFileBase_
00850       );
00851   }
00852 }
00853 
00854 void MoochoThyraSolver::writeParamsFile(
00855   const std::string &outputXmlFileName
00856   ) const
00857 {
00858   std::string xmlOutputFile
00859     = ( outputXmlFileName.length() ? outputXmlFileName : paramsUsedXmlOutFileName() );
00860   if( paramList_.get() && xmlOutputFile.length() ) {
00861     Teuchos::writeParameterListToXmlFile(*paramList_,xmlOutputFile);
00862   }
00863 }
00864 
00865 } // namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:09:53 2011 for MOOCHO/Thyra Adapter Software by  doxygen 1.6.3