newton.hpp
Go to the documentation of this file.
1 #ifndef STAN_SERVICES_OPTIMIZE_NEWTON_HPP
2 #define STAN_SERVICES_OPTIMIZE_NEWTON_HPP
3 
14 #include <cmath>
15 #include <limits>
16 #include <string>
17 #include <vector>
18 
19 namespace stan {
20  namespace services {
21  namespace optimize {
22 
23  /**
24  * Runs the Newton algorithm for a model.
25  *
26  * @tparam Model A model implementation
27  * @param[in] model the Stan model instantiated with data
28  * @param[in] init var context for initialization
29  * @param[in] random_seed random seed for the random number generator
30  * @param[in] chain chain id to advance the pseudo random number generator
31  * @param[in] init_radius radius to initialize
32  * @param[in] num_iterations maximum number of iterations
33  * @param[in] save_iterations indicates whether all the interations should
34  * be saved
35  * @param[in,out] interrupt callback to be called every iteration
36  * @param[in,out] logger Logger for messages
37  * @param[in,out] init_writer Writer callback for unconstrained inits
38  * @param[in,out] parameter_writer output for parameter values
39  * @return error_codes::OK if successful
40  */
41  template <class Model>
43  unsigned int random_seed, unsigned int chain,
44  double init_radius, int num_iterations,
45  bool save_iterations,
46  callbacks::interrupt& interrupt,
47  callbacks::logger& logger,
48  callbacks::writer& init_writer,
49  callbacks::writer& parameter_writer) {
50  boost::ecuyer1988 rng = util::create_rng(random_seed, chain);
51 
52  std::vector<int> disc_vector;
53  std::vector<double> cont_vector
54  = util::initialize<false>(model, init, rng, init_radius, false,
55  logger, init_writer);
56 
57 
58  double lp(0);
59  try {
60  std::stringstream message;
61  lp = model.template log_prob<false, false>(cont_vector, disc_vector,
62  &message);
63  logger.info(message);
64  } catch (const std::exception& e) {
65  logger.info("");
66  logger.info("Informational Message: The current Metropolis"
67  " proposal is about to be rejected because of"
68  " the following issue:");
69  logger.info(e.what());
70  logger.info("If this warning occurs sporadically, such as"
71  " for highly constrained variable types like"
72  " covariance matrices, then the sampler is fine,");
73  logger.info("but if this warning occurs often then your model"
74  " may be either severely ill-conditioned or"
75  " misspecified.");
76  lp = -std::numeric_limits<double>::infinity();
77  }
78 
79  std::stringstream msg;
80  msg << "Initial log joint probability = " << lp;
81  logger.info(msg);
82 
83  std::vector<std::string> names;
84  names.push_back("lp__");
85  model.constrained_param_names(names, true, true);
86  parameter_writer(names);
87 
88  double lastlp = lp;
89  for (int m = 0; m < num_iterations; m++) {
90  if (save_iterations) {
91  std::vector<double> values;
92  std::stringstream ss;
93  model.write_array(rng, cont_vector, disc_vector, values,
94  true, true, &ss);
95  if (ss.str().length() > 0)
96  logger.info(ss);
97  values.insert(values.begin(), lp);
98  parameter_writer(values);
99  }
100  interrupt();
101  lastlp = lp;
102  lp = stan::optimization::newton_step(model, cont_vector, disc_vector);
103 
104  std::stringstream msg2;
105  msg2 << "Iteration "
106  << std::setw(2) << (m + 1) << "."
107  << " Log joint probability = " << std::setw(10) << lp
108  << ". Improved by " << (lp - lastlp) << ".";
109  logger.info(msg2);
110 
111  if (std::fabs(lp - lastlp) <= 1e-8)
112  break;
113  }
114 
115  {
116  std::vector<double> values;
117  std::stringstream ss;
118  model.write_array(rng, cont_vector, disc_vector, values,
119  true, true, &ss);
120  if (ss.str().length() > 0)
121  logger.info(ss);
122  values.insert(values.begin(), lp);
123  parameter_writer(values);
124  }
125  return error_codes::OK;
126  }
127 
128  }
129  }
130 }
131 #endif
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t ss
Definition: plot.C:24
double newton_step(M &model, std::vector< double > &params_r, std::vector< int > &params_i, std::ostream *output_stream=0)
Definition: newton.hpp:32
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
rosenbrock_model_namespace::rosenbrock_model Model
chain
Check that an output directory exists.
int newton(Model &model, stan::io::var_context &init, unsigned int random_seed, unsigned int chain, double init_radius, int num_iterations, bool save_iterations, callbacks::interrupt &interrupt, callbacks::logger &logger, callbacks::writer &init_writer, callbacks::writer &parameter_writer)
Definition: newton.hpp:42
Float_t e
Definition: plot.C:35
const XML_Char XML_Content * model
Definition: expat.h:151
boost::ecuyer1988 create_rng(unsigned int seed, unsigned int chain)
Definition: create_rng.hpp:25