text: // Generated by rstantools.  Do not edit by hand.
  
  /*
      bmsr is free software: you can redistribute it and/or modify
      it under the terms of the GNU General Public License as published by
      the Free Software Foundation, either version 3 of the License, or
      (at your option) any later version.
  
      bmsr is distributed in the hope that it will be useful,
      but WITHOUT ANY WARRANTY; without even the implied warranty of
      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      GNU General Public License for more details.
  
      You should have received a copy of the GNU General Public License
      along with bmsr.  If not, see <http://www.gnu.org/licenses/>.
  */
  #ifndef MODELS_HPP
  #define MODELS_HPP
  #define STAN__SERVICES__COMMAND_HPP
  #include <rstan/rstaninc.hpp>
  // Code generated by Stan version 2.21.0
  #include <stan/model/model_header.hpp>
  namespace model_bmsmtr_namespace {
  using std::istream;
  using std::string;
  using std::stringstream;
  using std::vector;
  using stan::io::dump;
  using stan::math::lgamma;
  using stan::model::prob_grad;
  using namespace stan::math;
  static int current_statement_begin__;
  stan::io::program_reader prog_reader__() {
      stan::io::program_reader reader;
      reader.add_event(0, 0, "start", "model_bmsmtr");
      reader.add_event(71, 69, "end", "model_bmsmtr");
      return reader;
  }
  #include <stan_meta_header.hpp>
  class model_bmsmtr
    : public stan::model::model_base_crtp<model_bmsmtr> {
  private:
          int S;
          std::vector<int> nY;
          int dX;
          int dY;
          matrix_d Y;
          matrix_d X;
          int p0;
          int N;
          std::vector<std::vector<int> > iY;
  public:
      model_bmsmtr(stan::io::var_context& context__,
          std::ostream* pstream__ = 0)
          : model_base_crtp(0) {
          ctor_body(context__, 0, pstream__);
      }
      model_bmsmtr(stan::io::var_context& context__,
          unsigned int random_seed__,
          std::ostream* pstream__ = 0)
          : model_base_crtp(0) {
          ctor_body(context__, random_seed__, pstream__);
      }
      void ctor_body(stan::io::var_context& context__,
                     unsigned int random_seed__,
                     std::ostream* pstream__) {
          typedef double local_scalar_t__;
          boost::ecuyer1988 base_rng__ =
            stan::services::util::create_rng(random_seed__, 0);
          (void) base_rng__;  // suppress unused var warning
          current_statement_begin__ = -1;
          static const char* function__ = "model_bmsmtr_namespace::model_bmsmtr";
          (void) function__;  // dummy to suppress unused var warning
          size_t pos__;
          (void) pos__;  // dummy to suppress unused var warning
          std::vector<int> vals_i__;
          std::vector<double> vals_r__;
          local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
          (void) DUMMY_VAR__;  // suppress unused var warning
          try {
              // initialize data block variables from context__
              current_statement_begin__ = 10;
              context__.validate_dims("data initialization", "S", "int", context__.to_vec());
              S = int(0);
              vals_i__ = context__.vals_i("S");
              pos__ = 0;
              S = vals_i__[pos__++];
              check_greater_or_equal(function__, "S", S, 0);
              current_statement_begin__ = 11;
              validate_non_negative_index("nY", "S", S);
              context__.validate_dims("data initialization", "nY", "int", context__.to_vec(S));
              nY = std::vector<int>(S, int(0));
              vals_i__ = context__.vals_i("nY");
              pos__ = 0;
              size_t nY_k_0_max__ = S;
              for (size_t k_0__ = 0; k_0__ < nY_k_0_max__; ++k_0__) {
                  nY[k_0__] = vals_i__[pos__++];
              }
              size_t nY_i_0_max__ = S;
              for (size_t i_0__ = 0; i_0__ < nY_i_0_max__; ++i_0__) {
                  check_greater_or_equal(function__, "nY[i_0__]", nY[i_0__], 0);
              }
              current_statement_begin__ = 12;
              context__.validate_dims("data initialization", "dX", "int", context__.to_vec());
              dX = int(0);
              vals_i__ = context__.vals_i("dX");
              pos__ = 0;
              dX = vals_i__[pos__++];
              check_greater_or_equal(function__, "dX", dX, 0);
              current_statement_begin__ = 13;
              context__.validate_dims("data initialization", "dY", "int", context__.to_vec());
              dY = int(0);
              vals_i__ = context__.vals_i("dY");
              pos__ = 0;
              dY = vals_i__[pos__++];
              check_greater_or_equal(function__, "dY", dY, 0);
              current_statement_begin__ = 14;
              validate_non_negative_index("Y", "sum(nY)", sum(nY));
              validate_non_negative_index("Y", "dY", dY);
              context__.validate_dims("data initialization", "Y", "matrix_d", context__.to_vec(sum(nY),dY));
              Y = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(sum(nY), dY);
              vals_r__ = context__.vals_r("Y");
              pos__ = 0;
              size_t Y_j_2_max__ = dY;
              size_t Y_j_1_max__ = sum(nY);
              for (size_t j_2__ = 0; j_2__ < Y_j_2_max__; ++j_2__) {
                  for (size_t j_1__ = 0; j_1__ < Y_j_1_max__; ++j_1__) {
                      Y(j_1__, j_2__) = vals_r__[pos__++];
                  }
              }
              current_statement_begin__ = 15;
              validate_non_negative_index("X", "sum(nY)", sum(nY));
              validate_non_negative_index("X", "dX", dX);
              context__.validate_dims("data initialization", "X", "matrix_d", context__.to_vec(sum(nY),dX));
              X = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(sum(nY), dX);
              vals_r__ = context__.vals_r("X");
              pos__ = 0;
              size_t X_j_2_max__ = dX;
              size_t X_j_1_max__ = sum(nY);
              for (size_t j_2__ = 0; j_2__ < X_j_2_max__; ++j_2__) {
                  for (size_t j_1__ = 0; j_1__ < X_j_1_max__; ++j_1__) {
                      X(j_1__, j_2__) = vals_r__[pos__++];
                  }
              }
              current_statement_begin__ = 16;
              context__.validate_dims("data initialization", "p0", "int", context__.to_vec());
              p0 = int(0);
              vals_i__ = context__.vals_i("p0");
              pos__ = 0;
              p0 = vals_i__[pos__++];
              check_greater_or_equal(function__, "p0", p0, 0);
              // initialize transformed data variables
              current_statement_begin__ = 20;
              N = int(0);
              stan::math::fill(N, std::numeric_limits<int>::min());
              current_statement_begin__ = 21;
              validate_non_negative_index("iY", "S", S);
              validate_non_negative_index("iY", "2", 2);
              iY = std::vector<std::vector<int> >(S, std::vector<int>(2, int(0)));
              stan::math::fill(iY, std::numeric_limits<int>::min());
              // execute transformed data statements
              current_statement_begin__ = 23;
              stan::math::assign(N, sum(nY));
              current_statement_begin__ = 24;
              for (int s = 1; s <= S; ++s) {
                  current_statement_begin__ = 25;
                  if (as_bool(logical_eq(s, 1))) {
                      current_statement_begin__ = 26;
                      stan::model::assign(iY, 
                                  stan::model::cons_list(stan::model::index_uni(s), stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list())), 
                                  1, 
                                  "assigning variable iY");
                      current_statement_begin__ = 27;
                      stan::model::assign(iY, 
                                  stan::model::cons_list(stan::model::index_uni(s), stan::model::cons_list(stan::model::index_uni(2), stan::model::nil_index_list())), 
                                  get_base1(nY, s, "nY", 1), 
                                  "assigning variable iY");
                  } else {
                      current_statement_begin__ = 30;
                      stan::model::assign(iY, 
                                  stan::model::cons_list(stan::model::index_uni(s), stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list())), 
                                  (get_base1(get_base1(iY, (s - 1), "iY", 1), 2, "iY", 2) + 1), 
                                  "assigning variable iY");
                      current_statement_begin__ = 31;
                      stan::model::assign(iY, 
                                  stan::model::cons_list(stan::model::index_uni(s), stan::model::cons_list(stan::model::index_uni(2), stan::model::nil_index_list())), 
                                  (get_base1(get_base1(iY, (s - 1), "iY", 1), 2, "iY", 2) + get_base1(nY, s, "nY", 1)), 
                                  "assigning variable iY");
                  }
              }
              // validate transformed data
              current_statement_begin__ = 20;
              check_greater_or_equal(function__, "N", N, 0);
              current_statement_begin__ = 21;
              size_t iY_i_0_max__ = S;
              size_t iY_i_1_max__ = 2;
              for (size_t i_0__ = 0; i_0__ < iY_i_0_max__; ++i_0__) {
                  for (size_t i_1__ = 0; i_1__ < iY_i_1_max__; ++i_1__) {
                      check_greater_or_equal(function__, "iY[i_0__][i_1__]", iY[i_0__][i_1__], 0);
                  }
              }
              // validate, set parameter ranges
              num_params_r__ = 0U;
              param_ranges_i__.clear();
              current_statement_begin__ = 37;
              validate_non_negative_index("W", "dY", dY);
              num_params_r__ += dY;
              current_statement_begin__ = 38;
              validate_non_negative_index("betaM", "dX", dX);
              num_params_r__ += dX;
              current_statement_begin__ = 39;
              num_params_r__ += 1;
              current_statement_begin__ = 40;
              validate_non_negative_index("beta", "S", S);
              validate_non_negative_index("beta", "dX", dX);
              num_params_r__ += (S * dX);
              current_statement_begin__ = 41;
              validate_non_negative_index("lambda", "dX", dX);
              num_params_r__ += dX;
              current_statement_begin__ = 42;
              validate_non_negative_index("sigma", "S", S);
              num_params_r__ += (1 * S);
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
              // Next line prevents compiler griping about no return
              throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
          }
      }
      ~model_bmsmtr() { }
      void transform_inits(const stan::io::var_context& context__,
                           std::vector<int>& params_i__,
                           std::vector<double>& params_r__,
                           std::ostream* pstream__) const {
          typedef double local_scalar_t__;
          stan::io::writer<double> writer__(params_r__, params_i__);
          size_t pos__;
          (void) pos__; // dummy call to supress warning
          std::vector<double> vals_r__;
          std::vector<int> vals_i__;
          current_statement_begin__ = 37;
          if (!(context__.contains_r("W")))
              stan::lang::rethrow_located(std::runtime_error(std::string("Variable W missing")), current_statement_begin__, prog_reader__());
          vals_r__ = context__.vals_r("W");
          pos__ = 0U;
          validate_non_negative_index("W", "dY", dY);
          context__.validate_dims("parameter initialization", "W", "vector_d", context__.to_vec(dY));
          Eigen::Matrix<double, Eigen::Dynamic, 1> W(dY);
          size_t W_j_1_max__ = dY;
          for (size_t j_1__ = 0; j_1__ < W_j_1_max__; ++j_1__) {
              W(j_1__) = vals_r__[pos__++];
          }
          try {
              writer__.vector_unconstrain(W);
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable W: ") + e.what()), current_statement_begin__, prog_reader__());
          }
          current_statement_begin__ = 38;
          if (!(context__.contains_r("betaM")))
              stan::lang::rethrow_located(std::runtime_error(std::string("Variable betaM missing")), current_statement_begin__, prog_reader__());
          vals_r__ = context__.vals_r("betaM");
          pos__ = 0U;
          validate_non_negative_index("betaM", "dX", dX);
          context__.validate_dims("parameter initialization", "betaM", "vector_d", context__.to_vec(dX));
          Eigen::Matrix<double, Eigen::Dynamic, 1> betaM(dX);
          size_t betaM_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < betaM_j_1_max__; ++j_1__) {
              betaM(j_1__) = vals_r__[pos__++];
          }
          try {
              writer__.vector_unconstrain(betaM);
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable betaM: ") + e.what()), current_statement_begin__, prog_reader__());
          }
          current_statement_begin__ = 39;
          if (!(context__.contains_r("tauM")))
              stan::lang::rethrow_located(std::runtime_error(std::string("Variable tauM missing")), current_statement_begin__, prog_reader__());
          vals_r__ = context__.vals_r("tauM");
          pos__ = 0U;
          context__.validate_dims("parameter initialization", "tauM", "double", context__.to_vec());
          double tauM(0);
          tauM = vals_r__[pos__++];
          try {
              writer__.scalar_lb_unconstrain(0, tauM);
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable tauM: ") + e.what()), current_statement_begin__, prog_reader__());
          }
          current_statement_begin__ = 40;
          if (!(context__.contains_r("beta")))
              stan::lang::rethrow_located(std::runtime_error(std::string("Variable beta missing")), current_statement_begin__, prog_reader__());
          vals_r__ = context__.vals_r("beta");
          pos__ = 0U;
          validate_non_negative_index("beta", "S", S);
          validate_non_negative_index("beta", "dX", dX);
          context__.validate_dims("parameter initialization", "beta", "matrix_d", context__.to_vec(S,dX));
          Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> beta(S, dX);
          size_t beta_j_2_max__ = dX;
          size_t beta_j_1_max__ = S;
          for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) {
              for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) {
                  beta(j_1__, j_2__) = vals_r__[pos__++];
              }
          }
          try {
              writer__.matrix_unconstrain(beta);
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable beta: ") + e.what()), current_statement_begin__, prog_reader__());
          }
          current_statement_begin__ = 41;
          if (!(context__.contains_r("lambda")))
              stan::lang::rethrow_located(std::runtime_error(std::string("Variable lambda missing")), current_statement_begin__, prog_reader__());
          vals_r__ = context__.vals_r("lambda");
          pos__ = 0U;
          validate_non_negative_index("lambda", "dX", dX);
          context__.validate_dims("parameter initialization", "lambda", "vector_d", context__.to_vec(dX));
          Eigen::Matrix<double, Eigen::Dynamic, 1> lambda(dX);
          size_t lambda_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < lambda_j_1_max__; ++j_1__) {
              lambda(j_1__) = vals_r__[pos__++];
          }
          try {
              writer__.vector_lb_unconstrain(0, lambda);
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable lambda: ") + e.what()), current_statement_begin__, prog_reader__());
          }
          current_statement_begin__ = 42;
          if (!(context__.contains_r("sigma")))
              stan::lang::rethrow_located(std::runtime_error(std::string("Variable sigma missing")), current_statement_begin__, prog_reader__());
          vals_r__ = context__.vals_r("sigma");
          pos__ = 0U;
          validate_non_negative_index("sigma", "S", S);
          context__.validate_dims("parameter initialization", "sigma", "double", context__.to_vec(S));
          std::vector<double> sigma(S, double(0));
          size_t sigma_k_0_max__ = S;
          for (size_t k_0__ = 0; k_0__ < sigma_k_0_max__; ++k_0__) {
              sigma[k_0__] = vals_r__[pos__++];
          }
          size_t sigma_i_0_max__ = S;
          for (size_t i_0__ = 0; i_0__ < sigma_i_0_max__; ++i_0__) {
              try {
                  writer__.scalar_lb_unconstrain(0, sigma[i_0__]);
              } catch (const std::exception& e) {
                  stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable sigma: ") + e.what()), current_statement_begin__, prog_reader__());
              }
          }
          params_r__ = writer__.data_r();
          params_i__ = writer__.data_i();
      }
      void transform_inits(const stan::io::var_context& context,
                           Eigen::Matrix<double, Eigen::Dynamic, 1>& params_r,
                           std::ostream* pstream__) const {
        std::vector<double> params_r_vec;
        std::vector<int> params_i_vec;
        transform_inits(context, params_i_vec, params_r_vec, pstream__);
        params_r.resize(params_r_vec.size());
        for (int i = 0; i < params_r.size(); ++i)
          params_r(i) = params_r_vec[i];
      }
      template <bool propto__, bool jacobian__, typename T__>
      T__ log_prob(std::vector<T__>& params_r__,
                   std::vector<int>& params_i__,
                   std::ostream* pstream__ = 0) const {
          typedef T__ local_scalar_t__;
          local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
          (void) DUMMY_VAR__;  // dummy to suppress unused var warning
          T__ lp__(0.0);
          stan::math::accumulator<T__> lp_accum__;
          try {
              stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
              // model parameters
              current_statement_begin__ = 37;
              Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> W;
              (void) W;  // dummy to suppress unused var warning
              if (jacobian__)
                  W = in__.vector_constrain(dY, lp__);
              else
                  W = in__.vector_constrain(dY);
              current_statement_begin__ = 38;
              Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> betaM;
              (void) betaM;  // dummy to suppress unused var warning
              if (jacobian__)
                  betaM = in__.vector_constrain(dX, lp__);
              else
                  betaM = in__.vector_constrain(dX);
              current_statement_begin__ = 39;
              local_scalar_t__ tauM;
              (void) tauM;  // dummy to suppress unused var warning
              if (jacobian__)
                  tauM = in__.scalar_lb_constrain(0, lp__);
              else
                  tauM = in__.scalar_lb_constrain(0);
              current_statement_begin__ = 40;
              Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> beta;
              (void) beta;  // dummy to suppress unused var warning
              if (jacobian__)
                  beta = in__.matrix_constrain(S, dX, lp__);
              else
                  beta = in__.matrix_constrain(S, dX);
              current_statement_begin__ = 41;
              Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> lambda;
              (void) lambda;  // dummy to suppress unused var warning
              if (jacobian__)
                  lambda = in__.vector_lb_constrain(0, dX, lp__);
              else
                  lambda = in__.vector_lb_constrain(0, dX);
              current_statement_begin__ = 42;
              std::vector<local_scalar_t__> sigma;
              size_t sigma_d_0_max__ = S;
              sigma.reserve(sigma_d_0_max__);
              for (size_t d_0__ = 0; d_0__ < sigma_d_0_max__; ++d_0__) {
                  if (jacobian__)
                      sigma.push_back(in__.scalar_lb_constrain(0, lp__));
                  else
                      sigma.push_back(in__.scalar_lb_constrain(0));
              }
              // transformed parameters
              current_statement_begin__ = 46;
              local_scalar_t__ tau0M;
              (void) tau0M;  // dummy to suppress unused var warning
              stan::math::initialize(tau0M, DUMMY_VAR__);
              stan::math::fill(tau0M, DUMMY_VAR__);
              current_statement_begin__ = 47;
              validate_non_negative_index("betaT", "S", S);
              validate_non_negative_index("betaT", "dX", dX);
              Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> betaT(S, dX);
              stan::math::initialize(betaT, DUMMY_VAR__);
              stan::math::fill(betaT, DUMMY_VAR__);
              // transformed parameters block statements
              current_statement_begin__ = 48;
              stan::math::assign(tau0M, ((sum(sigma) * p0) / ((dX - p0) * stan::math::sqrt((N * 1.0)))));
              current_statement_begin__ = 49;
              for (int s = 1; s <= S; ++s) {
                  current_statement_begin__ = 51;
                  stan::model::assign(betaT, 
                              stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                              elt_multiply(add(elt_divide(get_base1(beta, s, "beta", 1), 2), 0.5), transpose(multiply(multiply(elt_multiply(betaM, lambda), tauM), tau0M))), 
                              "assigning variable betaT");
              }
              // validate transformed parameters
              const char* function__ = "validate transformed params";
              (void) function__;  // dummy to suppress unused var warning
              current_statement_begin__ = 46;
              if (stan::math::is_uninitialized(tau0M)) {
                  std::stringstream msg__;
                  msg__ << "Undefined transformed parameter: tau0M";
                  stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable tau0M: ") + msg__.str()), current_statement_begin__, prog_reader__());
              }
              check_greater_or_equal(function__, "tau0M", tau0M, 0);
              current_statement_begin__ = 47;
              size_t betaT_j_1_max__ = S;
              size_t betaT_j_2_max__ = dX;
              for (size_t j_1__ = 0; j_1__ < betaT_j_1_max__; ++j_1__) {
                  for (size_t j_2__ = 0; j_2__ < betaT_j_2_max__; ++j_2__) {
                      if (stan::math::is_uninitialized(betaT(j_1__, j_2__))) {
                          std::stringstream msg__;
                          msg__ << "Undefined transformed parameter: betaT" << "(" << j_1__ << ", " << j_2__ << ")";
                          stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable betaT: ") + msg__.str()), current_statement_begin__, prog_reader__());
                      }
                  }
              }
              // model body
              current_statement_begin__ = 56;
              lp_accum__.add(normal_log<propto__>(W, 0, 1));
              current_statement_begin__ = 57;
              lp_accum__.add(cauchy_log<propto__>(lambda, 0, 1));
              current_statement_begin__ = 58;
              lp_accum__.add(inv_gamma_log<propto__>(sigma, 1, 1));
              current_statement_begin__ = 59;
              lp_accum__.add(cauchy_log<propto__>(tauM, 0, 1));
              current_statement_begin__ = 60;
              lp_accum__.add(normal_log<propto__>(betaM, 0, 1));
              current_statement_begin__ = 62;
              for (int s = 1; s <= S; ++s) {
                  current_statement_begin__ = 64;
                  lp_accum__.add(normal_log<propto__>(get_base1(beta, s, "beta", 1), 0, 1));
                  current_statement_begin__ = 65;
                  for (int d = 1; d <= dY; ++d) {
                      current_statement_begin__ = 66;
                      lp_accum__.add(normal_log<propto__>(stan::model::rvalue(Y, stan::model::cons_list(stan::model::index_min_max(get_base1(get_base1(iY, s, "iY", 1), 1, "iY", 2), get_base1(get_base1(iY, s, "iY", 1), 2, "iY", 2)), stan::model::cons_list(stan::model::index_uni(d), stan::model::nil_index_list())), "Y"), multiply(multiply(stan::model::rvalue(X, stan::model::cons_list(stan::model::index_min_max(get_base1(get_base1(iY, s, "iY", 1), 1, "iY", 2), get_base1(get_base1(iY, s, "iY", 1), 2, "iY", 2)), stan::model::nil_index_list()), "X"), transpose(get_base1(betaT, s, "betaT", 1))), (0.5 + (get_base1(W, d, "W", 1) / 2))), get_base1(sigma, s, "sigma", 1)));
                  }
              }
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
              // Next line prevents compiler griping about no return
              throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
          }
          lp_accum__.add(lp__);
          return lp_accum__.sum();
      } // log_prob()
      template <bool propto, bool jacobian, typename T_>
      T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
                 std::ostream* pstream = 0) const {
        std::vector<T_> vec_params_r;
        vec_params_r.reserve(params_r.size());
        for (int i = 0; i < params_r.size(); ++i)
          vec_params_r.push_back(params_r(i));
        std::vector<int> vec_params_i;
        return log_prob<propto,jacobian,T_>(vec_params_r, vec_params_i, pstream);
      }
      void get_param_names(std::vector<std::string>& names__) const {
          names__.resize(0);
          names__.push_back("W");
          names__.push_back("betaM");
          names__.push_back("tauM");
          names__.push_back("beta");
          names__.push_back("lambda");
          names__.push_back("sigma");
          names__.push_back("tau0M");
          names__.push_back("betaT");
      }
      void get_dims(std::vector<std::vector<size_t> >& dimss__) const {
          dimss__.resize(0);
          std::vector<size_t> dims__;
          dims__.resize(0);
          dims__.push_back(dY);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dims__.push_back(dX);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dims__.push_back(S);
          dims__.push_back(dX);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dims__.push_back(dX);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dims__.push_back(S);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dimss__.push_back(dims__);
          dims__.resize(0);
          dims__.push_back(S);
          dims__.push_back(dX);
          dimss__.push_back(dims__);
      }
      template <typename RNG>
      void write_array(RNG& base_rng__,
                       std::vector<double>& params_r__,
                       std::vector<int>& params_i__,
                       std::vector<double>& vars__,
                       bool include_tparams__ = true,
                       bool include_gqs__ = true,
                       std::ostream* pstream__ = 0) const {
          typedef double local_scalar_t__;
          vars__.resize(0);
          stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
          static const char* function__ = "model_bmsmtr_namespace::write_array";
          (void) function__;  // dummy to suppress unused var warning
          // read-transform, write parameters
          Eigen::Matrix<double, Eigen::Dynamic, 1> W = in__.vector_constrain(dY);
          size_t W_j_1_max__ = dY;
          for (size_t j_1__ = 0; j_1__ < W_j_1_max__; ++j_1__) {
              vars__.push_back(W(j_1__));
          }
          Eigen::Matrix<double, Eigen::Dynamic, 1> betaM = in__.vector_constrain(dX);
          size_t betaM_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < betaM_j_1_max__; ++j_1__) {
              vars__.push_back(betaM(j_1__));
          }
          double tauM = in__.scalar_lb_constrain(0);
          vars__.push_back(tauM);
          Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> beta = in__.matrix_constrain(S, dX);
          size_t beta_j_2_max__ = dX;
          size_t beta_j_1_max__ = S;
          for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) {
              for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) {
                  vars__.push_back(beta(j_1__, j_2__));
              }
          }
          Eigen::Matrix<double, Eigen::Dynamic, 1> lambda = in__.vector_lb_constrain(0, dX);
          size_t lambda_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < lambda_j_1_max__; ++j_1__) {
              vars__.push_back(lambda(j_1__));
          }
          std::vector<double> sigma;
          size_t sigma_d_0_max__ = S;
          sigma.reserve(sigma_d_0_max__);
          for (size_t d_0__ = 0; d_0__ < sigma_d_0_max__; ++d_0__) {
              sigma.push_back(in__.scalar_lb_constrain(0));
          }
          size_t sigma_k_0_max__ = S;
          for (size_t k_0__ = 0; k_0__ < sigma_k_0_max__; ++k_0__) {
              vars__.push_back(sigma[k_0__]);
          }
          double lp__ = 0.0;
          (void) lp__;  // dummy to suppress unused var warning
          stan::math::accumulator<double> lp_accum__;
          local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
          (void) DUMMY_VAR__;  // suppress unused var warning
          if (!include_tparams__ && !include_gqs__) return;
          try {
              // declare and define transformed parameters
              current_statement_begin__ = 46;
              double tau0M;
              (void) tau0M;  // dummy to suppress unused var warning
              stan::math::initialize(tau0M, DUMMY_VAR__);
              stan::math::fill(tau0M, DUMMY_VAR__);
              current_statement_begin__ = 47;
              validate_non_negative_index("betaT", "S", S);
              validate_non_negative_index("betaT", "dX", dX);
              Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> betaT(S, dX);
              stan::math::initialize(betaT, DUMMY_VAR__);
              stan::math::fill(betaT, DUMMY_VAR__);
              // do transformed parameters statements
              current_statement_begin__ = 48;
              stan::math::assign(tau0M, ((sum(sigma) * p0) / ((dX - p0) * stan::math::sqrt((N * 1.0)))));
              current_statement_begin__ = 49;
              for (int s = 1; s <= S; ++s) {
                  current_statement_begin__ = 51;
                  stan::model::assign(betaT, 
                              stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                              elt_multiply(add(elt_divide(get_base1(beta, s, "beta", 1), 2), 0.5), transpose(multiply(multiply(elt_multiply(betaM, lambda), tauM), tau0M))), 
                              "assigning variable betaT");
              }
              if (!include_gqs__ && !include_tparams__) return;
              // validate transformed parameters
              const char* function__ = "validate transformed params";
              (void) function__;  // dummy to suppress unused var warning
              current_statement_begin__ = 46;
              check_greater_or_equal(function__, "tau0M", tau0M, 0);
              // write transformed parameters
              if (include_tparams__) {
                  vars__.push_back(tau0M);
                  size_t betaT_j_2_max__ = dX;
                  size_t betaT_j_1_max__ = S;
                  for (size_t j_2__ = 0; j_2__ < betaT_j_2_max__; ++j_2__) {
                      for (size_t j_1__ = 0; j_1__ < betaT_j_1_max__; ++j_1__) {
                          vars__.push_back(betaT(j_1__, j_2__));
                      }
                  }
              }
              if (!include_gqs__) return;
          } catch (const std::exception& e) {
              stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
              // Next line prevents compiler griping about no return
              throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
          }
      }
      template <typename RNG>
      void write_array(RNG& base_rng,
                       Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,
                       Eigen::Matrix<double,Eigen::Dynamic,1>& vars,
                       bool include_tparams = true,
                       bool include_gqs = true,
                       std::ostream* pstream = 0) const {
        std::vector<double> params_r_vec(params_r.size());
        for (int i = 0; i < params_r.size(); ++i)
          params_r_vec[i] = params_r(i);
        std::vector<double> vars_vec;
        std::vector<int> params_i_vec;
        write_array(base_rng, params_r_vec, params_i_vec, vars_vec, include_tparams, include_gqs, pstream);
        vars.resize(vars_vec.size());
        for (int i = 0; i < vars.size(); ++i)
          vars(i) = vars_vec[i];
      }
      std::string model_name() const {
          return "model_bmsmtr";
      }
      void constrained_param_names(std::vector<std::string>& param_names__,
                                   bool include_tparams__ = true,
                                   bool include_gqs__ = true) const {
          std::stringstream param_name_stream__;
          size_t W_j_1_max__ = dY;
          for (size_t j_1__ = 0; j_1__ < W_j_1_max__; ++j_1__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "W" << '.' << j_1__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          size_t betaM_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < betaM_j_1_max__; ++j_1__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "betaM" << '.' << j_1__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          param_name_stream__.str(std::string());
          param_name_stream__ << "tauM";
          param_names__.push_back(param_name_stream__.str());
          size_t beta_j_2_max__ = dX;
          size_t beta_j_1_max__ = S;
          for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) {
              for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) {
                  param_name_stream__.str(std::string());
                  param_name_stream__ << "beta" << '.' << j_1__ + 1 << '.' << j_2__ + 1;
                  param_names__.push_back(param_name_stream__.str());
              }
          }
          size_t lambda_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < lambda_j_1_max__; ++j_1__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "lambda" << '.' << j_1__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          size_t sigma_k_0_max__ = S;
          for (size_t k_0__ = 0; k_0__ < sigma_k_0_max__; ++k_0__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "sigma" << '.' << k_0__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          if (!include_gqs__ && !include_tparams__) return;
          if (include_tparams__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "tau0M";
              param_names__.push_back(param_name_stream__.str());
              size_t betaT_j_2_max__ = dX;
              size_t betaT_j_1_max__ = S;
              for (size_t j_2__ = 0; j_2__ < betaT_j_2_max__; ++j_2__) {
                  for (size_t j_1__ = 0; j_1__ < betaT_j_1_max__; ++j_1__) {
                      param_name_stream__.str(std::string());
                      param_name_stream__ << "betaT" << '.' << j_1__ + 1 << '.' << j_2__ + 1;
                      param_names__.push_back(param_name_stream__.str());
                  }
              }
          }
          if (!include_gqs__) return;
      }
      void unconstrained_param_names(std::vector<std::string>& param_names__,
                                     bool include_tparams__ = true,
                                     bool include_gqs__ = true) const {
          std::stringstream param_name_stream__;
          size_t W_j_1_max__ = dY;
          for (size_t j_1__ = 0; j_1__ < W_j_1_max__; ++j_1__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "W" << '.' << j_1__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          size_t betaM_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < betaM_j_1_max__; ++j_1__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "betaM" << '.' << j_1__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          param_name_stream__.str(std::string());
          param_name_stream__ << "tauM";
          param_names__.push_back(param_name_stream__.str());
          size_t beta_j_2_max__ = dX;
          size_t beta_j_1_max__ = S;
          for (size_t j_2__ = 0; j_2__ < beta_j_2_max__; ++j_2__) {
              for (size_t j_1__ = 0; j_1__ < beta_j_1_max__; ++j_1__) {
                  param_name_stream__.str(std::string());
                  param_name_stream__ << "beta" << '.' << j_1__ + 1 << '.' << j_2__ + 1;
                  param_names__.push_back(param_name_stream__.str());
              }
          }
          size_t lambda_j_1_max__ = dX;
          for (size_t j_1__ = 0; j_1__ < lambda_j_1_max__; ++j_1__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "lambda" << '.' << j_1__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          size_t sigma_k_0_max__ = S;
          for (size_t k_0__ = 0; k_0__ < sigma_k_0_max__; ++k_0__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "sigma" << '.' << k_0__ + 1;
              param_names__.push_back(param_name_stream__.str());
          }
          if (!include_gqs__ && !include_tparams__) return;
          if (include_tparams__) {
              param_name_stream__.str(std::string());
              param_name_stream__ << "tau0M";
              param_names__.push_back(param_name_stream__.str());
              size_t betaT_j_2_max__ = dX;
              size_t betaT_j_1_max__ = S;
              for (size_t j_2__ = 0; j_2__ < betaT_j_2_max__; ++j_2__) {
                  for (size_t j_1__ = 0; j_1__ < betaT_j_1_max__; ++j_1__) {
                      param_name_stream__.str(std::string());
                      param_name_stream__ << "betaT" << '.' << j_1__ + 1 << '.' << j_2__ + 1;
                      param_names__.push_back(param_name_stream__.str());
                  }
              }
          }
          if (!include_gqs__) return;
      }
  }; // model
  }  // namespace
  typedef model_bmsmtr_namespace::model_bmsmtr stan_model;
  #ifndef USING_R
  stan::model::model_base& new_model(
          stan::io::var_context& data_context,
          unsigned int seed,
          std::ostream* msg_stream) {
    stan_model* m = new stan_model(data_context, seed, msg_stream);
    return *m;
  }
  #endif
  #endif
  