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