Update log_prob method to take constrained_params in init format

src/cmdstan/arguments/arg_log_prob.hpp

@@ -8,6 +8,11 @@
 
 namespace cmdstan {
 
+/**
+ * Argument used for calling the log_prob method to calculate the log density
+ * and its gradients with respect to a user-provided set of parameter values
+ * on the constrained and/or unconstrained scale
+ */
 class arg_log_prob : public categorical_argument {
  public:
   arg_log_prob() {

src/cmdstan/arguments/arg_log_prob_constrained_params.hpp

@@ -6,9 +6,10 @@
 namespace cmdstan {
 /**
  * Argument for providing a file of parameters on the constrained scale
- * for use with the 'log_prob' method. The file can be in CSV or JSON format
- * and should contain a variable 'params_r' with either a vector or list/array
- * of vectors of constrained parameter values.
+ * for use with the 'log_prob' method. The file can be in JSON or R Dump format,
+ * using the same structure as the 'init' argument. Like the 'init' argument, if
+ * the file has a '.json' extension it is treated as a JSON file, otherwise it
+ * is treated as an RDump file.
  */
 class arg_log_prob_constrained_params : public string_argument {
  public:

src/cmdstan/arguments/arg_log_prob_output_file.hpp

@@ -4,7 +4,10 @@
 #include <cmdstan/arguments/singleton_argument.hpp>
 
 namespace cmdstan {
-
+/**
+ * Sub-argument to provide a custom path and filename for the output from
+ * the log_prob method
+ */
 class arg_log_prob_output_file : public string_argument {
  public:
   arg_log_prob_output_file() : string_argument() {

src/cmdstan/arguments/arg_log_prob_unconstrained_params.hpp

@@ -8,7 +8,9 @@ namespace cmdstan {
  * Argument for providing a file of parameters on the unconstrained scale
  * for use with the 'log_prob' method. The file can be in CSV or JSON format
  * and should contain a variable 'params_r' with either a vector or list/array
- * of vectors of unconstrained parameter values.
+ * of vectors of unconstrained parameter values. Like the 'init' argument, if
+ * the file has a '.json' extension it is treated as a JSON file, otherwise it
+ * is treated as an RDump file.
  */
 class arg_log_prob_unconstrained_params : public string_argument {
  public:

src/cmdstan/command.hpp

@@ -485,7 +485,7 @@ int command(int argc, const char *argv[]) {
       return std::string("_" + std::to_string(i + id));
     }
   };
-  for (int i = 0; i < num_chains; i++) {
+  for (int i = 0; i < num_chains; ++i) {
     auto output_filename = output_name + name_iterator(i) + output_ending;
     auto unique_fstream
         = std::make_unique<std::fstream>(output_filename, std::fstream::out);
@@ -504,7 +504,7 @@ int command(int argc, const char *argv[]) {
       diagnostic_writers.emplace_back(nullptr, "# ");
     }
   }
-  for (int i = 0; i < num_chains; i++) {
+  for (int i = 0; i < num_chains; ++i) {
     write_stan(sample_writers[i]);
     write_model(sample_writers[i], model.model_name());
     write_datetime(sample_writers[i]);
@@ -620,8 +620,13 @@ int command(int argc, const char *argv[]) {
       std::shared_ptr<stan::io::var_context> upars_context
           = get_var_context(u_fname);
 
-      u_params_r = (*upars_context).vals_r("params_r");
-      dims_u_params_r = (*upars_context).dims_r("params_r");
+      u_params_r = upars_context->vals_r("params_r");
+      if (u_params_r.size() == 0) {
+        msg << "Unconstrained parameters file has no variable 'params_r' with "
+               "unconstrained parameter values!";
+        throw std::invalid_argument(msg.str());
+      }
+      dims_u_params_r = upars_context->dims_r("params_r");
 
       // Detect whether multiple sets of parameter values have been passed
       // and set the sizes accordingly
@@ -644,50 +649,34 @@ int command(int argc, const char *argv[]) {
       throw std::invalid_argument(msg.str());
     }
 
-    size_t c_params_vec_size = 0;
-    size_t c_params_size = 0;
-    std::vector<double> c_params_r;
-    std::vector<size_t> dims_c_params_r;
-    if (!(cpars_file->is_default())) {
-      std::string c_fname(cpars_file->value());
-      std::ifstream c_stream(c_fname.c_str());
-
-      std::shared_ptr<stan::io::var_context> cpars_context
-          = get_var_context(c_fname);
-
-      c_params_r = (*cpars_context).vals_r("params_r");
-      dims_c_params_r = (*cpars_context).dims_r("params_r");
-
-      c_params_vec_size = dims_c_params_r.size() == 2 ? dims_c_params_r[0] : 1;
-      c_params_size = dims_c_params_r.size() == 2 ? dims_c_params_r[1]
-                                                  : dims_c_params_r[0];
-    }
-
-    if (c_params_size > 0 && c_params_size != param_names.size()) {
-      msg << "Incorrect number of constrained parameters provided! "
-             "Model has "
-          << param_names.size() << " parameters but " << c_params_size
-          << " were found.";
-      throw std::invalid_argument(msg.str());
-    }
-
+    int has_cpars = !cpars_file->is_default();
     // Store in single nested array to allow single loop for calc and print
-    size_t num_par_sets = c_params_vec_size + u_params_vec_size;
+    size_t num_par_sets = u_params_vec_size + has_cpars;
     std::vector<std::vector<double>> params_r_ind(num_par_sets);
-
-    // Use Map with inner stride to operate on all values from parameter set
-    using StrideT = Eigen::Stride<1, Eigen::Dynamic>;
     std::vector<int> dummy_params_i;
-    for (size_t i = 0; i < c_params_vec_size; i++) {
-      Eigen::Map<Eigen::VectorXd, 0, StrideT> map_r(
-          c_params_r.data() + i, c_params_size, StrideT(1, c_params_vec_size));
+    if (!cpars_file->is_default()) {
+      std::string cpars_filename(cpars_file->value());
 
-      stan::io::array_var_context context(param_names, map_r, param_dimss);
-      model.transform_inits(context, dummy_params_i, params_r_ind[i], &msg);
+      std::shared_ptr<stan::io::var_context> cpars_context
+          = get_var_context(cpars_filename);
+      std::vector<std::string> input_cpar_names;
+      cpars_context->names_r(input_cpar_names);
+
+      for (const std::string &m_param_name : param_names) {
+        if (!cpars_context->contains_r(m_param_name)) {
+          msg << "Constrained value(s) for parameter " << m_param_name
+              << " not found!";
+          throw std::invalid_argument(msg.str());
+        }
+      }
+      model.transform_inits((*cpars_context), dummy_params_i, params_r_ind[0],
+                            &msg);
     }
 
-    for (size_t i = c_params_vec_size; i < num_par_sets; i++) {
-      size_t iter = i - c_params_vec_size;
+    // Use Map with inner stride to operate on all values from parameter set
+    using StrideT = Eigen::Stride<1, Eigen::Dynamic>;
+    for (size_t i = has_cpars; i < num_par_sets; ++i) {
+      size_t iter = i - has_cpars;
       Eigen::Map<Eigen::VectorXd, 0, StrideT> map_r(
           u_params_r.data() + iter, u_params_size,
           StrideT(1, u_params_vec_size));
@@ -702,14 +691,14 @@ int command(int argc, const char *argv[]) {
 
     std::vector<std::string> p_names;
     model.constrained_param_names(p_names, false, false);
-    for (size_t i = 1; i < p_names.size(); i++) {
+    for (size_t i = 1; i < p_names.size(); ++i) {
       output_stream << "g_" << p_names[i] << ",";
     }
     output_stream << "g_" << p_names.back() << "\n";
     try {
       double lp;
       std::vector<double> gradients;
-      for (size_t i = 0; i < num_par_sets; i++) {
+      for (size_t i = 0; i < num_par_sets; ++i) {
         if (jacobian_adjust) {
           lp = stan::model::log_prob_grad<false, true>(
               model, params_r_ind[i], dummy_params_i, gradients);

src/test/test-models/bern_constrained_params.R

@@ -1,3 +1,3 @@
-params_r <-
-structure(c(0.5, 0.1, -1.5, -3.5, 2.6, 1.6),
-.Dim = c(2, 3))
+theta <- 0.1
+mu1 <- 2.0
+mu2 <- 4.0

src/test/test-models/bern_constrained_params.json

@@ -1 +1,5 @@
-{"params_r":[[0.5, -1.5, 2.6], [0.1, -3.5, 1.6]]}
+{
+  "theta": 0.1,
+  "mu1" : 2.0,
+  "mu2" : 4.0
+}

src/test/test-models/bern_constrained_params_short.json

@@ -1 +1,4 @@
-{"params_r":[[0.5, -1.5, 2.6], [0.1, -3.5]]}
+{
+  "theta": 0.1,
+  "mu1" : 2.0
+}