// Code generated by Stan version 2.17.0

#include <stan/model/model_header.hpp>

namespace model_count_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;

typedef Eigen::Matrix<double,Eigen::Dynamic,1> vector_d;
typedef Eigen::Matrix<double,1,Eigen::Dynamic> row_vector_d;
typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> matrix_d;

static int current_statement_begin__;

stan::io::program_reader prog_reader__() {
    stan::io::program_reader reader;
    reader.add_event(0, 0, "start", "model_count");
    reader.add_event(0, 0, "include", "/pre/Columbia_copyright.stan");
    reader.add_event(0, 0, "start", "/pre/Columbia_copyright.stan");
    reader.add_event(3, 3, "end", "/pre/Columbia_copyright.stan");
    reader.add_event(3, 1, "restart", "model_count");
    reader.add_event(3, 1, "include", "/pre/license.stan");
    reader.add_event(3, 0, "start", "/pre/license.stan");
    reader.add_event(17, 14, "end", "/pre/license.stan");
    reader.add_event(17, 2, "restart", "model_count");
    reader.add_event(20, 5, "include", "/functions/common_functions.stan");
    reader.add_event(20, 0, "start", "/functions/common_functions.stan");
    reader.add_event(318, 298, "end", "/functions/common_functions.stan");
    reader.add_event(318, 6, "restart", "model_count");
    reader.add_event(318, 6, "include", "/functions/count_likelihoods.stan");
    reader.add_event(318, 0, "start", "/functions/count_likelihoods.stan");
    reader.add_event(371, 53, "end", "/functions/count_likelihoods.stan");
    reader.add_event(371, 7, "restart", "model_count");
    reader.add_event(374, 10, "include", "/data/NKX.stan");
    reader.add_event(374, 0, "start", "/data/NKX.stan");
    reader.add_event(394, 20, "end", "/data/NKX.stan");
    reader.add_event(394, 11, "restart", "model_count");
    reader.add_event(396, 13, "include", "/data/data_glm.stan");
    reader.add_event(396, 0, "start", "/data/data_glm.stan");
    reader.add_event(415, 19, "end", "/data/data_glm.stan");
    reader.add_event(415, 14, "restart", "model_count");
    reader.add_event(416, 15, "include", "/data/weights_offset.stan");
    reader.add_event(416, 0, "start", "/data/weights_offset.stan");
    reader.add_event(423, 7, "end", "/data/weights_offset.stan");
    reader.add_event(423, 16, "restart", "model_count");
    reader.add_event(425, 18, "include", "/data/hyperparameters.stan");
    reader.add_event(425, 0, "start", "/data/hyperparameters.stan");
    reader.add_event(443, 18, "end", "/data/hyperparameters.stan");
    reader.add_event(443, 19, "restart", "model_count");
    reader.add_event(444, 20, "include", "/data/glmer_stuff.stan");
    reader.add_event(444, 0, "start", "/data/glmer_stuff.stan");
    reader.add_event(459, 15, "end", "/data/glmer_stuff.stan");
    reader.add_event(459, 21, "restart", "model_count");
    reader.add_event(460, 22, "include", "/data/glmer_stuff2.stan");
    reader.add_event(460, 0, "start", "/data/glmer_stuff2.stan");
    reader.add_event(465, 5, "end", "/data/glmer_stuff2.stan");
    reader.add_event(465, 23, "restart", "model_count");
    reader.add_event(470, 28, "include", "/tdata/tdata_glm.stan");
    reader.add_event(470, 0, "start", "/tdata/tdata_glm.stan");
    reader.add_event(492, 22, "end", "/tdata/tdata_glm.stan");
    reader.add_event(492, 29, "restart", "model_count");
    reader.add_event(496, 33, "include", "/parameters/parameters_glm.stan");
    reader.add_event(496, 0, "start", "/parameters/parameters_glm.stan");
    reader.add_event(509, 13, "end", "/parameters/parameters_glm.stan");
    reader.add_event(509, 34, "restart", "model_count");
    reader.add_event(515, 40, "include", "/tparameters/tparameters_glm.stan");
    reader.add_event(515, 0, "start", "/tparameters/tparameters_glm.stan");
    reader.add_event(559, 44, "end", "/tparameters/tparameters_glm.stan");
    reader.add_event(559, 41, "restart", "model_count");
    reader.add_event(591, 73, "include", "/model/make_eta.stan");
    reader.add_event(591, 0, "start", "/model/make_eta.stan");
    reader.add_event(599, 8, "end", "/model/make_eta.stan");
    reader.add_event(599, 74, "restart", "model_count");
    reader.add_event(600, 75, "include", "/model/eta_add_Zb.stan");
    reader.add_event(600, 0, "start", "/model/eta_add_Zb.stan");
    reader.add_event(602, 2, "end", "/model/eta_add_Zb.stan");
    reader.add_event(602, 76, "restart", "model_count");
    reader.add_event(608, 82, "include", "/model/eta_no_intercept.stan");
    reader.add_event(608, 0, "start", "/model/eta_no_intercept.stan");
    reader.add_event(610, 2, "end", "/model/eta_no_intercept.stan");
    reader.add_event(610, 83, "restart", "model_count");
    reader.add_event(646, 119, "include", "/model/priors_glm.stan");
    reader.add_event(646, 0, "start", "/model/priors_glm.stan");
    reader.add_event(706, 60, "end", "/model/priors_glm.stan");
    reader.add_event(706, 120, "restart", "model_count");
    reader.add_event(722, 136, "include", "/model/make_eta.stan");
    reader.add_event(722, 0, "start", "/model/make_eta.stan");
    reader.add_event(730, 8, "end", "/model/make_eta.stan");
    reader.add_event(730, 137, "restart", "model_count");
    reader.add_event(731, 138, "include", "/model/eta_add_Zb.stan");
    reader.add_event(731, 0, "start", "/model/eta_add_Zb.stan");
    reader.add_event(733, 2, "end", "/model/eta_add_Zb.stan");
    reader.add_event(733, 139, "restart", "model_count");
    reader.add_event(743, 149, "include", "/model/eta_no_intercept.stan");
    reader.add_event(743, 0, "start", "/model/eta_no_intercept.stan");
    reader.add_event(745, 2, "end", "/model/eta_no_intercept.stan");
    reader.add_event(745, 150, "restart", "model_count");
    reader.add_event(768, 173, "end", "model_count");
    return reader;
}

template <typename T2__, typename T3__, typename T4__, typename T5__, typename T6__, typename T7__>
Eigen::Matrix<typename boost::math::tools::promote_args<T2__, T3__, T4__, T5__, typename boost::math::tools::promote_args<T6__, T7__>::type>::type, Eigen::Dynamic,1>
make_theta_L(const int& len_theta_L,
                 const std::vector<int>& p,
                 const T2__& dispersion,
                 const Eigen::Matrix<T3__, Eigen::Dynamic,1>& tau,
                 const Eigen::Matrix<T4__, Eigen::Dynamic,1>& scale,
                 const Eigen::Matrix<T5__, Eigen::Dynamic,1>& zeta,
                 const Eigen::Matrix<T6__, Eigen::Dynamic,1>& rho,
                 const Eigen::Matrix<T7__, Eigen::Dynamic,1>& z_T, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T2__, T3__, T4__, T5__, typename boost::math::tools::promote_args<T6__, T7__>::type>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 42;
        validate_non_negative_index("theta_L", "len_theta_L", len_theta_L);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  theta_L(static_cast<Eigen::VectorXd::Index>(len_theta_L));
        (void) theta_L;  // dummy to suppress unused var warning

        stan::math::initialize(theta_L, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(theta_L,DUMMY_VAR__);
        current_statement_begin__ = 43;
        int zeta_mark(0);
        (void) zeta_mark;  // dummy to suppress unused var warning

        stan::math::fill(zeta_mark, std::numeric_limits<int>::min());
        stan::math::assign(zeta_mark,1);
        current_statement_begin__ = 44;
        int rho_mark(0);
        (void) rho_mark;  // dummy to suppress unused var warning

        stan::math::fill(rho_mark, std::numeric_limits<int>::min());
        stan::math::assign(rho_mark,1);
        current_statement_begin__ = 45;
        int z_T_mark(0);
        (void) z_T_mark;  // dummy to suppress unused var warning

        stan::math::fill(z_T_mark, std::numeric_limits<int>::min());
        stan::math::assign(z_T_mark,1);
        current_statement_begin__ = 46;
        int theta_L_mark(0);
        (void) theta_L_mark;  // dummy to suppress unused var warning

        stan::math::fill(theta_L_mark, std::numeric_limits<int>::min());
        stan::math::assign(theta_L_mark,1);


        current_statement_begin__ = 49;
        for (int i = 1; i <= size(p); ++i) {
            {
            current_statement_begin__ = 50;
            int nc(0);
            (void) nc;  // dummy to suppress unused var warning

            stan::math::fill(nc, std::numeric_limits<int>::min());
            stan::math::assign(nc,get_base1(p,i,"p",1));


            current_statement_begin__ = 51;
            if (as_bool(logical_eq(nc,1))) {

                current_statement_begin__ = 52;
                stan::math::assign(get_base1_lhs(theta_L,theta_L_mark,"theta_L",1), ((get_base1(tau,i,"tau",1) * get_base1(scale,i,"scale",1)) * dispersion));
                current_statement_begin__ = 54;
                stan::math::assign(theta_L_mark, (theta_L_mark + 1));
            } else {
                {
                current_statement_begin__ = 57;
                validate_non_negative_index("T_i", "nc", nc);
                validate_non_negative_index("T_i", "nc", nc);
                Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,Eigen::Dynamic>  T_i(static_cast<Eigen::VectorXd::Index>(nc),static_cast<Eigen::VectorXd::Index>(nc));
                (void) T_i;  // dummy to suppress unused var warning

                stan::math::initialize(T_i, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(T_i,DUMMY_VAR__);
                current_statement_begin__ = 58;
                fun_scalar_t__ std_dev;
                (void) std_dev;  // dummy to suppress unused var warning

                stan::math::initialize(std_dev, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(std_dev,DUMMY_VAR__);
                current_statement_begin__ = 59;
                fun_scalar_t__ T21;
                (void) T21;  // dummy to suppress unused var warning

                stan::math::initialize(T21, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(T21,DUMMY_VAR__);
                current_statement_begin__ = 60;
                fun_scalar_t__ trace_T_i;
                (void) trace_T_i;  // dummy to suppress unused var warning

                stan::math::initialize(trace_T_i, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(trace_T_i,DUMMY_VAR__);
                stan::math::assign(trace_T_i,(square(((get_base1(tau,i,"tau",1) * get_base1(scale,i,"scale",1)) * dispersion)) * nc));
                current_statement_begin__ = 61;
                validate_non_negative_index("pi", "nc", nc);
                Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  pi(static_cast<Eigen::VectorXd::Index>(nc));
                (void) pi;  // dummy to suppress unused var warning

                stan::math::initialize(pi, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(pi,DUMMY_VAR__);
                stan::math::assign(pi,segment(zeta,zeta_mark,nc));


                current_statement_begin__ = 62;
                stan::math::assign(pi, divide(pi,sum(pi)));
                current_statement_begin__ = 65;
                stan::math::assign(zeta_mark, (zeta_mark + nc));
                current_statement_begin__ = 66;
                stan::math::assign(std_dev, sqrt((get_base1(pi,1,"pi",1) * trace_T_i)));
                current_statement_begin__ = 67;
                stan::math::assign(get_base1_lhs(T_i,1,1,"T_i",1), std_dev);
                current_statement_begin__ = 70;
                stan::math::assign(std_dev, sqrt((get_base1(pi,2,"pi",1) * trace_T_i)));
                current_statement_begin__ = 71;
                stan::math::assign(T21, ((2.0 * get_base1(rho,rho_mark,"rho",1)) - 1.0));
                current_statement_begin__ = 72;
                stan::math::assign(rho_mark, (rho_mark + 1));
                current_statement_begin__ = 73;
                stan::math::assign(get_base1_lhs(T_i,2,2,"T_i",1), (std_dev * sqrt((1.0 - square(T21)))));
                current_statement_begin__ = 74;
                stan::math::assign(get_base1_lhs(T_i,2,1,"T_i",1), (std_dev * T21));
                current_statement_begin__ = 76;
                for (int r = 2; r <= (nc - 1); ++r) {
                    {
                    current_statement_begin__ = 77;
                    int rp1(0);
                    (void) rp1;  // dummy to suppress unused var warning

                    stan::math::fill(rp1, std::numeric_limits<int>::min());
                    stan::math::assign(rp1,(r + 1));
                    current_statement_begin__ = 78;
                    validate_non_negative_index("T_row", "r", r);
                    Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  T_row(static_cast<Eigen::VectorXd::Index>(r));
                    (void) T_row;  // dummy to suppress unused var warning

                    stan::math::initialize(T_row, std::numeric_limits<double>::quiet_NaN());
                    stan::math::fill(T_row,DUMMY_VAR__);
                    stan::math::assign(T_row,segment(z_T,z_T_mark,r));
                    current_statement_begin__ = 79;
                    fun_scalar_t__ scale_factor;
                    (void) scale_factor;  // dummy to suppress unused var warning

                    stan::math::initialize(scale_factor, std::numeric_limits<double>::quiet_NaN());
                    stan::math::fill(scale_factor,DUMMY_VAR__);
                    stan::math::assign(scale_factor,(sqrt((get_base1(rho,rho_mark,"rho",1) / dot_self(T_row))) * std_dev));


                    current_statement_begin__ = 80;
                    stan::math::assign(z_T_mark, (z_T_mark + r));
                    current_statement_begin__ = 81;
                    stan::math::assign(std_dev, sqrt((get_base1(pi,rp1,"pi",1) * trace_T_i)));
                    current_statement_begin__ = 82;
                    for (int c = 1; c <= r; ++c) {
                        current_statement_begin__ = 82;
                        stan::math::assign(get_base1_lhs(T_i,rp1,c,"T_i",1), (get_base1(T_row,c,"T_row",1) * scale_factor));
                    }
                    current_statement_begin__ = 83;
                    stan::math::assign(get_base1_lhs(T_i,rp1,rp1,"T_i",1), (sqrt((1.0 - get_base1(rho,rho_mark,"rho",1))) * std_dev));
                    current_statement_begin__ = 84;
                    stan::math::assign(rho_mark, (rho_mark + 1));
                    }
                }
                current_statement_begin__ = 88;
                for (int c = 1; c <= nc; ++c) {
                    current_statement_begin__ = 88;
                    for (int r = c; r <= nc; ++r) {

                        current_statement_begin__ = 89;
                        stan::math::assign(get_base1_lhs(theta_L,theta_L_mark,"theta_L",1), get_base1(T_i,r,c,"T_i",1));
                        current_statement_begin__ = 90;
                        stan::math::assign(theta_L_mark, (theta_L_mark + 1));
                    }
                }
                }
            }
            }
        }
        current_statement_begin__ = 94;
        return stan::math::promote_scalar<fun_return_scalar_t__>(theta_L);
        }
    } 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 ***");
    }
}


struct make_theta_L_functor__ {
    template <typename T2__, typename T3__, typename T4__, typename T5__, typename T6__, typename T7__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T2__, T3__, T4__, T5__, typename boost::math::tools::promote_args<T6__, T7__>::type>::type, Eigen::Dynamic,1>
    operator()(const int& len_theta_L,
                 const std::vector<int>& p,
                 const T2__& dispersion,
                 const Eigen::Matrix<T3__, Eigen::Dynamic,1>& tau,
                 const Eigen::Matrix<T4__, Eigen::Dynamic,1>& scale,
                 const Eigen::Matrix<T5__, Eigen::Dynamic,1>& zeta,
                 const Eigen::Matrix<T6__, Eigen::Dynamic,1>& rho,
                 const Eigen::Matrix<T7__, Eigen::Dynamic,1>& z_T, std::ostream* pstream__) const {
        return make_theta_L(len_theta_L, p, dispersion, tau, scale, zeta, rho, z_T, pstream__);
    }
};

template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic,1>
make_b(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_b,
           const Eigen::Matrix<T1__, Eigen::Dynamic,1>& theta_L,
           const std::vector<int>& p,
           const std::vector<int>& l, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 109;
        validate_non_negative_index("b", "rows(z_b)", rows(z_b));
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  b(static_cast<Eigen::VectorXd::Index>(rows(z_b)));
        (void) b;  // dummy to suppress unused var warning

        stan::math::initialize(b, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(b,DUMMY_VAR__);
        current_statement_begin__ = 110;
        int b_mark(0);
        (void) b_mark;  // dummy to suppress unused var warning

        stan::math::fill(b_mark, std::numeric_limits<int>::min());
        stan::math::assign(b_mark,1);
        current_statement_begin__ = 111;
        int theta_L_mark(0);
        (void) theta_L_mark;  // dummy to suppress unused var warning

        stan::math::fill(theta_L_mark, std::numeric_limits<int>::min());
        stan::math::assign(theta_L_mark,1);


        current_statement_begin__ = 112;
        for (int i = 1; i <= size(p); ++i) {
            {
            current_statement_begin__ = 113;
            int nc(0);
            (void) nc;  // dummy to suppress unused var warning

            stan::math::fill(nc, std::numeric_limits<int>::min());
            stan::math::assign(nc,get_base1(p,i,"p",1));


            current_statement_begin__ = 114;
            if (as_bool(logical_eq(nc,1))) {
                {
                current_statement_begin__ = 115;
                fun_scalar_t__ theta_L_start;
                (void) theta_L_start;  // dummy to suppress unused var warning

                stan::math::initialize(theta_L_start, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(theta_L_start,DUMMY_VAR__);
                stan::math::assign(theta_L_start,get_base1(theta_L,theta_L_mark,"theta_L",1));


                current_statement_begin__ = 116;
                for (int s = b_mark; s <= ((b_mark + get_base1(l,i,"l",1)) - 1); ++s) {
                    current_statement_begin__ = 117;
                    stan::math::assign(get_base1_lhs(b,s,"b",1), (theta_L_start * get_base1(z_b,s,"z_b",1)));
                }
                current_statement_begin__ = 118;
                stan::math::assign(b_mark, (b_mark + get_base1(l,i,"l",1)));
                current_statement_begin__ = 119;
                stan::math::assign(theta_L_mark, (theta_L_mark + 1));
                }
            } else {
                {
                current_statement_begin__ = 122;
                validate_non_negative_index("T_i", "nc", nc);
                validate_non_negative_index("T_i", "nc", nc);
                Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,Eigen::Dynamic>  T_i(static_cast<Eigen::VectorXd::Index>(nc),static_cast<Eigen::VectorXd::Index>(nc));
                (void) T_i;  // dummy to suppress unused var warning

                stan::math::initialize(T_i, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(T_i,DUMMY_VAR__);
                stan::math::assign(T_i,rep_matrix(0,nc,nc));


                current_statement_begin__ = 123;
                for (int c = 1; c <= nc; ++c) {

                    current_statement_begin__ = 124;
                    stan::math::assign(get_base1_lhs(T_i,c,c,"T_i",1), get_base1(theta_L,theta_L_mark,"theta_L",1));
                    current_statement_begin__ = 125;
                    stan::math::assign(theta_L_mark, (theta_L_mark + 1));
                    current_statement_begin__ = 126;
                    for (int r = (c + 1); r <= nc; ++r) {

                        current_statement_begin__ = 127;
                        stan::math::assign(get_base1_lhs(T_i,r,c,"T_i",1), get_base1(theta_L,theta_L_mark,"theta_L",1));
                        current_statement_begin__ = 128;
                        stan::math::assign(theta_L_mark, (theta_L_mark + 1));
                    }
                }
                current_statement_begin__ = 131;
                for (int j = 1; j <= get_base1(l,i,"l",1); ++j) {
                    {
                    current_statement_begin__ = 132;
                    validate_non_negative_index("temp", "nc", nc);
                    Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  temp(static_cast<Eigen::VectorXd::Index>(nc));
                    (void) temp;  // dummy to suppress unused var warning

                    stan::math::initialize(temp, std::numeric_limits<double>::quiet_NaN());
                    stan::math::fill(temp,DUMMY_VAR__);
                    stan::math::assign(temp,multiply(T_i,segment(z_b,b_mark,nc)));


                    current_statement_begin__ = 133;
                    stan::math::assign(b_mark, (b_mark - 1));
                    current_statement_begin__ = 134;
                    for (int s = 1; s <= nc; ++s) {
                        current_statement_begin__ = 134;
                        stan::math::assign(get_base1_lhs(b,(b_mark + s),"b",1), get_base1(temp,s,"temp",1));
                    }
                    current_statement_begin__ = 135;
                    stan::math::assign(b_mark, ((b_mark + nc) + 1));
                    }
                }
                }
            }
            }
        }
        current_statement_begin__ = 139;
        return stan::math::promote_scalar<fun_return_scalar_t__>(b);
        }
    } 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 ***");
    }
}


struct make_b_functor__ {
    template <typename T0__, typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic,1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_b,
           const Eigen::Matrix<T1__, Eigen::Dynamic,1>& theta_L,
           const std::vector<int>& p,
           const std::vector<int>& l, std::ostream* pstream__) const {
        return make_b(z_b, theta_L, p, l, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T6__, typename T7__, typename T_lp__, typename T_lp_accum__>
void
decov_lp(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_b,
             const Eigen::Matrix<T1__, Eigen::Dynamic,1>& z_T,
             const Eigen::Matrix<T2__, Eigen::Dynamic,1>& rho,
             const Eigen::Matrix<T3__, Eigen::Dynamic,1>& zeta,
             const Eigen::Matrix<T4__, Eigen::Dynamic,1>& tau,
             const std::vector<T5__>& regularization,
             const std::vector<T6__>& delta,
             const Eigen::Matrix<T7__, Eigen::Dynamic,1>& shape,
             const int& t,
             const std::vector<int>& p, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T6__, T7__, T_lp__>::type>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 160;
        int pos_reg(0);
        (void) pos_reg;  // dummy to suppress unused var warning

        stan::math::fill(pos_reg, std::numeric_limits<int>::min());
        stan::math::assign(pos_reg,1);
        current_statement_begin__ = 161;
        int pos_rho(0);
        (void) pos_rho;  // dummy to suppress unused var warning

        stan::math::fill(pos_rho, std::numeric_limits<int>::min());
        stan::math::assign(pos_rho,1);


        current_statement_begin__ = 162;
        lp_accum__.add(normal_log(z_b,0,1));
        current_statement_begin__ = 163;
        lp_accum__.add(normal_log(z_T,0,1));
        current_statement_begin__ = 164;
        for (int i = 1; i <= t; ++i) {
            current_statement_begin__ = 164;
            if (as_bool(logical_gt(get_base1(p,i,"p",1),1))) {
                {
                current_statement_begin__ = 165;
                validate_non_negative_index("shape1", "(get_base1(p,i,\"p\",1) - 1)", (get_base1(p,i,"p",1) - 1));
                Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  shape1(static_cast<Eigen::VectorXd::Index>((get_base1(p,i,"p",1) - 1)));
                (void) shape1;  // dummy to suppress unused var warning

                stan::math::initialize(shape1, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(shape1,DUMMY_VAR__);
                current_statement_begin__ = 166;
                validate_non_negative_index("shape2", "(get_base1(p,i,\"p\",1) - 1)", (get_base1(p,i,"p",1) - 1));
                Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  shape2(static_cast<Eigen::VectorXd::Index>((get_base1(p,i,"p",1) - 1)));
                (void) shape2;  // dummy to suppress unused var warning

                stan::math::initialize(shape2, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(shape2,DUMMY_VAR__);
                current_statement_begin__ = 167;
                fun_scalar_t__ nu;
                (void) nu;  // dummy to suppress unused var warning

                stan::math::initialize(nu, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(nu,DUMMY_VAR__);
                stan::math::assign(nu,(get_base1(regularization,pos_reg,"regularization",1) + (0.5 * (get_base1(p,i,"p",1) - 2))));


                current_statement_begin__ = 168;
                stan::math::assign(pos_reg, (pos_reg + 1));
                current_statement_begin__ = 169;
                stan::math::assign(get_base1_lhs(shape1,1,"shape1",1), nu);
                current_statement_begin__ = 170;
                stan::math::assign(get_base1_lhs(shape2,1,"shape2",1), nu);
                current_statement_begin__ = 171;
                for (int j = 2; j <= (get_base1(p,i,"p",1) - 1); ++j) {

                    current_statement_begin__ = 172;
                    stan::math::assign(nu, (nu - 0.5));
                    current_statement_begin__ = 173;
                    stan::math::assign(get_base1_lhs(shape1,j,"shape1",1), (0.5 * j));
                    current_statement_begin__ = 174;
                    stan::math::assign(get_base1_lhs(shape2,j,"shape2",1), nu);
                }
                current_statement_begin__ = 176;
                lp_accum__.add(beta_log(stan::model::rvalue(rho, stan::model::cons_list(stan::model::index_min_max(pos_rho, ((pos_rho + get_base1(p,i,"p",1)) - 2)), stan::model::nil_index_list()), "rho"),shape1,shape2));
                current_statement_begin__ = 177;
                stan::math::assign(pos_rho, ((pos_rho + get_base1(p,i,"p",1)) - 1));
                }
            }
        }
        current_statement_begin__ = 179;
        lp_accum__.add(gamma_log(zeta,delta,1));
        current_statement_begin__ = 180;
        lp_accum__.add(gamma_log(tau,shape,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 ***");
    }
}


struct decov_lp_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T6__, typename T7__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_b,
             const Eigen::Matrix<T1__, Eigen::Dynamic,1>& z_T,
             const Eigen::Matrix<T2__, Eigen::Dynamic,1>& rho,
             const Eigen::Matrix<T3__, Eigen::Dynamic,1>& zeta,
             const Eigen::Matrix<T4__, Eigen::Dynamic,1>& tau,
             const std::vector<T5__>& regularization,
             const std::vector<T6__>& delta,
             const Eigen::Matrix<T7__, Eigen::Dynamic,1>& shape,
             const int& t,
             const std::vector<int>& p, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return decov_lp(z_b, z_T, rho, zeta, tau, regularization, delta, shape, t, p, lp__, lp_accum__, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type, Eigen::Dynamic,1>
hs_prior(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_beta,
             const std::vector<T1__>& global,
             const std::vector<Eigen::Matrix<T2__, Eigen::Dynamic,1> >& local,
             const T3__& global_prior_scale,
             const T4__& error_scale,
             const T5__& c2, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 196;
        int K(0);
        (void) K;  // dummy to suppress unused var warning

        stan::math::fill(K, std::numeric_limits<int>::min());
        stan::math::assign(K,rows(z_beta));
        current_statement_begin__ = 197;
        validate_non_negative_index("lambda", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  lambda(static_cast<Eigen::VectorXd::Index>(K));
        (void) lambda;  // dummy to suppress unused var warning

        stan::math::initialize(lambda, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(lambda,DUMMY_VAR__);
        stan::math::assign(lambda,elt_multiply(get_base1(local,1,"local",1),sqrt(get_base1(local,2,"local",1))));
        current_statement_begin__ = 198;
        fun_scalar_t__ tau;
        (void) tau;  // dummy to suppress unused var warning

        stan::math::initialize(tau, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(tau,DUMMY_VAR__);
        stan::math::assign(tau,(((get_base1(global,1,"global",1) * sqrt(get_base1(global,2,"global",1))) * global_prior_scale) * error_scale));
        current_statement_begin__ = 199;
        validate_non_negative_index("lambda2", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  lambda2(static_cast<Eigen::VectorXd::Index>(K));
        (void) lambda2;  // dummy to suppress unused var warning

        stan::math::initialize(lambda2, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(lambda2,DUMMY_VAR__);
        stan::math::assign(lambda2,square(lambda));
        current_statement_begin__ = 200;
        validate_non_negative_index("lambda_tilde", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  lambda_tilde(static_cast<Eigen::VectorXd::Index>(K));
        (void) lambda_tilde;  // dummy to suppress unused var warning

        stan::math::initialize(lambda_tilde, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(lambda_tilde,DUMMY_VAR__);
        stan::math::assign(lambda_tilde,sqrt(elt_divide(multiply(c2,lambda2),add(c2,multiply(square(tau),lambda2)))));


        current_statement_begin__ = 201;
        return stan::math::promote_scalar<fun_return_scalar_t__>(multiply(elt_multiply(z_beta,lambda_tilde),tau));
        }
    } 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 ***");
    }
}


struct hs_prior_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type, Eigen::Dynamic,1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_beta,
             const std::vector<T1__>& global,
             const std::vector<Eigen::Matrix<T2__, Eigen::Dynamic,1> >& local,
             const T3__& global_prior_scale,
             const T4__& error_scale,
             const T5__& c2, std::ostream* pstream__) const {
        return hs_prior(z_beta, global, local, global_prior_scale, error_scale, c2, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type, Eigen::Dynamic,1>
hsplus_prior(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_beta,
                 const std::vector<T1__>& global,
                 const std::vector<Eigen::Matrix<T2__, Eigen::Dynamic,1> >& local,
                 const T3__& global_prior_scale,
                 const T4__& error_scale,
                 const T5__& c2, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 217;
        int K(0);
        (void) K;  // dummy to suppress unused var warning

        stan::math::fill(K, std::numeric_limits<int>::min());
        stan::math::assign(K,rows(z_beta));
        current_statement_begin__ = 218;
        validate_non_negative_index("lambda", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  lambda(static_cast<Eigen::VectorXd::Index>(K));
        (void) lambda;  // dummy to suppress unused var warning

        stan::math::initialize(lambda, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(lambda,DUMMY_VAR__);
        stan::math::assign(lambda,elt_multiply(get_base1(local,1,"local",1),sqrt(get_base1(local,2,"local",1))));
        current_statement_begin__ = 219;
        validate_non_negative_index("eta", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  eta(static_cast<Eigen::VectorXd::Index>(K));
        (void) eta;  // dummy to suppress unused var warning

        stan::math::initialize(eta, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(eta,DUMMY_VAR__);
        stan::math::assign(eta,elt_multiply(get_base1(local,3,"local",1),sqrt(get_base1(local,4,"local",1))));
        current_statement_begin__ = 220;
        fun_scalar_t__ tau;
        (void) tau;  // dummy to suppress unused var warning

        stan::math::initialize(tau, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(tau,DUMMY_VAR__);
        stan::math::assign(tau,(((get_base1(global,1,"global",1) * sqrt(get_base1(global,2,"global",1))) * global_prior_scale) * error_scale));
        current_statement_begin__ = 221;
        validate_non_negative_index("lambda_eta2", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  lambda_eta2(static_cast<Eigen::VectorXd::Index>(K));
        (void) lambda_eta2;  // dummy to suppress unused var warning

        stan::math::initialize(lambda_eta2, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(lambda_eta2,DUMMY_VAR__);
        stan::math::assign(lambda_eta2,square(elt_multiply(lambda,eta)));
        current_statement_begin__ = 222;
        validate_non_negative_index("lambda_tilde", "K", K);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  lambda_tilde(static_cast<Eigen::VectorXd::Index>(K));
        (void) lambda_tilde;  // dummy to suppress unused var warning

        stan::math::initialize(lambda_tilde, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(lambda_tilde,DUMMY_VAR__);
        stan::math::assign(lambda_tilde,sqrt(elt_divide(multiply(c2,lambda_eta2),add(c2,multiply(square(tau),lambda_eta2)))));


        current_statement_begin__ = 224;
        return stan::math::promote_scalar<fun_return_scalar_t__>(multiply(elt_multiply(z_beta,lambda_tilde),tau));
        }
    } 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 ***");
    }
}


struct hsplus_prior_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type, Eigen::Dynamic,1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& z_beta,
                 const std::vector<T1__>& global,
                 const std::vector<Eigen::Matrix<T2__, Eigen::Dynamic,1> >& local,
                 const T3__& global_prior_scale,
                 const T4__& error_scale,
                 const T5__& c2, std::ostream* pstream__) const {
        return hsplus_prior(z_beta, global, local, global_prior_scale, error_scale, c2, pstream__);
    }
};

template <typename T0__, typename T1__>
typename boost::math::tools::promote_args<T0__, T1__>::type
CFt(const T0__& z,
        const T1__& df, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 238;
        fun_scalar_t__ z2;
        (void) z2;  // dummy to suppress unused var warning

        stan::math::initialize(z2, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(z2,DUMMY_VAR__);
        stan::math::assign(z2,square(z));
        current_statement_begin__ = 239;
        fun_scalar_t__ z3;
        (void) z3;  // dummy to suppress unused var warning

        stan::math::initialize(z3, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(z3,DUMMY_VAR__);
        stan::math::assign(z3,(z2 * z));
        current_statement_begin__ = 240;
        fun_scalar_t__ z5;
        (void) z5;  // dummy to suppress unused var warning

        stan::math::initialize(z5, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(z5,DUMMY_VAR__);
        stan::math::assign(z5,(z2 * z3));
        current_statement_begin__ = 241;
        fun_scalar_t__ z7;
        (void) z7;  // dummy to suppress unused var warning

        stan::math::initialize(z7, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(z7,DUMMY_VAR__);
        stan::math::assign(z7,(z2 * z5));
        current_statement_begin__ = 242;
        fun_scalar_t__ z9;
        (void) z9;  // dummy to suppress unused var warning

        stan::math::initialize(z9, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(z9,DUMMY_VAR__);
        stan::math::assign(z9,(z2 * z7));
        current_statement_begin__ = 243;
        fun_scalar_t__ df2;
        (void) df2;  // dummy to suppress unused var warning

        stan::math::initialize(df2, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(df2,DUMMY_VAR__);
        stan::math::assign(df2,square(df));
        current_statement_begin__ = 244;
        fun_scalar_t__ df3;
        (void) df3;  // dummy to suppress unused var warning

        stan::math::initialize(df3, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(df3,DUMMY_VAR__);
        stan::math::assign(df3,(df2 * df));
        current_statement_begin__ = 245;
        fun_scalar_t__ df4;
        (void) df4;  // dummy to suppress unused var warning

        stan::math::initialize(df4, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(df4,DUMMY_VAR__);
        stan::math::assign(df4,(df2 * df2));


        current_statement_begin__ = 246;
        return stan::math::promote_scalar<fun_return_scalar_t__>(((((z + ((z3 + z) / (4 * df))) + ((((5 * z5) + (16 * z3)) + (3 * z)) / (96 * df2))) + (((((3 * z7) + (19 * z5)) + (17 * z3)) - (15 * z)) / (384 * df3))) + ((((((79 * z9) + (776 * z7)) + (1482 * z5)) - (1920 * z3)) - (945 * z)) / (92160 * df4))));
        }
    } 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 ***");
    }
}


struct CFt_functor__ {
    template <typename T0__, typename T1__>
        typename boost::math::tools::promote_args<T0__, T1__>::type
    operator()(const T0__& z,
        const T1__& df, std::ostream* pstream__) const {
        return CFt(z, df, pstream__);
    }
};

std::vector<std::vector<int> >
make_V(const int& N,
           const int& t,
           const std::vector<int>& v, std::ostream* pstream__) {
    typedef double fun_scalar_t__;
    typedef int fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 260;
        validate_non_negative_index("V", "t", t);
        validate_non_negative_index("V", "N", N);
        vector<vector<int> > V(t, (vector<int>(N, 0)));
        stan::math::fill(V, std::numeric_limits<int>::min());
        current_statement_begin__ = 261;
        int pos(0);
        (void) pos;  // dummy to suppress unused var warning

        stan::math::fill(pos, std::numeric_limits<int>::min());
        stan::math::assign(pos,1);


        current_statement_begin__ = 262;
        if (as_bool(logical_gt(t,0))) {
            current_statement_begin__ = 262;
            for (int j = 1; j <= N; ++j) {
                current_statement_begin__ = 262;
                for (int i = 1; i <= t; ++i) {

                    current_statement_begin__ = 263;
                    stan::math::assign(get_base1_lhs(get_base1_lhs(V,i,"V",1),j,"V",2), (get_base1(v,pos,"v",1) + 1));
                    current_statement_begin__ = 264;
                    stan::math::assign(pos, (pos + 1));
                }
            }
        }
        current_statement_begin__ = 266;
        return stan::math::promote_scalar<fun_return_scalar_t__>(V);
        }
    } 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 ***");
    }
}


struct make_V_functor__ {
            std::vector<std::vector<int> >
    operator()(const int& N,
           const int& t,
           const std::vector<int>& v, std::ostream* pstream__) const {
        return make_V(N, t, v, pstream__);
    }
};

template <typename T2__, typename T5__>
Eigen::Matrix<typename boost::math::tools::promote_args<T2__, T5__>::type, Eigen::Dynamic,1>
csr_matrix_times_vector2(const int& m,
                             const int& n,
                             const Eigen::Matrix<T2__, Eigen::Dynamic,1>& w,
                             const std::vector<int>& v,
                             const std::vector<int>& u,
                             const Eigen::Matrix<T5__, Eigen::Dynamic,1>& b, std::ostream* pstream__);

double
make_lower(const int& family,
               const int& link, std::ostream* pstream__) {
    typedef double fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {

        current_statement_begin__ = 300;
        if (as_bool(logical_eq(family,1))) {
            current_statement_begin__ = 300;
            return stan::math::promote_scalar<fun_return_scalar_t__>(stan::math::negative_infinity());
        }
        current_statement_begin__ = 301;
        if (as_bool(logical_lte(family,3))) {

            current_statement_begin__ = 302;
            if (as_bool(logical_eq(link,2))) {
                current_statement_begin__ = 302;
                return stan::math::promote_scalar<fun_return_scalar_t__>(stan::math::negative_infinity());
            }
            current_statement_begin__ = 303;
            return stan::math::promote_scalar<fun_return_scalar_t__>(0);
        }
        current_statement_begin__ = 305;
        return stan::math::promote_scalar<fun_return_scalar_t__>(stan::math::negative_infinity());
    } 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 ***");
    }
}


struct make_lower_functor__ {
            double
    operator()(const int& family,
               const int& link, std::ostream* pstream__) const {
        return make_lower(family, link, pstream__);
    }
};

double
make_upper(const int& family,
               const int& link, std::ostream* pstream__) {
    typedef double fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {

        current_statement_begin__ = 316;
        if (as_bool((primitive_value(logical_eq(family,4)) && primitive_value(logical_eq(link,5))))) {
            current_statement_begin__ = 316;
            return stan::math::promote_scalar<fun_return_scalar_t__>(0);
        }
        current_statement_begin__ = 317;
        return stan::math::promote_scalar<fun_return_scalar_t__>(stan::math::positive_infinity());
    } 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 ***");
    }
}


struct make_upper_functor__ {
            double
    operator()(const int& family,
               const int& link, std::ostream* pstream__) const {
        return make_upper(family, link, pstream__);
    }
};

template <typename T0__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__>::type, Eigen::Dynamic,1>
linkinv_count(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& eta,
                  const int& link, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {

        current_statement_begin__ = 328;
        if (as_bool(logical_eq(link,1))) {
            current_statement_begin__ = 328;
            return stan::math::promote_scalar<fun_return_scalar_t__>(exp(eta));
        } else if (as_bool(logical_eq(link,2))) {
            current_statement_begin__ = 329;
            return stan::math::promote_scalar<fun_return_scalar_t__>(eta);
        } else if (as_bool(logical_eq(link,3))) {
            current_statement_begin__ = 330;
            return stan::math::promote_scalar<fun_return_scalar_t__>(square(eta));
        } else {
            current_statement_begin__ = 331;
            std::stringstream errmsg_stream__;
            errmsg_stream__ << "Invalid link";
            throw std::domain_error(errmsg_stream__.str());
        }
        current_statement_begin__ = 332;
        return stan::math::promote_scalar<fun_return_scalar_t__>(eta);
    } 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 ***");
    }
}


struct linkinv_count_functor__ {
    template <typename T0__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__>::type, Eigen::Dynamic,1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic,1>& eta,
                  const int& link, std::ostream* pstream__) const {
        return linkinv_count(eta, link, pstream__);
    }
};

template <typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__>::type, Eigen::Dynamic,1>
pw_pois(const std::vector<int>& y,
            const Eigen::Matrix<T1__, Eigen::Dynamic,1>& eta,
            const int& link, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 344;
        int N(0);
        (void) N;  // dummy to suppress unused var warning

        stan::math::fill(N, std::numeric_limits<int>::min());
        stan::math::assign(N,rows(eta));
        current_statement_begin__ = 345;
        validate_non_negative_index("ll", "N", N);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  ll(static_cast<Eigen::VectorXd::Index>(N));
        (void) ll;  // dummy to suppress unused var warning

        stan::math::initialize(ll, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(ll,DUMMY_VAR__);


        current_statement_begin__ = 346;
        if (as_bool(logical_eq(link,1))) {
            current_statement_begin__ = 347;
            for (int n = 1; n <= N; ++n) {
                current_statement_begin__ = 347;
                stan::math::assign(get_base1_lhs(ll,n,"ll",1), poisson_log_log(get_base1(y,n,"y",1),get_base1(eta,n,"eta",1)));
            }
        } else if (as_bool(logical_lte(link,3))) {
            {
            current_statement_begin__ = 349;
            validate_non_negative_index("phi", "N", N);
            Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  phi(static_cast<Eigen::VectorXd::Index>(N));
            (void) phi;  // dummy to suppress unused var warning

            stan::math::initialize(phi, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(phi,DUMMY_VAR__);
            stan::math::assign(phi,linkinv_count(eta,link, pstream__));


            current_statement_begin__ = 350;
            for (int n = 1; n <= N; ++n) {
                current_statement_begin__ = 350;
                stan::math::assign(get_base1_lhs(ll,n,"ll",1), poisson_log(get_base1(y,n,"y",1),get_base1(phi,n,"phi",1)));
            }
            }
        } else {
            current_statement_begin__ = 352;
            std::stringstream errmsg_stream__;
            errmsg_stream__ << "Invalid link";
            throw std::domain_error(errmsg_stream__.str());
        }
        current_statement_begin__ = 353;
        return stan::math::promote_scalar<fun_return_scalar_t__>(ll);
        }
    } 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 ***");
    }
}


struct pw_pois_functor__ {
    template <typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__>::type, Eigen::Dynamic,1>
    operator()(const std::vector<int>& y,
            const Eigen::Matrix<T1__, Eigen::Dynamic,1>& eta,
            const int& link, std::ostream* pstream__) const {
        return pw_pois(y, eta, link, pstream__);
    }
};

template <typename T1__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__>::type, Eigen::Dynamic,1>
pw_nb(const std::vector<int>& y,
          const Eigen::Matrix<T1__, Eigen::Dynamic,1>& eta,
          const T2__& theta,
          const int& link, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__, T2__>::type fun_scalar_t__;
    typedef fun_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        fun_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 366;
        int N(0);
        (void) N;  // dummy to suppress unused var warning

        stan::math::fill(N, std::numeric_limits<int>::min());
        stan::math::assign(N,rows(eta));
        current_statement_begin__ = 367;
        validate_non_negative_index("rho", "N", N);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  rho(static_cast<Eigen::VectorXd::Index>(N));
        (void) rho;  // dummy to suppress unused var warning

        stan::math::initialize(rho, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(rho,DUMMY_VAR__);
        stan::math::assign(rho,linkinv_count(eta,link, pstream__));
        current_statement_begin__ = 368;
        validate_non_negative_index("ll", "N", N);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,1>  ll(static_cast<Eigen::VectorXd::Index>(N));
        (void) ll;  // dummy to suppress unused var warning

        stan::math::initialize(ll, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(ll,DUMMY_VAR__);


        current_statement_begin__ = 369;
        for (int n = 1; n <= N; ++n) {
            current_statement_begin__ = 369;
            stan::math::assign(get_base1_lhs(ll,n,"ll",1), neg_binomial_2_log(get_base1(y,n,"y",1),get_base1(rho,n,"rho",1),theta));
        }
        current_statement_begin__ = 370;
        return stan::math::promote_scalar<fun_return_scalar_t__>(ll);
        }
    } 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 ***");
    }
}


struct pw_nb_functor__ {
    template <typename T1__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__>::type, Eigen::Dynamic,1>
    operator()(const std::vector<int>& y,
          const Eigen::Matrix<T1__, Eigen::Dynamic,1>& eta,
          const T2__& theta,
          const int& link, std::ostream* pstream__) const {
        return pw_nb(y, eta, theta, link, pstream__);
    }
};

class model_count : public prob_grad {
private:
    int N;
    int K;
    vector_d xbar;
    int dense_X;
    vector<matrix_d> X;
    int nnz_X;
    vector_d w_X;
    vector<int> v_X;
    vector<int> u_X;
    int K_smooth;
    matrix_d S;
    vector<int> smooth_map;
    vector<int> y;
    int prior_PD;
    int has_intercept;
    int link;
    int prior_dist;
    int prior_dist_for_intercept;
    int prior_dist_for_aux;
    int prior_dist_for_smooth;
    int has_weights;
    vector_d weights;
    int has_offset;
    vector_d offset;
    int family;
    vector_d prior_scale;
    double prior_scale_for_intercept;
    double prior_scale_for_aux;
    vector_d prior_scale_for_smooth;
    vector_d prior_mean;
    double prior_mean_for_intercept;
    double prior_mean_for_aux;
    vector_d prior_mean_for_smooth;
    vector_d prior_df;
    double prior_df_for_intercept;
    double prior_df_for_aux;
    vector_d prior_df_for_smooth;
    double global_prior_df;
    double global_prior_scale;
    double slab_df;
    double slab_scale;
    vector<int> num_normals;
    int t;
    vector<int> p;
    vector<int> l;
    int q;
    int len_theta_L;
    vector_d shape;
    vector_d scale;
    int len_concentration;
    vector<double> concentration;
    int len_regularization;
    vector<double> regularization;
    int num_non_zero;
    vector_d w;
    vector<int> v;
    vector<int> u;
    int special_case;
    double poisson_max;
    vector<vector<int> > V;
    int len_z_T;
    int len_var_group;
    int len_rho;
    int is_continuous;
    int pos;
    vector<double> delta;
    int hs;
public:
    model_count(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, 0, pstream__);
    }

    model_count(stan::io::var_context& context__,
        unsigned int random_seed__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, random_seed__, pstream__);
    }

    void ctor_body(stan::io::var_context& context__,
                   unsigned int random_seed__,
                   std::ostream* pstream__) {
        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_count_namespace::model_count";
        (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__;
        double DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

        // initialize member variables
        try {
            current_statement_begin__ = 376;
            context__.validate_dims("data initialization", "N", "int", context__.to_vec());
            N = int(0);
            vals_i__ = context__.vals_i("N");
            pos__ = 0;
            N = vals_i__[pos__++];
            current_statement_begin__ = 377;
            context__.validate_dims("data initialization", "K", "int", context__.to_vec());
            K = int(0);
            vals_i__ = context__.vals_i("K");
            pos__ = 0;
            K = vals_i__[pos__++];
            current_statement_begin__ = 380;
            validate_non_negative_index("xbar", "K", K);
            context__.validate_dims("data initialization", "xbar", "vector_d", context__.to_vec(K));
            validate_non_negative_index("xbar", "K", K);
            xbar = vector_d(static_cast<Eigen::VectorXd::Index>(K));
            vals_r__ = context__.vals_r("xbar");
            pos__ = 0;
            size_t xbar_i_vec_lim__ = K;
            for (size_t i_vec__ = 0; i_vec__ < xbar_i_vec_lim__; ++i_vec__) {
                xbar[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 381;
            context__.validate_dims("data initialization", "dense_X", "int", context__.to_vec());
            dense_X = int(0);
            vals_i__ = context__.vals_i("dense_X");
            pos__ = 0;
            dense_X = vals_i__[pos__++];
            current_statement_begin__ = 382;
            validate_non_negative_index("X", "dense_X", dense_X);
            validate_non_negative_index("X", "N", N);
            validate_non_negative_index("X", "K", K);
            context__.validate_dims("data initialization", "X", "matrix_d", context__.to_vec(dense_X,N,K));
            validate_non_negative_index("X", "dense_X", dense_X);
            validate_non_negative_index("X", "N", N);
            validate_non_negative_index("X", "K", K);
            X = std::vector<matrix_d>(dense_X,matrix_d(static_cast<Eigen::VectorXd::Index>(N),static_cast<Eigen::VectorXd::Index>(K)));
            vals_r__ = context__.vals_r("X");
            pos__ = 0;
            size_t X_m_mat_lim__ = N;
            size_t X_n_mat_lim__ = K;
            for (size_t n_mat__ = 0; n_mat__ < X_n_mat_lim__; ++n_mat__) {
                for (size_t m_mat__ = 0; m_mat__ < X_m_mat_lim__; ++m_mat__) {
                    size_t X_limit_0__ = dense_X;
                    for (size_t i_0__ = 0; i_0__ < X_limit_0__; ++i_0__) {
                        X[i_0__](m_mat__,n_mat__) = vals_r__[pos__++];
            }
                }
            }
            current_statement_begin__ = 385;
            context__.validate_dims("data initialization", "nnz_X", "int", context__.to_vec());
            nnz_X = int(0);
            vals_i__ = context__.vals_i("nnz_X");
            pos__ = 0;
            nnz_X = vals_i__[pos__++];
            current_statement_begin__ = 386;
            validate_non_negative_index("w_X", "nnz_X", nnz_X);
            context__.validate_dims("data initialization", "w_X", "vector_d", context__.to_vec(nnz_X));
            validate_non_negative_index("w_X", "nnz_X", nnz_X);
            w_X = vector_d(static_cast<Eigen::VectorXd::Index>(nnz_X));
            vals_r__ = context__.vals_r("w_X");
            pos__ = 0;
            size_t w_X_i_vec_lim__ = nnz_X;
            for (size_t i_vec__ = 0; i_vec__ < w_X_i_vec_lim__; ++i_vec__) {
                w_X[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 387;
            validate_non_negative_index("v_X", "nnz_X", nnz_X);
            context__.validate_dims("data initialization", "v_X", "int", context__.to_vec(nnz_X));
            validate_non_negative_index("v_X", "nnz_X", nnz_X);
            v_X = std::vector<int>(nnz_X,int(0));
            vals_i__ = context__.vals_i("v_X");
            pos__ = 0;
            size_t v_X_limit_0__ = nnz_X;
            for (size_t i_0__ = 0; i_0__ < v_X_limit_0__; ++i_0__) {
                v_X[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 389;
            validate_non_negative_index("u_X", "(dense_X ? 0 : (N + 1) )", (dense_X ? 0 : (N + 1) ));
            context__.validate_dims("data initialization", "u_X", "int", context__.to_vec((dense_X ? 0 : (N + 1) )));
            validate_non_negative_index("u_X", "(dense_X ? 0 : (N + 1) )", (dense_X ? 0 : (N + 1) ));
            u_X = std::vector<int>((dense_X ? 0 : (N + 1) ),int(0));
            vals_i__ = context__.vals_i("u_X");
            pos__ = 0;
            size_t u_X_limit_0__ = (dense_X ? 0 : (N + 1) );
            for (size_t i_0__ = 0; i_0__ < u_X_limit_0__; ++i_0__) {
                u_X[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 392;
            context__.validate_dims("data initialization", "K_smooth", "int", context__.to_vec());
            K_smooth = int(0);
            vals_i__ = context__.vals_i("K_smooth");
            pos__ = 0;
            K_smooth = vals_i__[pos__++];
            current_statement_begin__ = 393;
            validate_non_negative_index("S", "N", N);
            validate_non_negative_index("S", "K_smooth", K_smooth);
            context__.validate_dims("data initialization", "S", "matrix_d", context__.to_vec(N,K_smooth));
            validate_non_negative_index("S", "N", N);
            validate_non_negative_index("S", "K_smooth", K_smooth);
            S = matrix_d(static_cast<Eigen::VectorXd::Index>(N),static_cast<Eigen::VectorXd::Index>(K_smooth));
            vals_r__ = context__.vals_r("S");
            pos__ = 0;
            size_t S_m_mat_lim__ = N;
            size_t S_n_mat_lim__ = K_smooth;
            for (size_t n_mat__ = 0; n_mat__ < S_n_mat_lim__; ++n_mat__) {
                for (size_t m_mat__ = 0; m_mat__ < S_m_mat_lim__; ++m_mat__) {
                    S(m_mat__,n_mat__) = vals_r__[pos__++];
                }
            }
            current_statement_begin__ = 394;
            validate_non_negative_index("smooth_map", "K_smooth", K_smooth);
            context__.validate_dims("data initialization", "smooth_map", "int", context__.to_vec(K_smooth));
            validate_non_negative_index("smooth_map", "K_smooth", K_smooth);
            smooth_map = std::vector<int>(K_smooth,int(0));
            vals_i__ = context__.vals_i("smooth_map");
            pos__ = 0;
            size_t smooth_map_limit_0__ = K_smooth;
            for (size_t i_0__ = 0; i_0__ < smooth_map_limit_0__; ++i_0__) {
                smooth_map[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 395;
            validate_non_negative_index("y", "N", N);
            context__.validate_dims("data initialization", "y", "int", context__.to_vec(N));
            validate_non_negative_index("y", "N", N);
            y = std::vector<int>(N,int(0));
            vals_i__ = context__.vals_i("y");
            pos__ = 0;
            size_t y_limit_0__ = N;
            for (size_t i_0__ = 0; i_0__ < y_limit_0__; ++i_0__) {
                y[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 398;
            context__.validate_dims("data initialization", "prior_PD", "int", context__.to_vec());
            prior_PD = int(0);
            vals_i__ = context__.vals_i("prior_PD");
            pos__ = 0;
            prior_PD = vals_i__[pos__++];
            current_statement_begin__ = 401;
            context__.validate_dims("data initialization", "has_intercept", "int", context__.to_vec());
            has_intercept = int(0);
            vals_i__ = context__.vals_i("has_intercept");
            pos__ = 0;
            has_intercept = vals_i__[pos__++];
            current_statement_begin__ = 404;
            context__.validate_dims("data initialization", "link", "int", context__.to_vec());
            link = int(0);
            vals_i__ = context__.vals_i("link");
            pos__ = 0;
            link = vals_i__[pos__++];
            current_statement_begin__ = 408;
            context__.validate_dims("data initialization", "prior_dist", "int", context__.to_vec());
            prior_dist = int(0);
            vals_i__ = context__.vals_i("prior_dist");
            pos__ = 0;
            prior_dist = vals_i__[pos__++];
            current_statement_begin__ = 409;
            context__.validate_dims("data initialization", "prior_dist_for_intercept", "int", context__.to_vec());
            prior_dist_for_intercept = int(0);
            vals_i__ = context__.vals_i("prior_dist_for_intercept");
            pos__ = 0;
            prior_dist_for_intercept = vals_i__[pos__++];
            current_statement_begin__ = 412;
            context__.validate_dims("data initialization", "prior_dist_for_aux", "int", context__.to_vec());
            prior_dist_for_aux = int(0);
            vals_i__ = context__.vals_i("prior_dist_for_aux");
            pos__ = 0;
            prior_dist_for_aux = vals_i__[pos__++];
            current_statement_begin__ = 415;
            context__.validate_dims("data initialization", "prior_dist_for_smooth", "int", context__.to_vec());
            prior_dist_for_smooth = int(0);
            vals_i__ = context__.vals_i("prior_dist_for_smooth");
            pos__ = 0;
            prior_dist_for_smooth = vals_i__[pos__++];
            current_statement_begin__ = 418;
            context__.validate_dims("data initialization", "has_weights", "int", context__.to_vec());
            has_weights = int(0);
            vals_i__ = context__.vals_i("has_weights");
            pos__ = 0;
            has_weights = vals_i__[pos__++];
            current_statement_begin__ = 419;
            validate_non_negative_index("weights", "(has_weights ? N : 0 )", (has_weights ? N : 0 ));
            context__.validate_dims("data initialization", "weights", "vector_d", context__.to_vec((has_weights ? N : 0 )));
            validate_non_negative_index("weights", "(has_weights ? N : 0 )", (has_weights ? N : 0 ));
            weights = vector_d(static_cast<Eigen::VectorXd::Index>((has_weights ? N : 0 )));
            vals_r__ = context__.vals_r("weights");
            pos__ = 0;
            size_t weights_i_vec_lim__ = (has_weights ? N : 0 );
            for (size_t i_vec__ = 0; i_vec__ < weights_i_vec_lim__; ++i_vec__) {
                weights[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 422;
            context__.validate_dims("data initialization", "has_offset", "int", context__.to_vec());
            has_offset = int(0);
            vals_i__ = context__.vals_i("has_offset");
            pos__ = 0;
            has_offset = vals_i__[pos__++];
            current_statement_begin__ = 423;
            validate_non_negative_index("offset", "(has_offset ? N : 0 )", (has_offset ? N : 0 ));
            context__.validate_dims("data initialization", "offset", "vector_d", context__.to_vec((has_offset ? N : 0 )));
            validate_non_negative_index("offset", "(has_offset ? N : 0 )", (has_offset ? N : 0 ));
            offset = vector_d(static_cast<Eigen::VectorXd::Index>((has_offset ? N : 0 )));
            vals_r__ = context__.vals_r("offset");
            pos__ = 0;
            size_t offset_i_vec_lim__ = (has_offset ? N : 0 );
            for (size_t i_vec__ = 0; i_vec__ < offset_i_vec_lim__; ++i_vec__) {
                offset[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 424;
            context__.validate_dims("data initialization", "family", "int", context__.to_vec());
            family = int(0);
            vals_i__ = context__.vals_i("family");
            pos__ = 0;
            family = vals_i__[pos__++];
            current_statement_begin__ = 427;
            validate_non_negative_index("prior_scale", "K", K);
            context__.validate_dims("data initialization", "prior_scale", "vector_d", context__.to_vec(K));
            validate_non_negative_index("prior_scale", "K", K);
            prior_scale = vector_d(static_cast<Eigen::VectorXd::Index>(K));
            vals_r__ = context__.vals_r("prior_scale");
            pos__ = 0;
            size_t prior_scale_i_vec_lim__ = K;
            for (size_t i_vec__ = 0; i_vec__ < prior_scale_i_vec_lim__; ++i_vec__) {
                prior_scale[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 428;
            context__.validate_dims("data initialization", "prior_scale_for_intercept", "double", context__.to_vec());
            prior_scale_for_intercept = double(0);
            vals_r__ = context__.vals_r("prior_scale_for_intercept");
            pos__ = 0;
            prior_scale_for_intercept = vals_r__[pos__++];
            current_statement_begin__ = 429;
            context__.validate_dims("data initialization", "prior_scale_for_aux", "double", context__.to_vec());
            prior_scale_for_aux = double(0);
            vals_r__ = context__.vals_r("prior_scale_for_aux");
            pos__ = 0;
            prior_scale_for_aux = vals_r__[pos__++];
            current_statement_begin__ = 430;
            validate_non_negative_index("prior_scale_for_smooth", "(logical_gt(K_smooth,0) ? max(smooth_map) : 0 )", (logical_gt(K_smooth,0) ? max(smooth_map) : 0 ));
            context__.validate_dims("data initialization", "prior_scale_for_smooth", "vector_d", context__.to_vec((logical_gt(K_smooth,0) ? max(smooth_map) : 0 )));
            validate_non_negative_index("prior_scale_for_smooth", "(logical_gt(K_smooth,0) ? max(smooth_map) : 0 )", (logical_gt(K_smooth,0) ? max(smooth_map) : 0 ));
            prior_scale_for_smooth = vector_d(static_cast<Eigen::VectorXd::Index>((logical_gt(K_smooth,0) ? max(smooth_map) : 0 )));
            vals_r__ = context__.vals_r("prior_scale_for_smooth");
            pos__ = 0;
            size_t prior_scale_for_smooth_i_vec_lim__ = (logical_gt(K_smooth,0) ? max(smooth_map) : 0 );
            for (size_t i_vec__ = 0; i_vec__ < prior_scale_for_smooth_i_vec_lim__; ++i_vec__) {
                prior_scale_for_smooth[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 431;
            validate_non_negative_index("prior_mean", "K", K);
            context__.validate_dims("data initialization", "prior_mean", "vector_d", context__.to_vec(K));
            validate_non_negative_index("prior_mean", "K", K);
            prior_mean = vector_d(static_cast<Eigen::VectorXd::Index>(K));
            vals_r__ = context__.vals_r("prior_mean");
            pos__ = 0;
            size_t prior_mean_i_vec_lim__ = K;
            for (size_t i_vec__ = 0; i_vec__ < prior_mean_i_vec_lim__; ++i_vec__) {
                prior_mean[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 432;
            context__.validate_dims("data initialization", "prior_mean_for_intercept", "double", context__.to_vec());
            prior_mean_for_intercept = double(0);
            vals_r__ = context__.vals_r("prior_mean_for_intercept");
            pos__ = 0;
            prior_mean_for_intercept = vals_r__[pos__++];
            current_statement_begin__ = 433;
            context__.validate_dims("data initialization", "prior_mean_for_aux", "double", context__.to_vec());
            prior_mean_for_aux = double(0);
            vals_r__ = context__.vals_r("prior_mean_for_aux");
            pos__ = 0;
            prior_mean_for_aux = vals_r__[pos__++];
            current_statement_begin__ = 434;
            validate_non_negative_index("prior_mean_for_smooth", "(logical_gt(K_smooth,0) ? max(smooth_map) : 0 )", (logical_gt(K_smooth,0) ? max(smooth_map) : 0 ));
            context__.validate_dims("data initialization", "prior_mean_for_smooth", "vector_d", context__.to_vec((logical_gt(K_smooth,0) ? max(smooth_map) : 0 )));
            validate_non_negative_index("prior_mean_for_smooth", "(logical_gt(K_smooth,0) ? max(smooth_map) : 0 )", (logical_gt(K_smooth,0) ? max(smooth_map) : 0 ));
            prior_mean_for_smooth = vector_d(static_cast<Eigen::VectorXd::Index>((logical_gt(K_smooth,0) ? max(smooth_map) : 0 )));
            vals_r__ = context__.vals_r("prior_mean_for_smooth");
            pos__ = 0;
            size_t prior_mean_for_smooth_i_vec_lim__ = (logical_gt(K_smooth,0) ? max(smooth_map) : 0 );
            for (size_t i_vec__ = 0; i_vec__ < prior_mean_for_smooth_i_vec_lim__; ++i_vec__) {
                prior_mean_for_smooth[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 435;
            validate_non_negative_index("prior_df", "K", K);
            context__.validate_dims("data initialization", "prior_df", "vector_d", context__.to_vec(K));
            validate_non_negative_index("prior_df", "K", K);
            prior_df = vector_d(static_cast<Eigen::VectorXd::Index>(K));
            vals_r__ = context__.vals_r("prior_df");
            pos__ = 0;
            size_t prior_df_i_vec_lim__ = K;
            for (size_t i_vec__ = 0; i_vec__ < prior_df_i_vec_lim__; ++i_vec__) {
                prior_df[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 436;
            context__.validate_dims("data initialization", "prior_df_for_intercept", "double", context__.to_vec());
            prior_df_for_intercept = double(0);
            vals_r__ = context__.vals_r("prior_df_for_intercept");
            pos__ = 0;
            prior_df_for_intercept = vals_r__[pos__++];
            current_statement_begin__ = 437;
            context__.validate_dims("data initialization", "prior_df_for_aux", "double", context__.to_vec());
            prior_df_for_aux = double(0);
            vals_r__ = context__.vals_r("prior_df_for_aux");
            pos__ = 0;
            prior_df_for_aux = vals_r__[pos__++];
            current_statement_begin__ = 438;
            validate_non_negative_index("prior_df_for_smooth", "(logical_gt(K_smooth,0) ? max(smooth_map) : 0 )", (logical_gt(K_smooth,0) ? max(smooth_map) : 0 ));
            context__.validate_dims("data initialization", "prior_df_for_smooth", "vector_d", context__.to_vec((logical_gt(K_smooth,0) ? max(smooth_map) : 0 )));
            validate_non_negative_index("prior_df_for_smooth", "(logical_gt(K_smooth,0) ? max(smooth_map) : 0 )", (logical_gt(K_smooth,0) ? max(smooth_map) : 0 ));
            prior_df_for_smooth = vector_d(static_cast<Eigen::VectorXd::Index>((logical_gt(K_smooth,0) ? max(smooth_map) : 0 )));
            vals_r__ = context__.vals_r("prior_df_for_smooth");
            pos__ = 0;
            size_t prior_df_for_smooth_i_vec_lim__ = (logical_gt(K_smooth,0) ? max(smooth_map) : 0 );
            for (size_t i_vec__ = 0; i_vec__ < prior_df_for_smooth_i_vec_lim__; ++i_vec__) {
                prior_df_for_smooth[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 439;
            context__.validate_dims("data initialization", "global_prior_df", "double", context__.to_vec());
            global_prior_df = double(0);
            vals_r__ = context__.vals_r("global_prior_df");
            pos__ = 0;
            global_prior_df = vals_r__[pos__++];
            current_statement_begin__ = 440;
            context__.validate_dims("data initialization", "global_prior_scale", "double", context__.to_vec());
            global_prior_scale = double(0);
            vals_r__ = context__.vals_r("global_prior_scale");
            pos__ = 0;
            global_prior_scale = vals_r__[pos__++];
            current_statement_begin__ = 441;
            context__.validate_dims("data initialization", "slab_df", "double", context__.to_vec());
            slab_df = double(0);
            vals_r__ = context__.vals_r("slab_df");
            pos__ = 0;
            slab_df = vals_r__[pos__++];
            current_statement_begin__ = 442;
            context__.validate_dims("data initialization", "slab_scale", "double", context__.to_vec());
            slab_scale = double(0);
            vals_r__ = context__.vals_r("slab_scale");
            pos__ = 0;
            slab_scale = vals_r__[pos__++];
            current_statement_begin__ = 443;
            validate_non_negative_index("num_normals", "(logical_eq(prior_dist,7) ? K : 0 )", (logical_eq(prior_dist,7) ? K : 0 ));
            context__.validate_dims("data initialization", "num_normals", "int", context__.to_vec((logical_eq(prior_dist,7) ? K : 0 )));
            validate_non_negative_index("num_normals", "(logical_eq(prior_dist,7) ? K : 0 )", (logical_eq(prior_dist,7) ? K : 0 ));
            num_normals = std::vector<int>((logical_eq(prior_dist,7) ? K : 0 ),int(0));
            vals_i__ = context__.vals_i("num_normals");
            pos__ = 0;
            size_t num_normals_limit_0__ = (logical_eq(prior_dist,7) ? K : 0 );
            for (size_t i_0__ = 0; i_0__ < num_normals_limit_0__; ++i_0__) {
                num_normals[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 447;
            context__.validate_dims("data initialization", "t", "int", context__.to_vec());
            t = int(0);
            vals_i__ = context__.vals_i("t");
            pos__ = 0;
            t = vals_i__[pos__++];
            current_statement_begin__ = 448;
            validate_non_negative_index("p", "t", t);
            context__.validate_dims("data initialization", "p", "int", context__.to_vec(t));
            validate_non_negative_index("p", "t", t);
            p = std::vector<int>(t,int(0));
            vals_i__ = context__.vals_i("p");
            pos__ = 0;
            size_t p_limit_0__ = t;
            for (size_t i_0__ = 0; i_0__ < p_limit_0__; ++i_0__) {
                p[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 449;
            validate_non_negative_index("l", "t", t);
            context__.validate_dims("data initialization", "l", "int", context__.to_vec(t));
            validate_non_negative_index("l", "t", t);
            l = std::vector<int>(t,int(0));
            vals_i__ = context__.vals_i("l");
            pos__ = 0;
            size_t l_limit_0__ = t;
            for (size_t i_0__ = 0; i_0__ < l_limit_0__; ++i_0__) {
                l[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 450;
            context__.validate_dims("data initialization", "q", "int", context__.to_vec());
            q = int(0);
            vals_i__ = context__.vals_i("q");
            pos__ = 0;
            q = vals_i__[pos__++];
            current_statement_begin__ = 451;
            context__.validate_dims("data initialization", "len_theta_L", "int", context__.to_vec());
            len_theta_L = int(0);
            vals_i__ = context__.vals_i("len_theta_L");
            pos__ = 0;
            len_theta_L = vals_i__[pos__++];
            current_statement_begin__ = 454;
            validate_non_negative_index("shape", "t", t);
            context__.validate_dims("data initialization", "shape", "vector_d", context__.to_vec(t));
            validate_non_negative_index("shape", "t", t);
            shape = vector_d(static_cast<Eigen::VectorXd::Index>(t));
            vals_r__ = context__.vals_r("shape");
            pos__ = 0;
            size_t shape_i_vec_lim__ = t;
            for (size_t i_vec__ = 0; i_vec__ < shape_i_vec_lim__; ++i_vec__) {
                shape[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 455;
            validate_non_negative_index("scale", "t", t);
            context__.validate_dims("data initialization", "scale", "vector_d", context__.to_vec(t));
            validate_non_negative_index("scale", "t", t);
            scale = vector_d(static_cast<Eigen::VectorXd::Index>(t));
            vals_r__ = context__.vals_r("scale");
            pos__ = 0;
            size_t scale_i_vec_lim__ = t;
            for (size_t i_vec__ = 0; i_vec__ < scale_i_vec_lim__; ++i_vec__) {
                scale[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 456;
            context__.validate_dims("data initialization", "len_concentration", "int", context__.to_vec());
            len_concentration = int(0);
            vals_i__ = context__.vals_i("len_concentration");
            pos__ = 0;
            len_concentration = vals_i__[pos__++];
            current_statement_begin__ = 457;
            validate_non_negative_index("concentration", "len_concentration", len_concentration);
            context__.validate_dims("data initialization", "concentration", "double", context__.to_vec(len_concentration));
            validate_non_negative_index("concentration", "len_concentration", len_concentration);
            concentration = std::vector<double>(len_concentration,double(0));
            vals_r__ = context__.vals_r("concentration");
            pos__ = 0;
            size_t concentration_limit_0__ = len_concentration;
            for (size_t i_0__ = 0; i_0__ < concentration_limit_0__; ++i_0__) {
                concentration[i_0__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 458;
            context__.validate_dims("data initialization", "len_regularization", "int", context__.to_vec());
            len_regularization = int(0);
            vals_i__ = context__.vals_i("len_regularization");
            pos__ = 0;
            len_regularization = vals_i__[pos__++];
            current_statement_begin__ = 459;
            validate_non_negative_index("regularization", "len_regularization", len_regularization);
            context__.validate_dims("data initialization", "regularization", "double", context__.to_vec(len_regularization));
            validate_non_negative_index("regularization", "len_regularization", len_regularization);
            regularization = std::vector<double>(len_regularization,double(0));
            vals_r__ = context__.vals_r("regularization");
            pos__ = 0;
            size_t regularization_limit_0__ = len_regularization;
            for (size_t i_0__ = 0; i_0__ < regularization_limit_0__; ++i_0__) {
                regularization[i_0__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 461;
            context__.validate_dims("data initialization", "num_non_zero", "int", context__.to_vec());
            num_non_zero = int(0);
            vals_i__ = context__.vals_i("num_non_zero");
            pos__ = 0;
            num_non_zero = vals_i__[pos__++];
            current_statement_begin__ = 462;
            validate_non_negative_index("w", "num_non_zero", num_non_zero);
            context__.validate_dims("data initialization", "w", "vector_d", context__.to_vec(num_non_zero));
            validate_non_negative_index("w", "num_non_zero", num_non_zero);
            w = vector_d(static_cast<Eigen::VectorXd::Index>(num_non_zero));
            vals_r__ = context__.vals_r("w");
            pos__ = 0;
            size_t w_i_vec_lim__ = num_non_zero;
            for (size_t i_vec__ = 0; i_vec__ < w_i_vec_lim__; ++i_vec__) {
                w[i_vec__] = vals_r__[pos__++];
            }
            current_statement_begin__ = 463;
            validate_non_negative_index("v", "num_non_zero", num_non_zero);
            context__.validate_dims("data initialization", "v", "int", context__.to_vec(num_non_zero));
            validate_non_negative_index("v", "num_non_zero", num_non_zero);
            v = std::vector<int>(num_non_zero,int(0));
            vals_i__ = context__.vals_i("v");
            pos__ = 0;
            size_t v_limit_0__ = num_non_zero;
            for (size_t i_0__ = 0; i_0__ < v_limit_0__; ++i_0__) {
                v[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 464;
            validate_non_negative_index("u", "(logical_gt(t,0) ? (N + 1) : 0 )", (logical_gt(t,0) ? (N + 1) : 0 ));
            context__.validate_dims("data initialization", "u", "int", context__.to_vec((logical_gt(t,0) ? (N + 1) : 0 )));
            validate_non_negative_index("u", "(logical_gt(t,0) ? (N + 1) : 0 )", (logical_gt(t,0) ? (N + 1) : 0 ));
            u = std::vector<int>((logical_gt(t,0) ? (N + 1) : 0 ),int(0));
            vals_i__ = context__.vals_i("u");
            pos__ = 0;
            size_t u_limit_0__ = (logical_gt(t,0) ? (N + 1) : 0 );
            for (size_t i_0__ = 0; i_0__ < u_limit_0__; ++i_0__) {
                u[i_0__] = vals_i__[pos__++];
            }
            current_statement_begin__ = 465;
            context__.validate_dims("data initialization", "special_case", "int", context__.to_vec());
            special_case = int(0);
            vals_i__ = context__.vals_i("special_case");
            pos__ = 0;
            special_case = vals_i__[pos__++];

            // validate, data variables
            current_statement_begin__ = 376;
            check_greater_or_equal(function__,"N",N,0);
            current_statement_begin__ = 377;
            check_greater_or_equal(function__,"K",K,0);
            current_statement_begin__ = 380;
            current_statement_begin__ = 381;
            check_greater_or_equal(function__,"dense_X",dense_X,0);
            check_less_or_equal(function__,"dense_X",dense_X,1);
            current_statement_begin__ = 382;
            current_statement_begin__ = 385;
            check_greater_or_equal(function__,"nnz_X",nnz_X,0);
            current_statement_begin__ = 386;
            current_statement_begin__ = 387;
            for (int k0__ = 0; k0__ < nnz_X; ++k0__) {
                check_greater_or_equal(function__,"v_X[k0__]",v_X[k0__],0);
                check_less_or_equal(function__,"v_X[k0__]",v_X[k0__],(K - 1));
            }
            current_statement_begin__ = 389;
            for (int k0__ = 0; k0__ < (dense_X ? 0 : (N + 1) ); ++k0__) {
                check_greater_or_equal(function__,"u_X[k0__]",u_X[k0__],0);
                check_less_or_equal(function__,"u_X[k0__]",u_X[k0__],(rows(w_X) + 1));
            }
            current_statement_begin__ = 392;
            check_greater_or_equal(function__,"K_smooth",K_smooth,0);
            current_statement_begin__ = 393;
            current_statement_begin__ = 394;
            for (int k0__ = 0; k0__ < K_smooth; ++k0__) {
                check_greater_or_equal(function__,"smooth_map[k0__]",smooth_map[k0__],1);
            }
            current_statement_begin__ = 395;
            for (int k0__ = 0; k0__ < N; ++k0__) {
                check_greater_or_equal(function__,"y[k0__]",y[k0__],0);
            }
            current_statement_begin__ = 398;
            check_greater_or_equal(function__,"prior_PD",prior_PD,0);
            check_less_or_equal(function__,"prior_PD",prior_PD,1);
            current_statement_begin__ = 401;
            check_greater_or_equal(function__,"has_intercept",has_intercept,0);
            check_less_or_equal(function__,"has_intercept",has_intercept,1);
            current_statement_begin__ = 404;
            check_greater_or_equal(function__,"link",link,1);
            current_statement_begin__ = 408;
            check_greater_or_equal(function__,"prior_dist",prior_dist,0);
            check_less_or_equal(function__,"prior_dist",prior_dist,7);
            current_statement_begin__ = 409;
            check_greater_or_equal(function__,"prior_dist_for_intercept",prior_dist_for_intercept,0);
            check_less_or_equal(function__,"prior_dist_for_intercept",prior_dist_for_intercept,2);
            current_statement_begin__ = 412;
            check_greater_or_equal(function__,"prior_dist_for_aux",prior_dist_for_aux,0);
            check_less_or_equal(function__,"prior_dist_for_aux",prior_dist_for_aux,3);
            current_statement_begin__ = 415;
            check_greater_or_equal(function__,"prior_dist_for_smooth",prior_dist_for_smooth,0);
            check_less_or_equal(function__,"prior_dist_for_smooth",prior_dist_for_smooth,3);
            current_statement_begin__ = 418;
            check_greater_or_equal(function__,"has_weights",has_weights,0);
            check_less_or_equal(function__,"has_weights",has_weights,1);
            current_statement_begin__ = 419;
            current_statement_begin__ = 422;
            check_greater_or_equal(function__,"has_offset",has_offset,0);
            check_less_or_equal(function__,"has_offset",has_offset,1);
            current_statement_begin__ = 423;
            current_statement_begin__ = 424;
            check_greater_or_equal(function__,"family",family,6);
            check_less_or_equal(function__,"family",family,7);
            current_statement_begin__ = 427;
            check_greater_or_equal(function__,"prior_scale",prior_scale,0);
            current_statement_begin__ = 428;
            check_greater_or_equal(function__,"prior_scale_for_intercept",prior_scale_for_intercept,0);
            current_statement_begin__ = 429;
            check_greater_or_equal(function__,"prior_scale_for_aux",prior_scale_for_aux,0);
            current_statement_begin__ = 430;
            check_greater_or_equal(function__,"prior_scale_for_smooth",prior_scale_for_smooth,0);
            current_statement_begin__ = 431;
            current_statement_begin__ = 432;
            current_statement_begin__ = 433;
            check_greater_or_equal(function__,"prior_mean_for_aux",prior_mean_for_aux,0);
            current_statement_begin__ = 434;
            check_greater_or_equal(function__,"prior_mean_for_smooth",prior_mean_for_smooth,0);
            current_statement_begin__ = 435;
            check_greater_or_equal(function__,"prior_df",prior_df,0);
            current_statement_begin__ = 436;
            check_greater_or_equal(function__,"prior_df_for_intercept",prior_df_for_intercept,0);
            current_statement_begin__ = 437;
            check_greater_or_equal(function__,"prior_df_for_aux",prior_df_for_aux,0);
            current_statement_begin__ = 438;
            check_greater_or_equal(function__,"prior_df_for_smooth",prior_df_for_smooth,0);
            current_statement_begin__ = 439;
            check_greater_or_equal(function__,"global_prior_df",global_prior_df,0);
            current_statement_begin__ = 440;
            check_greater_or_equal(function__,"global_prior_scale",global_prior_scale,0);
            current_statement_begin__ = 441;
            check_greater_or_equal(function__,"slab_df",slab_df,0);
            current_statement_begin__ = 442;
            check_greater_or_equal(function__,"slab_scale",slab_scale,0);
            current_statement_begin__ = 443;
            for (int k0__ = 0; k0__ < (logical_eq(prior_dist,7) ? K : 0 ); ++k0__) {
                check_greater_or_equal(function__,"num_normals[k0__]",num_normals[k0__],2);
            }
            current_statement_begin__ = 447;
            check_greater_or_equal(function__,"t",t,0);
            current_statement_begin__ = 448;
            for (int k0__ = 0; k0__ < t; ++k0__) {
                check_greater_or_equal(function__,"p[k0__]",p[k0__],1);
            }
            current_statement_begin__ = 449;
            for (int k0__ = 0; k0__ < t; ++k0__) {
                check_greater_or_equal(function__,"l[k0__]",l[k0__],1);
            }
            current_statement_begin__ = 450;
            check_greater_or_equal(function__,"q",q,0);
            current_statement_begin__ = 451;
            check_greater_or_equal(function__,"len_theta_L",len_theta_L,0);
            current_statement_begin__ = 454;
            check_greater_or_equal(function__,"shape",shape,0);
            current_statement_begin__ = 455;
            check_greater_or_equal(function__,"scale",scale,0);
            current_statement_begin__ = 456;
            check_greater_or_equal(function__,"len_concentration",len_concentration,0);
            current_statement_begin__ = 457;
            for (int k0__ = 0; k0__ < len_concentration; ++k0__) {
                check_greater_or_equal(function__,"concentration[k0__]",concentration[k0__],0);
            }
            current_statement_begin__ = 458;
            check_greater_or_equal(function__,"len_regularization",len_regularization,0);
            current_statement_begin__ = 459;
            for (int k0__ = 0; k0__ < len_regularization; ++k0__) {
                check_greater_or_equal(function__,"regularization[k0__]",regularization[k0__],0);
            }
            current_statement_begin__ = 461;
            check_greater_or_equal(function__,"num_non_zero",num_non_zero,0);
            current_statement_begin__ = 462;
            current_statement_begin__ = 463;
            for (int k0__ = 0; k0__ < num_non_zero; ++k0__) {
                check_greater_or_equal(function__,"v[k0__]",v[k0__],0);
                check_less_or_equal(function__,"v[k0__]",v[k0__],(q - 1));
            }
            current_statement_begin__ = 464;
            for (int k0__ = 0; k0__ < (logical_gt(t,0) ? (N + 1) : 0 ); ++k0__) {
                check_greater_or_equal(function__,"u[k0__]",u[k0__],0);
                check_less_or_equal(function__,"u[k0__]",u[k0__],(rows(w) + 1));
            }
            current_statement_begin__ = 465;
            check_greater_or_equal(function__,"special_case",special_case,0);
            check_less_or_equal(function__,"special_case",special_case,1);
            // initialize data variables
            current_statement_begin__ = 468;
            poisson_max = double(0);
            stan::math::fill(poisson_max,DUMMY_VAR__);
            stan::math::assign(poisson_max,pow(2.0,30.0));
            current_statement_begin__ = 469;
            validate_non_negative_index("V", "(special_case ? t : 0 )", (special_case ? t : 0 ));
            validate_non_negative_index("V", "N", N);
            V = std::vector<std::vector<int> >((special_case ? t : 0 ),std::vector<int>(N,int(0)));
            stan::math::fill(V, std::numeric_limits<int>::min());
            stan::math::assign(V,make_V(N,(special_case ? t : 0 ),v, pstream__));
            current_statement_begin__ = 471;
            len_z_T = int(0);
            stan::math::fill(len_z_T, std::numeric_limits<int>::min());
            stan::math::assign(len_z_T,0);
            current_statement_begin__ = 472;
            len_var_group = int(0);
            stan::math::fill(len_var_group, std::numeric_limits<int>::min());
            stan::math::assign(len_var_group,(sum(p) * logical_gt(t,0)));
            current_statement_begin__ = 473;
            len_rho = int(0);
            stan::math::fill(len_rho, std::numeric_limits<int>::min());
            stan::math::assign(len_rho,(sum(p) - t));
            current_statement_begin__ = 474;
            is_continuous = int(0);
            stan::math::fill(is_continuous, std::numeric_limits<int>::min());
            stan::math::assign(is_continuous,0);
            current_statement_begin__ = 475;
            pos = int(0);
            stan::math::fill(pos, std::numeric_limits<int>::min());
            stan::math::assign(pos,1);
            current_statement_begin__ = 476;
            validate_non_negative_index("delta", "len_concentration", len_concentration);
            delta = std::vector<double>(len_concentration,double(0));
            stan::math::fill(delta,DUMMY_VAR__);
            current_statement_begin__ = 477;
            hs = int(0);
            stan::math::fill(hs, std::numeric_limits<int>::min());

            current_statement_begin__ = 478;
            if (as_bool(logical_lte(prior_dist,2))) {
                current_statement_begin__ = 478;
                stan::math::assign(hs, 0);
            } else if (as_bool(logical_eq(prior_dist,3))) {
                current_statement_begin__ = 479;
                stan::math::assign(hs, 2);
            } else if (as_bool(logical_eq(prior_dist,4))) {
                current_statement_begin__ = 480;
                stan::math::assign(hs, 4);
            } else {
                current_statement_begin__ = 481;
                stan::math::assign(hs, 0);
            }
            current_statement_begin__ = 483;
            stan::math::assign(pos, 1);
            current_statement_begin__ = 484;
            for (int i = 1; i <= t; ++i) {

                current_statement_begin__ = 485;
                if (as_bool(logical_gt(get_base1(p,i,"p",1),1))) {

                    current_statement_begin__ = 486;
                    for (int j = 1; j <= get_base1(p,i,"p",1); ++j) {

                        current_statement_begin__ = 487;
                        stan::math::assign(get_base1_lhs(delta,pos,"delta",1), get_base1(concentration,j,"concentration",1));
                        current_statement_begin__ = 488;
                        stan::math::assign(pos, (pos + 1));
                    }
                }
                current_statement_begin__ = 491;
                for (int j = 3; j <= get_base1(p,i,"p",1); ++j) {
                    current_statement_begin__ = 491;
                    stan::math::assign(len_z_T, ((len_z_T + get_base1(p,i,"p",1)) - 1));
                }
            }

            // validate transformed data
            current_statement_begin__ = 468;
            current_statement_begin__ = 469;
            for (int k0__ = 0; k0__ < (special_case ? t : 0 ); ++k0__) {
                for (int k1__ = 0; k1__ < N; ++k1__) {
                    check_greater_or_equal(function__,"V[k0__][k1__]",V[k0__][k1__],1);
                }
            }
            current_statement_begin__ = 471;
            check_greater_or_equal(function__,"len_z_T",len_z_T,0);
            current_statement_begin__ = 472;
            check_greater_or_equal(function__,"len_var_group",len_var_group,0);
            current_statement_begin__ = 473;
            check_greater_or_equal(function__,"len_rho",len_rho,0);
            current_statement_begin__ = 474;
            check_greater_or_equal(function__,"is_continuous",is_continuous,0);
            check_less_or_equal(function__,"is_continuous",is_continuous,1);
            current_statement_begin__ = 475;
            check_greater_or_equal(function__,"pos",pos,1);
            current_statement_begin__ = 476;
            for (int k0__ = 0; k0__ < len_concentration; ++k0__) {
                check_greater_or_equal(function__,"delta[k0__]",delta[k0__],0);
            }
            current_statement_begin__ = 477;
            check_greater_or_equal(function__,"hs",hs,0);

            // validate, set parameter ranges
            num_params_r__ = 0U;
            param_ranges_i__.clear();
            current_statement_begin__ = 495;
            validate_non_negative_index("gamma", "has_intercept", has_intercept);
            num_params_r__ += has_intercept;
            current_statement_begin__ = 497;
            validate_non_negative_index("z_beta", "(logical_eq(prior_dist,7) ? sum(num_normals) : K )", (logical_eq(prior_dist,7) ? sum(num_normals) : K ));
            num_params_r__ += (logical_eq(prior_dist,7) ? sum(num_normals) : K );
            current_statement_begin__ = 498;
            validate_non_negative_index("z_beta_smooth", "K_smooth", K_smooth);
            num_params_r__ += K_smooth;
            current_statement_begin__ = 499;
            validate_non_negative_index("smooth_sd_raw", "(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,\"smooth_map\",1) : 0 )", (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
            num_params_r__ += (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 );
            current_statement_begin__ = 500;
            validate_non_negative_index("global", "hs", hs);
            num_params_r__ += hs;
            current_statement_begin__ = 501;
            validate_non_negative_index("local", "K", K);
            validate_non_negative_index("local", "hs", hs);
            num_params_r__ += K * hs;
            current_statement_begin__ = 502;
            validate_non_negative_index("caux", "logical_gt(hs,0)", logical_gt(hs,0));
            num_params_r__ += logical_gt(hs,0);
            current_statement_begin__ = 503;
            validate_non_negative_index("mix", "K", K);
            validate_non_negative_index("mix", "(primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6)))", (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))));
            num_params_r__ += K * (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6)));
            current_statement_begin__ = 504;
            validate_non_negative_index("one_over_lambda", "logical_eq(prior_dist,6)", logical_eq(prior_dist,6));
            num_params_r__ += logical_eq(prior_dist,6);
            current_statement_begin__ = 505;
            validate_non_negative_index("z_b", "q", q);
            num_params_r__ += q;
            current_statement_begin__ = 506;
            validate_non_negative_index("z_T", "len_z_T", len_z_T);
            num_params_r__ += len_z_T;
            current_statement_begin__ = 507;
            validate_non_negative_index("rho", "len_rho", len_rho);
            num_params_r__ += len_rho;
            current_statement_begin__ = 508;
            validate_non_negative_index("zeta", "len_concentration", len_concentration);
            num_params_r__ += len_concentration;
            current_statement_begin__ = 509;
            validate_non_negative_index("tau", "t", t);
            num_params_r__ += t;
            current_statement_begin__ = 510;
            validate_non_negative_index("aux_unscaled", "logical_gt(family,6)", logical_gt(family,6));
            num_params_r__ += logical_gt(family,6);
            current_statement_begin__ = 511;
            validate_non_negative_index("noise", "N", N);
            validate_non_negative_index("noise", "logical_eq(family,8)", logical_eq(family,8));
            num_params_r__ += N * logical_eq(family,8);
        } 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_count() { }


    void transform_inits(const stan::io::var_context& context__,
                         std::vector<int>& params_i__,
                         std::vector<double>& params_r__,
                         std::ostream* pstream__) const {
        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__;

        if (!(context__.contains_r("gamma")))
            throw std::runtime_error("variable gamma missing");
        vals_r__ = context__.vals_r("gamma");
        pos__ = 0U;
        validate_non_negative_index("gamma", "has_intercept", has_intercept);
        context__.validate_dims("initialization", "gamma", "double", context__.to_vec(has_intercept));
        std::vector<double> gamma(has_intercept,double(0));
        for (int i0__ = 0U; i0__ < has_intercept; ++i0__)
            gamma[i0__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < has_intercept; ++i0__)
            try {
            writer__.scalar_lb_unconstrain((logical_eq(link,1) ? stan::math::negative_infinity() : 0.0 ),gamma[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable gamma: ") + e.what());
        }

        if (!(context__.contains_r("z_beta")))
            throw std::runtime_error("variable z_beta missing");
        vals_r__ = context__.vals_r("z_beta");
        pos__ = 0U;
        validate_non_negative_index("z_beta", "(logical_eq(prior_dist,7) ? sum(num_normals) : K )", (logical_eq(prior_dist,7) ? sum(num_normals) : K ));
        context__.validate_dims("initialization", "z_beta", "vector_d", context__.to_vec((logical_eq(prior_dist,7) ? sum(num_normals) : K )));
        vector_d z_beta(static_cast<Eigen::VectorXd::Index>((logical_eq(prior_dist,7) ? sum(num_normals) : K )));
        for (int j1__ = 0U; j1__ < (logical_eq(prior_dist,7) ? sum(num_normals) : K ); ++j1__)
            z_beta(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_unconstrain(z_beta);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable z_beta: ") + e.what());
        }

        if (!(context__.contains_r("z_beta_smooth")))
            throw std::runtime_error("variable z_beta_smooth missing");
        vals_r__ = context__.vals_r("z_beta_smooth");
        pos__ = 0U;
        validate_non_negative_index("z_beta_smooth", "K_smooth", K_smooth);
        context__.validate_dims("initialization", "z_beta_smooth", "vector_d", context__.to_vec(K_smooth));
        vector_d z_beta_smooth(static_cast<Eigen::VectorXd::Index>(K_smooth));
        for (int j1__ = 0U; j1__ < K_smooth; ++j1__)
            z_beta_smooth(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_unconstrain(z_beta_smooth);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable z_beta_smooth: ") + e.what());
        }

        if (!(context__.contains_r("smooth_sd_raw")))
            throw std::runtime_error("variable smooth_sd_raw missing");
        vals_r__ = context__.vals_r("smooth_sd_raw");
        pos__ = 0U;
        validate_non_negative_index("smooth_sd_raw", "(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,\"smooth_map\",1) : 0 )", (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
        context__.validate_dims("initialization", "smooth_sd_raw", "vector_d", context__.to_vec((logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 )));
        vector_d smooth_sd_raw(static_cast<Eigen::VectorXd::Index>((logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 )));
        for (int j1__ = 0U; j1__ < (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++j1__)
            smooth_sd_raw(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_lb_unconstrain(0,smooth_sd_raw);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable smooth_sd_raw: ") + e.what());
        }

        if (!(context__.contains_r("global")))
            throw std::runtime_error("variable global missing");
        vals_r__ = context__.vals_r("global");
        pos__ = 0U;
        validate_non_negative_index("global", "hs", hs);
        context__.validate_dims("initialization", "global", "double", context__.to_vec(hs));
        std::vector<double> global(hs,double(0));
        for (int i0__ = 0U; i0__ < hs; ++i0__)
            global[i0__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < hs; ++i0__)
            try {
            writer__.scalar_lb_unconstrain(0,global[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable global: ") + e.what());
        }

        if (!(context__.contains_r("local")))
            throw std::runtime_error("variable local missing");
        vals_r__ = context__.vals_r("local");
        pos__ = 0U;
        validate_non_negative_index("local", "hs", hs);
        validate_non_negative_index("local", "K", K);
        context__.validate_dims("initialization", "local", "vector_d", context__.to_vec(hs,K));
        std::vector<vector_d> local(hs,vector_d(static_cast<Eigen::VectorXd::Index>(K)));
        for (int j1__ = 0U; j1__ < K; ++j1__)
            for (int i0__ = 0U; i0__ < hs; ++i0__)
                local[i0__](j1__) = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < hs; ++i0__)
            try {
            writer__.vector_lb_unconstrain(0,local[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable local: ") + e.what());
        }

        if (!(context__.contains_r("caux")))
            throw std::runtime_error("variable caux missing");
        vals_r__ = context__.vals_r("caux");
        pos__ = 0U;
        validate_non_negative_index("caux", "logical_gt(hs,0)", logical_gt(hs,0));
        context__.validate_dims("initialization", "caux", "double", context__.to_vec(logical_gt(hs,0)));
        std::vector<double> caux(logical_gt(hs,0),double(0));
        for (int i0__ = 0U; i0__ < logical_gt(hs,0); ++i0__)
            caux[i0__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < logical_gt(hs,0); ++i0__)
            try {
            writer__.scalar_lb_unconstrain(0,caux[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable caux: ") + e.what());
        }

        if (!(context__.contains_r("mix")))
            throw std::runtime_error("variable mix missing");
        vals_r__ = context__.vals_r("mix");
        pos__ = 0U;
        validate_non_negative_index("mix", "(primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6)))", (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))));
        validate_non_negative_index("mix", "K", K);
        context__.validate_dims("initialization", "mix", "vector_d", context__.to_vec((primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))),K));
        std::vector<vector_d> mix((primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))),vector_d(static_cast<Eigen::VectorXd::Index>(K)));
        for (int j1__ = 0U; j1__ < K; ++j1__)
            for (int i0__ = 0U; i0__ < (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))); ++i0__)
                mix[i0__](j1__) = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))); ++i0__)
            try {
            writer__.vector_lb_unconstrain(0,mix[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable mix: ") + e.what());
        }

        if (!(context__.contains_r("one_over_lambda")))
            throw std::runtime_error("variable one_over_lambda missing");
        vals_r__ = context__.vals_r("one_over_lambda");
        pos__ = 0U;
        validate_non_negative_index("one_over_lambda", "logical_eq(prior_dist,6)", logical_eq(prior_dist,6));
        context__.validate_dims("initialization", "one_over_lambda", "double", context__.to_vec(logical_eq(prior_dist,6)));
        std::vector<double> one_over_lambda(logical_eq(prior_dist,6),double(0));
        for (int i0__ = 0U; i0__ < logical_eq(prior_dist,6); ++i0__)
            one_over_lambda[i0__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < logical_eq(prior_dist,6); ++i0__)
            try {
            writer__.scalar_lb_unconstrain(0,one_over_lambda[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable one_over_lambda: ") + e.what());
        }

        if (!(context__.contains_r("z_b")))
            throw std::runtime_error("variable z_b missing");
        vals_r__ = context__.vals_r("z_b");
        pos__ = 0U;
        validate_non_negative_index("z_b", "q", q);
        context__.validate_dims("initialization", "z_b", "vector_d", context__.to_vec(q));
        vector_d z_b(static_cast<Eigen::VectorXd::Index>(q));
        for (int j1__ = 0U; j1__ < q; ++j1__)
            z_b(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_unconstrain(z_b);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable z_b: ") + e.what());
        }

        if (!(context__.contains_r("z_T")))
            throw std::runtime_error("variable z_T missing");
        vals_r__ = context__.vals_r("z_T");
        pos__ = 0U;
        validate_non_negative_index("z_T", "len_z_T", len_z_T);
        context__.validate_dims("initialization", "z_T", "vector_d", context__.to_vec(len_z_T));
        vector_d z_T(static_cast<Eigen::VectorXd::Index>(len_z_T));
        for (int j1__ = 0U; j1__ < len_z_T; ++j1__)
            z_T(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_unconstrain(z_T);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable z_T: ") + e.what());
        }

        if (!(context__.contains_r("rho")))
            throw std::runtime_error("variable rho missing");
        vals_r__ = context__.vals_r("rho");
        pos__ = 0U;
        validate_non_negative_index("rho", "len_rho", len_rho);
        context__.validate_dims("initialization", "rho", "vector_d", context__.to_vec(len_rho));
        vector_d rho(static_cast<Eigen::VectorXd::Index>(len_rho));
        for (int j1__ = 0U; j1__ < len_rho; ++j1__)
            rho(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_lub_unconstrain(0,1,rho);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable rho: ") + e.what());
        }

        if (!(context__.contains_r("zeta")))
            throw std::runtime_error("variable zeta missing");
        vals_r__ = context__.vals_r("zeta");
        pos__ = 0U;
        validate_non_negative_index("zeta", "len_concentration", len_concentration);
        context__.validate_dims("initialization", "zeta", "vector_d", context__.to_vec(len_concentration));
        vector_d zeta(static_cast<Eigen::VectorXd::Index>(len_concentration));
        for (int j1__ = 0U; j1__ < len_concentration; ++j1__)
            zeta(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_lb_unconstrain(0,zeta);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable zeta: ") + e.what());
        }

        if (!(context__.contains_r("tau")))
            throw std::runtime_error("variable tau missing");
        vals_r__ = context__.vals_r("tau");
        pos__ = 0U;
        validate_non_negative_index("tau", "t", t);
        context__.validate_dims("initialization", "tau", "vector_d", context__.to_vec(t));
        vector_d tau(static_cast<Eigen::VectorXd::Index>(t));
        for (int j1__ = 0U; j1__ < t; ++j1__)
            tau(j1__) = vals_r__[pos__++];
        try {
            writer__.vector_lb_unconstrain(0,tau);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable tau: ") + e.what());
        }

        if (!(context__.contains_r("aux_unscaled")))
            throw std::runtime_error("variable aux_unscaled missing");
        vals_r__ = context__.vals_r("aux_unscaled");
        pos__ = 0U;
        validate_non_negative_index("aux_unscaled", "logical_gt(family,6)", logical_gt(family,6));
        context__.validate_dims("initialization", "aux_unscaled", "double", context__.to_vec(logical_gt(family,6)));
        std::vector<double> aux_unscaled(logical_gt(family,6),double(0));
        for (int i0__ = 0U; i0__ < logical_gt(family,6); ++i0__)
            aux_unscaled[i0__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < logical_gt(family,6); ++i0__)
            try {
            writer__.scalar_lb_unconstrain(0,aux_unscaled[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable aux_unscaled: ") + e.what());
        }

        if (!(context__.contains_r("noise")))
            throw std::runtime_error("variable noise missing");
        vals_r__ = context__.vals_r("noise");
        pos__ = 0U;
        validate_non_negative_index("noise", "logical_eq(family,8)", logical_eq(family,8));
        validate_non_negative_index("noise", "N", N);
        context__.validate_dims("initialization", "noise", "vector_d", context__.to_vec(logical_eq(family,8),N));
        std::vector<vector_d> noise(logical_eq(family,8),vector_d(static_cast<Eigen::VectorXd::Index>(N)));
        for (int j1__ = 0U; j1__ < N; ++j1__)
            for (int i0__ = 0U; i0__ < logical_eq(family,8); ++i0__)
                noise[i0__](j1__) = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < logical_eq(family,8); ++i0__)
            try {
            writer__.vector_lb_unconstrain(0,noise[i0__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable noise: ") + e.what());
        }

        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(vector<T__>& params_r__,
                 vector<int>& params_i__,
                 std::ostream* pstream__ = 0) const {

        T__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

        T__ lp__(0.0);
        stan::math::accumulator<T__> lp_accum__;

        try {
            // model parameters
            stan::io::reader<T__> in__(params_r__,params_i__);

            vector<T__> gamma;
            size_t dim_gamma_0__ = has_intercept;
            gamma.reserve(dim_gamma_0__);
            for (size_t k_0__ = 0; k_0__ < dim_gamma_0__; ++k_0__) {
                if (jacobian__)
                    gamma.push_back(in__.scalar_lb_constrain((logical_eq(link,1) ? stan::math::negative_infinity() : 0.0 ),lp__));
                else
                    gamma.push_back(in__.scalar_lb_constrain((logical_eq(link,1) ? stan::math::negative_infinity() : 0.0 )));
            }

            Eigen::Matrix<T__,Eigen::Dynamic,1>  z_beta;
            (void) z_beta;  // dummy to suppress unused var warning
            if (jacobian__)
                z_beta = in__.vector_constrain((logical_eq(prior_dist,7) ? sum(num_normals) : K ),lp__);
            else
                z_beta = in__.vector_constrain((logical_eq(prior_dist,7) ? sum(num_normals) : K ));

            Eigen::Matrix<T__,Eigen::Dynamic,1>  z_beta_smooth;
            (void) z_beta_smooth;  // dummy to suppress unused var warning
            if (jacobian__)
                z_beta_smooth = in__.vector_constrain(K_smooth,lp__);
            else
                z_beta_smooth = in__.vector_constrain(K_smooth);

            Eigen::Matrix<T__,Eigen::Dynamic,1>  smooth_sd_raw;
            (void) smooth_sd_raw;  // dummy to suppress unused var warning
            if (jacobian__)
                smooth_sd_raw = in__.vector_lb_constrain(0,(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ),lp__);
            else
                smooth_sd_raw = in__.vector_lb_constrain(0,(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));

            vector<T__> global;
            size_t dim_global_0__ = hs;
            global.reserve(dim_global_0__);
            for (size_t k_0__ = 0; k_0__ < dim_global_0__; ++k_0__) {
                if (jacobian__)
                    global.push_back(in__.scalar_lb_constrain(0,lp__));
                else
                    global.push_back(in__.scalar_lb_constrain(0));
            }

            vector<Eigen::Matrix<T__,Eigen::Dynamic,1> > local;
            size_t dim_local_0__ = hs;
            local.reserve(dim_local_0__);
            for (size_t k_0__ = 0; k_0__ < dim_local_0__; ++k_0__) {
                if (jacobian__)
                    local.push_back(in__.vector_lb_constrain(0,K,lp__));
                else
                    local.push_back(in__.vector_lb_constrain(0,K));
            }

            vector<T__> caux;
            size_t dim_caux_0__ = logical_gt(hs,0);
            caux.reserve(dim_caux_0__);
            for (size_t k_0__ = 0; k_0__ < dim_caux_0__; ++k_0__) {
                if (jacobian__)
                    caux.push_back(in__.scalar_lb_constrain(0,lp__));
                else
                    caux.push_back(in__.scalar_lb_constrain(0));
            }

            vector<Eigen::Matrix<T__,Eigen::Dynamic,1> > mix;
            size_t dim_mix_0__ = (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6)));
            mix.reserve(dim_mix_0__);
            for (size_t k_0__ = 0; k_0__ < dim_mix_0__; ++k_0__) {
                if (jacobian__)
                    mix.push_back(in__.vector_lb_constrain(0,K,lp__));
                else
                    mix.push_back(in__.vector_lb_constrain(0,K));
            }

            vector<T__> one_over_lambda;
            size_t dim_one_over_lambda_0__ = logical_eq(prior_dist,6);
            one_over_lambda.reserve(dim_one_over_lambda_0__);
            for (size_t k_0__ = 0; k_0__ < dim_one_over_lambda_0__; ++k_0__) {
                if (jacobian__)
                    one_over_lambda.push_back(in__.scalar_lb_constrain(0,lp__));
                else
                    one_over_lambda.push_back(in__.scalar_lb_constrain(0));
            }

            Eigen::Matrix<T__,Eigen::Dynamic,1>  z_b;
            (void) z_b;  // dummy to suppress unused var warning
            if (jacobian__)
                z_b = in__.vector_constrain(q,lp__);
            else
                z_b = in__.vector_constrain(q);

            Eigen::Matrix<T__,Eigen::Dynamic,1>  z_T;
            (void) z_T;  // dummy to suppress unused var warning
            if (jacobian__)
                z_T = in__.vector_constrain(len_z_T,lp__);
            else
                z_T = in__.vector_constrain(len_z_T);

            Eigen::Matrix<T__,Eigen::Dynamic,1>  rho;
            (void) rho;  // dummy to suppress unused var warning
            if (jacobian__)
                rho = in__.vector_lub_constrain(0,1,len_rho,lp__);
            else
                rho = in__.vector_lub_constrain(0,1,len_rho);

            Eigen::Matrix<T__,Eigen::Dynamic,1>  zeta;
            (void) zeta;  // dummy to suppress unused var warning
            if (jacobian__)
                zeta = in__.vector_lb_constrain(0,len_concentration,lp__);
            else
                zeta = in__.vector_lb_constrain(0,len_concentration);

            Eigen::Matrix<T__,Eigen::Dynamic,1>  tau;
            (void) tau;  // dummy to suppress unused var warning
            if (jacobian__)
                tau = in__.vector_lb_constrain(0,t,lp__);
            else
                tau = in__.vector_lb_constrain(0,t);

            vector<T__> aux_unscaled;
            size_t dim_aux_unscaled_0__ = logical_gt(family,6);
            aux_unscaled.reserve(dim_aux_unscaled_0__);
            for (size_t k_0__ = 0; k_0__ < dim_aux_unscaled_0__; ++k_0__) {
                if (jacobian__)
                    aux_unscaled.push_back(in__.scalar_lb_constrain(0,lp__));
                else
                    aux_unscaled.push_back(in__.scalar_lb_constrain(0));
            }

            vector<Eigen::Matrix<T__,Eigen::Dynamic,1> > noise;
            size_t dim_noise_0__ = logical_eq(family,8);
            noise.reserve(dim_noise_0__);
            for (size_t k_0__ = 0; k_0__ < dim_noise_0__; ++k_0__) {
                if (jacobian__)
                    noise.push_back(in__.vector_lb_constrain(0,N,lp__));
                else
                    noise.push_back(in__.vector_lb_constrain(0,N));
            }


            // transformed parameters
            current_statement_begin__ = 514;
            T__ aux;
            (void) aux;  // dummy to suppress unused var warning

            stan::math::initialize(aux, DUMMY_VAR__);
            stan::math::fill(aux,DUMMY_VAR__);
            stan::math::assign(aux,stan::math::negative_infinity());
            current_statement_begin__ = 516;
            validate_non_negative_index("beta", "K", K);
            Eigen::Matrix<T__,Eigen::Dynamic,1>  beta(static_cast<Eigen::VectorXd::Index>(K));
            (void) beta;  // dummy to suppress unused var warning

            stan::math::initialize(beta, DUMMY_VAR__);
            stan::math::fill(beta,DUMMY_VAR__);
            current_statement_begin__ = 517;
            validate_non_negative_index("beta_smooth", "K_smooth", K_smooth);
            Eigen::Matrix<T__,Eigen::Dynamic,1>  beta_smooth(static_cast<Eigen::VectorXd::Index>(K_smooth));
            (void) beta_smooth;  // dummy to suppress unused var warning

            stan::math::initialize(beta_smooth, DUMMY_VAR__);
            stan::math::fill(beta_smooth,DUMMY_VAR__);
            current_statement_begin__ = 518;
            validate_non_negative_index("smooth_sd", "(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,\"smooth_map\",1) : 0 )", (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
            Eigen::Matrix<T__,Eigen::Dynamic,1>  smooth_sd(static_cast<Eigen::VectorXd::Index>((logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 )));
            (void) smooth_sd;  // dummy to suppress unused var warning

            stan::math::initialize(smooth_sd, DUMMY_VAR__);
            stan::math::fill(smooth_sd,DUMMY_VAR__);
            current_statement_begin__ = 519;
            validate_non_negative_index("b", "q", q);
            Eigen::Matrix<T__,Eigen::Dynamic,1>  b(static_cast<Eigen::VectorXd::Index>(q));
            (void) b;  // dummy to suppress unused var warning

            stan::math::initialize(b, DUMMY_VAR__);
            stan::math::fill(b,DUMMY_VAR__);
            current_statement_begin__ = 520;
            validate_non_negative_index("theta_L", "len_theta_L", len_theta_L);
            Eigen::Matrix<T__,Eigen::Dynamic,1>  theta_L(static_cast<Eigen::VectorXd::Index>(len_theta_L));
            (void) theta_L;  // dummy to suppress unused var warning

            stan::math::initialize(theta_L, DUMMY_VAR__);
            stan::math::fill(theta_L,DUMMY_VAR__);


            current_statement_begin__ = 521;
            if (as_bool(logical_eq(prior_dist,0))) {
                current_statement_begin__ = 521;
                stan::math::assign(beta, z_beta);
            } else if (as_bool(logical_eq(prior_dist,1))) {
                current_statement_begin__ = 522;
                stan::math::assign(beta, add(elt_multiply(z_beta,prior_scale),prior_mean));
            } else if (as_bool(logical_eq(prior_dist,2))) {
                current_statement_begin__ = 523;
                for (int k = 1; k <= K; ++k) {

                    current_statement_begin__ = 524;
                    stan::math::assign(get_base1_lhs(beta,k,"beta",1), ((CFt(get_base1(z_beta,k,"z_beta",1),get_base1(prior_df,k,"prior_df",1), pstream__) * get_base1(prior_scale,k,"prior_scale",1)) + get_base1(prior_mean,k,"prior_mean",1)));
                }
            } else if (as_bool(logical_eq(prior_dist,3))) {
                {
                current_statement_begin__ = 527;
                T__ c2;
                (void) c2;  // dummy to suppress unused var warning

                stan::math::initialize(c2, DUMMY_VAR__);
                stan::math::fill(c2,DUMMY_VAR__);
                stan::math::assign(c2,(square(slab_scale) * get_base1(caux,1,"caux",1)));


                current_statement_begin__ = 528;
                if (as_bool((primitive_value(logical_eq(is_continuous,1)) && primitive_value(logical_eq(family,1))))) {
                    current_statement_begin__ = 529;
                    stan::math::assign(beta, hs_prior(z_beta,global,local,global_prior_scale,aux,c2, pstream__));
                } else {
                    current_statement_begin__ = 530;
                    stan::math::assign(beta, hs_prior(z_beta,global,local,global_prior_scale,1,c2, pstream__));
                }
                }
            } else if (as_bool(logical_eq(prior_dist,4))) {
                {
                current_statement_begin__ = 533;
                T__ c2;
                (void) c2;  // dummy to suppress unused var warning

                stan::math::initialize(c2, DUMMY_VAR__);
                stan::math::fill(c2,DUMMY_VAR__);
                stan::math::assign(c2,(square(slab_scale) * get_base1(caux,1,"caux",1)));


                current_statement_begin__ = 534;
                if (as_bool((primitive_value(logical_eq(is_continuous,1)) && primitive_value(logical_eq(family,1))))) {
                    current_statement_begin__ = 535;
                    stan::math::assign(beta, hsplus_prior(z_beta,global,local,global_prior_scale,aux,c2, pstream__));
                } else {
                    current_statement_begin__ = 536;
                    stan::math::assign(beta, hsplus_prior(z_beta,global,local,global_prior_scale,1,c2, pstream__));
                }
                }
            } else if (as_bool(logical_eq(prior_dist,5))) {
                current_statement_begin__ = 539;
                stan::math::assign(beta, add(prior_mean,elt_multiply(elt_multiply(prior_scale,sqrt(multiply(2,get_base1(mix,1,"mix",1)))),z_beta)));
            } else if (as_bool(logical_eq(prior_dist,6))) {
                current_statement_begin__ = 541;
                stan::math::assign(beta, add(prior_mean,elt_multiply(elt_multiply(multiply(get_base1(one_over_lambda,1,"one_over_lambda",1),prior_scale),sqrt(multiply(2,get_base1(mix,1,"mix",1)))),z_beta)));
            } else if (as_bool(logical_eq(prior_dist,7))) {
                {
                current_statement_begin__ = 543;
                int z_pos(0);
                (void) z_pos;  // dummy to suppress unused var warning

                stan::math::fill(z_pos, std::numeric_limits<int>::min());
                stan::math::assign(z_pos,1);


                current_statement_begin__ = 544;
                for (int k = 1; k <= K; ++k) {

                    current_statement_begin__ = 545;
                    stan::math::assign(get_base1_lhs(beta,k,"beta",1), get_base1(z_beta,z_pos,"z_beta",1));
                    current_statement_begin__ = 546;
                    stan::math::assign(z_pos, (z_pos + 1));
                    current_statement_begin__ = 547;
                    for (int n = 2; n <= get_base1(num_normals,k,"num_normals",1); ++n) {

                        current_statement_begin__ = 548;
                        stan::math::assign(get_base1_lhs(beta,k,"beta",1), (get_base1(beta,k,"beta",1) * get_base1(z_beta,z_pos,"z_beta",1)));
                        current_statement_begin__ = 549;
                        stan::math::assign(z_pos, (z_pos + 1));
                    }
                    current_statement_begin__ = 551;
                    stan::math::assign(get_base1_lhs(beta,k,"beta",1), ((get_base1(beta,k,"beta",1) * pow(get_base1(prior_scale,k,"prior_scale",1),get_base1(num_normals,k,"num_normals",1))) + get_base1(prior_mean,k,"prior_mean",1)));
                }
                }
            }
            current_statement_begin__ = 555;
            if (as_bool(K_smooth)) {

                current_statement_begin__ = 556;
                stan::math::assign(smooth_sd, add(prior_mean_for_smooth,elt_multiply(prior_scale_for_smooth,smooth_sd_raw)));
                current_statement_begin__ = 557;
                if (as_bool((primitive_value(is_continuous) && primitive_value(logical_eq(family,1))))) {
                    current_statement_begin__ = 557;
                    stan::math::assign(smooth_sd, multiply(smooth_sd,aux));
                }
                current_statement_begin__ = 558;
                stan::math::assign(beta_smooth, elt_multiply(z_beta_smooth,stan::model::rvalue(smooth_sd, stan::model::cons_list(stan::model::index_multi(smooth_map), stan::model::nil_index_list()), "smooth_sd")));
            }
            current_statement_begin__ = 561;
            if (as_bool((primitive_value(logical_gt(family,6)) && primitive_value((primitive_value(logical_eq(prior_dist_for_aux,0)) || primitive_value(logical_lte(prior_scale_for_aux,0))))))) {
                current_statement_begin__ = 562;
                stan::math::assign(aux, get_base1(aux_unscaled,1,"aux_unscaled",1));
            } else if (as_bool(logical_gt(family,6))) {

                current_statement_begin__ = 564;
                stan::math::assign(aux, (prior_scale_for_aux * get_base1(aux_unscaled,1,"aux_unscaled",1)));
                current_statement_begin__ = 565;
                if (as_bool(logical_lte(prior_dist_for_aux,2))) {
                    current_statement_begin__ = 566;
                    stan::math::assign(aux, (aux + prior_mean_for_aux));
                }
            }
            current_statement_begin__ = 569;
            if (as_bool(logical_gt(t,0))) {

                current_statement_begin__ = 570;
                if (as_bool(logical_eq(special_case,1))) {
                    {
                    current_statement_begin__ = 571;
                    int start(0);
                    (void) start;  // dummy to suppress unused var warning

                    stan::math::fill(start, std::numeric_limits<int>::min());
                    stan::math::assign(start,1);


                    current_statement_begin__ = 572;
                    stan::math::assign(theta_L, elt_multiply(scale,(logical_eq(family,6) ? stan::math::promote_scalar<T__>(tau) : stan::math::promote_scalar<T__>(multiply(tau,aux)) )));
                    current_statement_begin__ = 573;
                    if (as_bool(logical_eq(t,1))) {
                        current_statement_begin__ = 573;
                        stan::math::assign(b, multiply(get_base1(theta_L,1,"theta_L",1),z_b));
                    } else {
                        current_statement_begin__ = 574;
                        for (int i = 1; i <= t; ++i) {
                            {
                            current_statement_begin__ = 575;
                            int end(0);
                            (void) end;  // dummy to suppress unused var warning

                            stan::math::fill(end, std::numeric_limits<int>::min());
                            stan::math::assign(end,((start + get_base1(l,i,"l",1)) - 1));


                            current_statement_begin__ = 576;
                            stan::model::assign(b, 
                                        stan::model::cons_list(stan::model::index_min_max(start, end), stan::model::nil_index_list()), 
                                        multiply(get_base1(theta_L,i,"theta_L",1),stan::model::rvalue(z_b, stan::model::cons_list(stan::model::index_min_max(start, end), stan::model::nil_index_list()), "z_b")), 
                                        "assigning variable b");
                            current_statement_begin__ = 577;
                            stan::math::assign(start, (end + 1));
                            }
                        }
                    }
                    }
                } else {

                    current_statement_begin__ = 581;
                    if (as_bool(logical_eq(family,6))) {
                        current_statement_begin__ = 582;
                        stan::math::assign(theta_L, make_theta_L(len_theta_L,p,1.0,tau,scale,zeta,rho,z_T, pstream__));
                    } else {
                        current_statement_begin__ = 585;
                        stan::math::assign(theta_L, make_theta_L(len_theta_L,p,aux,tau,scale,zeta,rho,z_T, pstream__));
                    }
                    current_statement_begin__ = 587;
                    stan::math::assign(b, make_b(z_b,theta_L,p,l, pstream__));
                }
            }

            // validate transformed parameters
            if (stan::math::is_uninitialized(aux)) {
                std::stringstream msg__;
                msg__ << "Undefined transformed parameter: aux";
                throw std::runtime_error(msg__.str());
            }
            for (int i0__ = 0; i0__ < K; ++i0__) {
                if (stan::math::is_uninitialized(beta(i0__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: beta" << '[' << i0__ << ']';
                    throw std::runtime_error(msg__.str());
                }
            }
            for (int i0__ = 0; i0__ < K_smooth; ++i0__) {
                if (stan::math::is_uninitialized(beta_smooth(i0__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: beta_smooth" << '[' << i0__ << ']';
                    throw std::runtime_error(msg__.str());
                }
            }
            for (int i0__ = 0; i0__ < (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++i0__) {
                if (stan::math::is_uninitialized(smooth_sd(i0__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: smooth_sd" << '[' << i0__ << ']';
                    throw std::runtime_error(msg__.str());
                }
            }
            for (int i0__ = 0; i0__ < q; ++i0__) {
                if (stan::math::is_uninitialized(b(i0__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: b" << '[' << i0__ << ']';
                    throw std::runtime_error(msg__.str());
                }
            }
            for (int i0__ = 0; i0__ < len_theta_L; ++i0__) {
                if (stan::math::is_uninitialized(theta_L(i0__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: theta_L" << '[' << i0__ << ']';
                    throw std::runtime_error(msg__.str());
                }
            }

            const char* function__ = "validate transformed params";
            (void) function__;  // dummy to suppress unused var warning
            current_statement_begin__ = 514;
            current_statement_begin__ = 516;
            current_statement_begin__ = 517;
            current_statement_begin__ = 518;
            current_statement_begin__ = 519;
            current_statement_begin__ = 520;

            // model body
            {
            current_statement_begin__ = 592;
            validate_non_negative_index("eta", "N", N);
            Eigen::Matrix<T__,Eigen::Dynamic,1>  eta(static_cast<Eigen::VectorXd::Index>(N));
            (void) eta;  // dummy to suppress unused var warning

            stan::math::initialize(eta, DUMMY_VAR__);
            stan::math::fill(eta,DUMMY_VAR__);


            current_statement_begin__ = 593;
            if (as_bool(logical_gt(K,0))) {

                current_statement_begin__ = 594;
                if (as_bool(dense_X)) {
                    current_statement_begin__ = 594;
                    stan::math::assign(eta, multiply(get_base1(X,1,"X",1),beta));
                } else {
                    current_statement_begin__ = 595;
                    stan::math::assign(eta, csr_matrix_times_vector2(N,K,w_X,v_X,u_X,beta, pstream__));
                }
            } else {
                current_statement_begin__ = 597;
                stan::math::assign(eta, rep_vector(0.0,N));
            }
            current_statement_begin__ = 598;
            if (as_bool(logical_eq(has_offset,1))) {
                current_statement_begin__ = 598;
                stan::math::assign(eta, add(eta,offset));
            }
            current_statement_begin__ = 599;
            if (as_bool(K_smooth)) {
                current_statement_begin__ = 599;
                stan::math::assign(eta, add(eta,multiply(S,beta_smooth)));
            }
            current_statement_begin__ = 600;
            if (as_bool(logical_gt(t,0))) {

                current_statement_begin__ = 601;
                if (as_bool(special_case)) {
                    current_statement_begin__ = 601;
                    for (int i = 1; i <= t; ++i) {
                        current_statement_begin__ = 601;
                        stan::math::assign(eta, add(eta,stan::model::rvalue(b, stan::model::cons_list(stan::model::index_multi(get_base1(V,i,"V",1)), stan::model::nil_index_list()), "b")));
                    }
                } else {
                    current_statement_begin__ = 602;
                    stan::math::assign(eta, add(eta,csr_matrix_times_vector2(N,q,w,v,u,b, pstream__)));
                }
            }
            current_statement_begin__ = 604;
            if (as_bool(logical_eq(has_intercept,1))) {

                current_statement_begin__ = 605;
                if (as_bool(logical_eq(link,1))) {
                    current_statement_begin__ = 605;
                    stan::math::assign(eta, add(eta,get_base1(gamma,1,"gamma",1)));
                } else {
                    current_statement_begin__ = 606;
                    stan::math::assign(eta, add(subtract(eta,min(eta)),get_base1(gamma,1,"gamma",1)));
                }
            } else {

                current_statement_begin__ = 610;
                stan::math::assign(eta, add(eta,dot_product(xbar,beta)));
            }
            current_statement_begin__ = 613;
            if (as_bool(logical_eq(family,8))) {

                current_statement_begin__ = 614;
                if (as_bool(logical_eq(link,1))) {
                    current_statement_begin__ = 614;
                    stan::math::assign(eta, add(add(eta,log(aux)),log(get_base1(noise,1,"noise",1))));
                } else if (as_bool(logical_eq(link,2))) {
                    current_statement_begin__ = 615;
                    stan::math::assign(eta, elt_multiply(multiply(eta,aux),get_base1(noise,1,"noise",1)));
                } else {
                    current_statement_begin__ = 616;
                    stan::math::assign(eta, add(add(eta,sqrt(aux)),sqrt(get_base1(noise,1,"noise",1))));
                }
            }
            current_statement_begin__ = 620;
            if (as_bool((primitive_value(logical_eq(has_weights,0)) && primitive_value(logical_eq(prior_PD,0))))) {

                current_statement_begin__ = 621;
                if (as_bool(logical_neq(family,7))) {

                    current_statement_begin__ = 622;
                    if (as_bool(logical_eq(link,1))) {
                        current_statement_begin__ = 622;
                        lp_accum__.add(poisson_log_log(y,eta));
                    } else {
                        current_statement_begin__ = 623;
                        lp_accum__.add(poisson_log(y,linkinv_count(eta,link, pstream__)));
                    }
                } else {

                    current_statement_begin__ = 626;
                    if (as_bool(logical_eq(link,1))) {
                        current_statement_begin__ = 626;
                        lp_accum__.add(neg_binomial_2_log_log(y,eta,aux));
                    } else {
                        current_statement_begin__ = 627;
                        lp_accum__.add(neg_binomial_2_log(y,linkinv_count(eta,link, pstream__),aux));
                    }
                }
            } else if (as_bool((primitive_value(logical_neq(family,7)) && primitive_value(logical_eq(prior_PD,0))))) {
                current_statement_begin__ = 631;
                lp_accum__.add(dot_product(weights,pw_pois(y,eta,link, pstream__)));
            } else if (as_bool(logical_eq(prior_PD,0))) {
                current_statement_begin__ = 633;
                lp_accum__.add(dot_product(weights,pw_nb(y,eta,aux,link, pstream__)));
            }
            current_statement_begin__ = 636;
            if (as_bool((primitive_value((primitive_value(logical_gt(family,6)) && primitive_value(logical_gt(prior_dist_for_aux,0)))) && primitive_value(logical_gt(prior_scale_for_aux,0))))) {
                {
                current_statement_begin__ = 638;
                T__ log_half;
                (void) log_half;  // dummy to suppress unused var warning

                stan::math::initialize(log_half, DUMMY_VAR__);
                stan::math::fill(log_half,DUMMY_VAR__);
                stan::math::assign(log_half,-(0.6931471805599454));


                current_statement_begin__ = 639;
                if (as_bool(logical_eq(prior_dist_for_aux,1))) {
                    current_statement_begin__ = 640;
                    lp_accum__.add((normal_log(aux_unscaled,0,1) - log_half));
                } else if (as_bool(logical_eq(prior_dist_for_aux,2))) {
                    current_statement_begin__ = 642;
                    lp_accum__.add((student_t_log(aux_unscaled,prior_df_for_aux,0,1) - log_half));
                } else {
                    current_statement_begin__ = 644;
                    lp_accum__.add(exponential_log(aux_unscaled,1));
                }
                }
            }
            current_statement_begin__ = 648;
            if (as_bool(logical_eq(prior_dist,1))) {
                current_statement_begin__ = 648;
                lp_accum__.add(normal_log(z_beta,0,1));
            } else if (as_bool(logical_eq(prior_dist,2))) {
                current_statement_begin__ = 649;
                lp_accum__.add(normal_log(z_beta,0,1));
            } else if (as_bool(logical_eq(prior_dist,3))) {
                {
                current_statement_begin__ = 651;
                T__ log_half;
                (void) log_half;  // dummy to suppress unused var warning

                stan::math::initialize(log_half, DUMMY_VAR__);
                stan::math::fill(log_half,DUMMY_VAR__);
                stan::math::assign(log_half,-(0.6931471805599454));


                current_statement_begin__ = 652;
                lp_accum__.add(normal_log(z_beta,0,1));
                current_statement_begin__ = 653;
                lp_accum__.add((normal_log(get_base1(local,1,"local",1),0,1) - log_half));
                current_statement_begin__ = 654;
                lp_accum__.add(inv_gamma_log(get_base1(local,2,"local",1),multiply(0.5,prior_df),multiply(0.5,prior_df)));
                current_statement_begin__ = 655;
                lp_accum__.add((normal_log(get_base1(global,1,"global",1),0,1) - log_half));
                current_statement_begin__ = 656;
                lp_accum__.add(inv_gamma_log(get_base1(global,2,"global",1),(0.5 * global_prior_df),(0.5 * global_prior_df)));
                current_statement_begin__ = 657;
                lp_accum__.add(inv_gamma_log(caux,(0.5 * slab_df),(0.5 * slab_df)));
                }
            } else if (as_bool(logical_eq(prior_dist,4))) {
                {
                current_statement_begin__ = 660;
                T__ log_half;
                (void) log_half;  // dummy to suppress unused var warning

                stan::math::initialize(log_half, DUMMY_VAR__);
                stan::math::fill(log_half,DUMMY_VAR__);
                stan::math::assign(log_half,-(0.6931471805599454));


                current_statement_begin__ = 661;
                lp_accum__.add(normal_log(z_beta,0,1));
                current_statement_begin__ = 662;
                lp_accum__.add((normal_log(get_base1(local,1,"local",1),0,1) - log_half));
                current_statement_begin__ = 663;
                lp_accum__.add(inv_gamma_log(get_base1(local,2,"local",1),multiply(0.5,prior_df),multiply(0.5,prior_df)));
                current_statement_begin__ = 664;
                lp_accum__.add((normal_log(get_base1(local,3,"local",1),0,1) - log_half));
                current_statement_begin__ = 666;
                lp_accum__.add(inv_gamma_log(get_base1(local,4,"local",1),multiply(0.5,prior_scale),multiply(0.5,prior_scale)));
                current_statement_begin__ = 667;
                lp_accum__.add((normal_log(get_base1(global,1,"global",1),0,1) - log_half));
                current_statement_begin__ = 668;
                lp_accum__.add(inv_gamma_log(get_base1(global,2,"global",1),(0.5 * global_prior_df),(0.5 * global_prior_df)));
                current_statement_begin__ = 669;
                lp_accum__.add(inv_gamma_log(caux,(0.5 * slab_df),(0.5 * slab_df)));
                }
            } else if (as_bool(logical_eq(prior_dist,5))) {

                current_statement_begin__ = 672;
                lp_accum__.add(normal_log(z_beta,0,1));
                current_statement_begin__ = 673;
                lp_accum__.add(exponential_log(get_base1(mix,1,"mix",1),1));
            } else if (as_bool(logical_eq(prior_dist,6))) {

                current_statement_begin__ = 676;
                lp_accum__.add(normal_log(z_beta,0,1));
                current_statement_begin__ = 677;
                lp_accum__.add(exponential_log(get_base1(mix,1,"mix",1),1));
                current_statement_begin__ = 678;
                lp_accum__.add(chi_square_log(get_base1(one_over_lambda,1,"one_over_lambda",1),get_base1(prior_df,1,"prior_df",1)));
            } else if (as_bool(logical_eq(prior_dist,7))) {

                current_statement_begin__ = 681;
                lp_accum__.add(normal_log(z_beta,0,1));
            }
            current_statement_begin__ = 686;
            if (as_bool(logical_eq(has_intercept,1))) {

                current_statement_begin__ = 687;
                if (as_bool(logical_eq(prior_dist_for_intercept,1))) {
                    current_statement_begin__ = 688;
                    lp_accum__.add(normal_log(gamma,prior_mean_for_intercept,prior_scale_for_intercept));
                } else if (as_bool(logical_eq(prior_dist_for_intercept,2))) {
                    current_statement_begin__ = 690;
                    lp_accum__.add(student_t_log(gamma,prior_df_for_intercept,prior_mean_for_intercept,prior_scale_for_intercept));
                }
            }
            current_statement_begin__ = 695;
            if (as_bool(K_smooth)) {

                current_statement_begin__ = 696;
                lp_accum__.add(normal_log(z_beta_smooth,0,1));
                current_statement_begin__ = 697;
                if (as_bool(logical_gt(prior_dist_for_smooth,0))) {
                    {
                    current_statement_begin__ = 698;
                    T__ log_half;
                    (void) log_half;  // dummy to suppress unused var warning

                    stan::math::initialize(log_half, DUMMY_VAR__);
                    stan::math::fill(log_half,DUMMY_VAR__);
                    stan::math::assign(log_half,-(0.6931471805599454));


                    current_statement_begin__ = 699;
                    if (as_bool(logical_eq(prior_dist_for_smooth,1))) {
                        current_statement_begin__ = 700;
                        lp_accum__.add((normal_log(smooth_sd_raw,0,1) - log_half));
                    } else if (as_bool(logical_eq(prior_dist_for_smooth,2))) {
                        current_statement_begin__ = 702;
                        lp_accum__.add((student_t_log(smooth_sd_raw,prior_df_for_smooth,0,1) - log_half));
                    } else if (as_bool(logical_eq(prior_dist_for_smooth,3))) {
                        current_statement_begin__ = 704;
                        lp_accum__.add(exponential_log(smooth_sd_raw,1));
                    }
                    }
                }
            }
            current_statement_begin__ = 709;
            if (as_bool(logical_eq(family,8))) {
                current_statement_begin__ = 709;
                lp_accum__.add(gamma_log(get_base1(noise,1,"noise",1),aux,1));
            }
            current_statement_begin__ = 711;
            if (as_bool(logical_gt(t,0))) {
                current_statement_begin__ = 711;
                decov_lp(z_b,z_T,rho,zeta,tau,regularization,delta,shape,t,p, lp__, lp_accum__, pstream__);
            }
            }

        } 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("gamma");
        names__.push_back("z_beta");
        names__.push_back("z_beta_smooth");
        names__.push_back("smooth_sd_raw");
        names__.push_back("global");
        names__.push_back("local");
        names__.push_back("caux");
        names__.push_back("mix");
        names__.push_back("one_over_lambda");
        names__.push_back("z_b");
        names__.push_back("z_T");
        names__.push_back("rho");
        names__.push_back("zeta");
        names__.push_back("tau");
        names__.push_back("aux_unscaled");
        names__.push_back("noise");
        names__.push_back("aux");
        names__.push_back("beta");
        names__.push_back("beta_smooth");
        names__.push_back("smooth_sd");
        names__.push_back("b");
        names__.push_back("theta_L");
        names__.push_back("alpha");
        names__.push_back("mean_PPD");
    }


    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(has_intercept);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_eq(prior_dist,7) ? sum(num_normals) : K ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(K_smooth);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(hs);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(hs);
        dims__.push_back(K);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(logical_gt(hs,0));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))));
        dims__.push_back(K);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(logical_eq(prior_dist,6));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(q);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(len_z_T);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(len_rho);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(len_concentration);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(t);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(logical_gt(family,6));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(logical_eq(family,8));
        dims__.push_back(N);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(K);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(K_smooth);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(q);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(len_theta_L);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(has_intercept);
        dimss__.push_back(dims__);
        dims__.resize(0);
        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 {
        vars__.resize(0);
        stan::io::reader<double> in__(params_r__,params_i__);
        static const char* function__ = "model_count_namespace::write_array";
        (void) function__;  // dummy to suppress unused var warning
        // read-transform, write parameters
        vector<double> gamma;
        size_t dim_gamma_0__ = has_intercept;
        for (size_t k_0__ = 0; k_0__ < dim_gamma_0__; ++k_0__) {
            gamma.push_back(in__.scalar_lb_constrain((logical_eq(link,1) ? stan::math::negative_infinity() : 0.0 )));
        }
        vector_d z_beta = in__.vector_constrain((logical_eq(prior_dist,7) ? sum(num_normals) : K ));
        vector_d z_beta_smooth = in__.vector_constrain(K_smooth);
        vector_d smooth_sd_raw = in__.vector_lb_constrain(0,(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
        vector<double> global;
        size_t dim_global_0__ = hs;
        for (size_t k_0__ = 0; k_0__ < dim_global_0__; ++k_0__) {
            global.push_back(in__.scalar_lb_constrain(0));
        }
        vector<vector_d> local;
        size_t dim_local_0__ = hs;
        for (size_t k_0__ = 0; k_0__ < dim_local_0__; ++k_0__) {
            local.push_back(in__.vector_lb_constrain(0,K));
        }
        vector<double> caux;
        size_t dim_caux_0__ = logical_gt(hs,0);
        for (size_t k_0__ = 0; k_0__ < dim_caux_0__; ++k_0__) {
            caux.push_back(in__.scalar_lb_constrain(0));
        }
        vector<vector_d> mix;
        size_t dim_mix_0__ = (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6)));
        for (size_t k_0__ = 0; k_0__ < dim_mix_0__; ++k_0__) {
            mix.push_back(in__.vector_lb_constrain(0,K));
        }
        vector<double> one_over_lambda;
        size_t dim_one_over_lambda_0__ = logical_eq(prior_dist,6);
        for (size_t k_0__ = 0; k_0__ < dim_one_over_lambda_0__; ++k_0__) {
            one_over_lambda.push_back(in__.scalar_lb_constrain(0));
        }
        vector_d z_b = in__.vector_constrain(q);
        vector_d z_T = in__.vector_constrain(len_z_T);
        vector_d rho = in__.vector_lub_constrain(0,1,len_rho);
        vector_d zeta = in__.vector_lb_constrain(0,len_concentration);
        vector_d tau = in__.vector_lb_constrain(0,t);
        vector<double> aux_unscaled;
        size_t dim_aux_unscaled_0__ = logical_gt(family,6);
        for (size_t k_0__ = 0; k_0__ < dim_aux_unscaled_0__; ++k_0__) {
            aux_unscaled.push_back(in__.scalar_lb_constrain(0));
        }
        vector<vector_d> noise;
        size_t dim_noise_0__ = logical_eq(family,8);
        for (size_t k_0__ = 0; k_0__ < dim_noise_0__; ++k_0__) {
            noise.push_back(in__.vector_lb_constrain(0,N));
        }
            for (int k_0__ = 0; k_0__ < has_intercept; ++k_0__) {
            vars__.push_back(gamma[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < (logical_eq(prior_dist,7) ? sum(num_normals) : K ); ++k_0__) {
            vars__.push_back(z_beta[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < K_smooth; ++k_0__) {
            vars__.push_back(z_beta_smooth[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++k_0__) {
            vars__.push_back(smooth_sd_raw[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < hs; ++k_0__) {
            vars__.push_back(global[k_0__]);
            }
            for (int k_1__ = 0; k_1__ < K; ++k_1__) {
                for (int k_0__ = 0; k_0__ < hs; ++k_0__) {
                vars__.push_back(local[k_0__][k_1__]);
                }
            }
            for (int k_0__ = 0; k_0__ < logical_gt(hs,0); ++k_0__) {
            vars__.push_back(caux[k_0__]);
            }
            for (int k_1__ = 0; k_1__ < K; ++k_1__) {
                for (int k_0__ = 0; k_0__ < (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))); ++k_0__) {
                vars__.push_back(mix[k_0__][k_1__]);
                }
            }
            for (int k_0__ = 0; k_0__ < logical_eq(prior_dist,6); ++k_0__) {
            vars__.push_back(one_over_lambda[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < q; ++k_0__) {
            vars__.push_back(z_b[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < len_z_T; ++k_0__) {
            vars__.push_back(z_T[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < len_rho; ++k_0__) {
            vars__.push_back(rho[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < len_concentration; ++k_0__) {
            vars__.push_back(zeta[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < t; ++k_0__) {
            vars__.push_back(tau[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < logical_gt(family,6); ++k_0__) {
            vars__.push_back(aux_unscaled[k_0__]);
            }
            for (int k_1__ = 0; k_1__ < N; ++k_1__) {
                for (int k_0__ = 0; k_0__ < logical_eq(family,8); ++k_0__) {
                vars__.push_back(noise[k_0__][k_1__]);
                }
            }

        if (!include_tparams__) return;
        // declare and define transformed parameters
        double lp__ = 0.0;
        (void) lp__;  // dummy to suppress unused var warning
        stan::math::accumulator<double> lp_accum__;

        double DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

        try {
            current_statement_begin__ = 514;
            double aux(0.0);
            (void) aux;  // dummy to suppress unused var warning

            stan::math::initialize(aux, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(aux,DUMMY_VAR__);
            stan::math::assign(aux,stan::math::negative_infinity());
            current_statement_begin__ = 516;
            validate_non_negative_index("beta", "K", K);
            vector_d beta(static_cast<Eigen::VectorXd::Index>(K));
            (void) beta;  // dummy to suppress unused var warning

            stan::math::initialize(beta, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(beta,DUMMY_VAR__);
            current_statement_begin__ = 517;
            validate_non_negative_index("beta_smooth", "K_smooth", K_smooth);
            vector_d beta_smooth(static_cast<Eigen::VectorXd::Index>(K_smooth));
            (void) beta_smooth;  // dummy to suppress unused var warning

            stan::math::initialize(beta_smooth, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(beta_smooth,DUMMY_VAR__);
            current_statement_begin__ = 518;
            validate_non_negative_index("smooth_sd", "(logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,\"smooth_map\",1) : 0 )", (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ));
            vector_d smooth_sd(static_cast<Eigen::VectorXd::Index>((logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 )));
            (void) smooth_sd;  // dummy to suppress unused var warning

            stan::math::initialize(smooth_sd, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(smooth_sd,DUMMY_VAR__);
            current_statement_begin__ = 519;
            validate_non_negative_index("b", "q", q);
            vector_d b(static_cast<Eigen::VectorXd::Index>(q));
            (void) b;  // dummy to suppress unused var warning

            stan::math::initialize(b, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(b,DUMMY_VAR__);
            current_statement_begin__ = 520;
            validate_non_negative_index("theta_L", "len_theta_L", len_theta_L);
            vector_d theta_L(static_cast<Eigen::VectorXd::Index>(len_theta_L));
            (void) theta_L;  // dummy to suppress unused var warning

            stan::math::initialize(theta_L, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(theta_L,DUMMY_VAR__);


            current_statement_begin__ = 521;
            if (as_bool(logical_eq(prior_dist,0))) {
                current_statement_begin__ = 521;
                stan::math::assign(beta, z_beta);
            } else if (as_bool(logical_eq(prior_dist,1))) {
                current_statement_begin__ = 522;
                stan::math::assign(beta, add(elt_multiply(z_beta,prior_scale),prior_mean));
            } else if (as_bool(logical_eq(prior_dist,2))) {
                current_statement_begin__ = 523;
                for (int k = 1; k <= K; ++k) {

                    current_statement_begin__ = 524;
                    stan::math::assign(get_base1_lhs(beta,k,"beta",1), ((CFt(get_base1(z_beta,k,"z_beta",1),get_base1(prior_df,k,"prior_df",1), pstream__) * get_base1(prior_scale,k,"prior_scale",1)) + get_base1(prior_mean,k,"prior_mean",1)));
                }
            } else if (as_bool(logical_eq(prior_dist,3))) {
                {
                current_statement_begin__ = 527;
                double c2(0.0);
                (void) c2;  // dummy to suppress unused var warning

                stan::math::initialize(c2, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(c2,DUMMY_VAR__);
                stan::math::assign(c2,(square(slab_scale) * get_base1(caux,1,"caux",1)));


                current_statement_begin__ = 528;
                if (as_bool((primitive_value(logical_eq(is_continuous,1)) && primitive_value(logical_eq(family,1))))) {
                    current_statement_begin__ = 529;
                    stan::math::assign(beta, hs_prior(z_beta,global,local,global_prior_scale,aux,c2, pstream__));
                } else {
                    current_statement_begin__ = 530;
                    stan::math::assign(beta, hs_prior(z_beta,global,local,global_prior_scale,1,c2, pstream__));
                }
                }
            } else if (as_bool(logical_eq(prior_dist,4))) {
                {
                current_statement_begin__ = 533;
                double c2(0.0);
                (void) c2;  // dummy to suppress unused var warning

                stan::math::initialize(c2, std::numeric_limits<double>::quiet_NaN());
                stan::math::fill(c2,DUMMY_VAR__);
                stan::math::assign(c2,(square(slab_scale) * get_base1(caux,1,"caux",1)));


                current_statement_begin__ = 534;
                if (as_bool((primitive_value(logical_eq(is_continuous,1)) && primitive_value(logical_eq(family,1))))) {
                    current_statement_begin__ = 535;
                    stan::math::assign(beta, hsplus_prior(z_beta,global,local,global_prior_scale,aux,c2, pstream__));
                } else {
                    current_statement_begin__ = 536;
                    stan::math::assign(beta, hsplus_prior(z_beta,global,local,global_prior_scale,1,c2, pstream__));
                }
                }
            } else if (as_bool(logical_eq(prior_dist,5))) {
                current_statement_begin__ = 539;
                stan::math::assign(beta, add(prior_mean,elt_multiply(elt_multiply(prior_scale,sqrt(multiply(2,get_base1(mix,1,"mix",1)))),z_beta)));
            } else if (as_bool(logical_eq(prior_dist,6))) {
                current_statement_begin__ = 541;
                stan::math::assign(beta, add(prior_mean,elt_multiply(elt_multiply(multiply(get_base1(one_over_lambda,1,"one_over_lambda",1),prior_scale),sqrt(multiply(2,get_base1(mix,1,"mix",1)))),z_beta)));
            } else if (as_bool(logical_eq(prior_dist,7))) {
                {
                current_statement_begin__ = 543;
                int z_pos(0);
                (void) z_pos;  // dummy to suppress unused var warning

                stan::math::fill(z_pos, std::numeric_limits<int>::min());
                stan::math::assign(z_pos,1);


                current_statement_begin__ = 544;
                for (int k = 1; k <= K; ++k) {

                    current_statement_begin__ = 545;
                    stan::math::assign(get_base1_lhs(beta,k,"beta",1), get_base1(z_beta,z_pos,"z_beta",1));
                    current_statement_begin__ = 546;
                    stan::math::assign(z_pos, (z_pos + 1));
                    current_statement_begin__ = 547;
                    for (int n = 2; n <= get_base1(num_normals,k,"num_normals",1); ++n) {

                        current_statement_begin__ = 548;
                        stan::math::assign(get_base1_lhs(beta,k,"beta",1), (get_base1(beta,k,"beta",1) * get_base1(z_beta,z_pos,"z_beta",1)));
                        current_statement_begin__ = 549;
                        stan::math::assign(z_pos, (z_pos + 1));
                    }
                    current_statement_begin__ = 551;
                    stan::math::assign(get_base1_lhs(beta,k,"beta",1), ((get_base1(beta,k,"beta",1) * pow(get_base1(prior_scale,k,"prior_scale",1),get_base1(num_normals,k,"num_normals",1))) + get_base1(prior_mean,k,"prior_mean",1)));
                }
                }
            }
            current_statement_begin__ = 555;
            if (as_bool(K_smooth)) {

                current_statement_begin__ = 556;
                stan::math::assign(smooth_sd, add(prior_mean_for_smooth,elt_multiply(prior_scale_for_smooth,smooth_sd_raw)));
                current_statement_begin__ = 557;
                if (as_bool((primitive_value(is_continuous) && primitive_value(logical_eq(family,1))))) {
                    current_statement_begin__ = 557;
                    stan::math::assign(smooth_sd, multiply(smooth_sd,aux));
                }
                current_statement_begin__ = 558;
                stan::math::assign(beta_smooth, elt_multiply(z_beta_smooth,stan::model::rvalue(smooth_sd, stan::model::cons_list(stan::model::index_multi(smooth_map), stan::model::nil_index_list()), "smooth_sd")));
            }
            current_statement_begin__ = 561;
            if (as_bool((primitive_value(logical_gt(family,6)) && primitive_value((primitive_value(logical_eq(prior_dist_for_aux,0)) || primitive_value(logical_lte(prior_scale_for_aux,0))))))) {
                current_statement_begin__ = 562;
                stan::math::assign(aux, get_base1(aux_unscaled,1,"aux_unscaled",1));
            } else if (as_bool(logical_gt(family,6))) {

                current_statement_begin__ = 564;
                stan::math::assign(aux, (prior_scale_for_aux * get_base1(aux_unscaled,1,"aux_unscaled",1)));
                current_statement_begin__ = 565;
                if (as_bool(logical_lte(prior_dist_for_aux,2))) {
                    current_statement_begin__ = 566;
                    stan::math::assign(aux, (aux + prior_mean_for_aux));
                }
            }
            current_statement_begin__ = 569;
            if (as_bool(logical_gt(t,0))) {

                current_statement_begin__ = 570;
                if (as_bool(logical_eq(special_case,1))) {
                    {
                    current_statement_begin__ = 571;
                    int start(0);
                    (void) start;  // dummy to suppress unused var warning

                    stan::math::fill(start, std::numeric_limits<int>::min());
                    stan::math::assign(start,1);


                    current_statement_begin__ = 572;
                    stan::math::assign(theta_L, elt_multiply(scale,(logical_eq(family,6) ? stan::math::promote_scalar<double>(tau) : stan::math::promote_scalar<double>(multiply(tau,aux)) )));
                    current_statement_begin__ = 573;
                    if (as_bool(logical_eq(t,1))) {
                        current_statement_begin__ = 573;
                        stan::math::assign(b, multiply(get_base1(theta_L,1,"theta_L",1),z_b));
                    } else {
                        current_statement_begin__ = 574;
                        for (int i = 1; i <= t; ++i) {
                            {
                            current_statement_begin__ = 575;
                            int end(0);
                            (void) end;  // dummy to suppress unused var warning

                            stan::math::fill(end, std::numeric_limits<int>::min());
                            stan::math::assign(end,((start + get_base1(l,i,"l",1)) - 1));


                            current_statement_begin__ = 576;
                            stan::model::assign(b, 
                                        stan::model::cons_list(stan::model::index_min_max(start, end), stan::model::nil_index_list()), 
                                        multiply(get_base1(theta_L,i,"theta_L",1),stan::model::rvalue(z_b, stan::model::cons_list(stan::model::index_min_max(start, end), stan::model::nil_index_list()), "z_b")), 
                                        "assigning variable b");
                            current_statement_begin__ = 577;
                            stan::math::assign(start, (end + 1));
                            }
                        }
                    }
                    }
                } else {

                    current_statement_begin__ = 581;
                    if (as_bool(logical_eq(family,6))) {
                        current_statement_begin__ = 582;
                        stan::math::assign(theta_L, make_theta_L(len_theta_L,p,1.0,tau,scale,zeta,rho,z_T, pstream__));
                    } else {
                        current_statement_begin__ = 585;
                        stan::math::assign(theta_L, make_theta_L(len_theta_L,p,aux,tau,scale,zeta,rho,z_T, pstream__));
                    }
                    current_statement_begin__ = 587;
                    stan::math::assign(b, make_b(z_b,theta_L,p,l, pstream__));
                }
            }

            // validate transformed parameters
            current_statement_begin__ = 514;
            current_statement_begin__ = 516;
            current_statement_begin__ = 517;
            current_statement_begin__ = 518;
            current_statement_begin__ = 519;
            current_statement_begin__ = 520;

            // write transformed parameters
        vars__.push_back(aux);
            for (int k_0__ = 0; k_0__ < K; ++k_0__) {
            vars__.push_back(beta[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < K_smooth; ++k_0__) {
            vars__.push_back(beta_smooth[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++k_0__) {
            vars__.push_back(smooth_sd[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < q; ++k_0__) {
            vars__.push_back(b[k_0__]);
            }
            for (int k_0__ = 0; k_0__ < len_theta_L; ++k_0__) {
            vars__.push_back(theta_L[k_0__]);
            }

            if (!include_gqs__) return;
            // declare and define generated quantities
            current_statement_begin__ = 715;
            validate_non_negative_index("alpha", "has_intercept", has_intercept);
            vector<double> alpha(has_intercept, 0.0);
            stan::math::initialize(alpha, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(alpha,DUMMY_VAR__);
            current_statement_begin__ = 716;
            double mean_PPD(0.0);
            (void) mean_PPD;  // dummy to suppress unused var warning

            stan::math::initialize(mean_PPD, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(mean_PPD,DUMMY_VAR__);
            stan::math::assign(mean_PPD,0);


            current_statement_begin__ = 717;
            if (as_bool(logical_eq(has_intercept,1))) {

                current_statement_begin__ = 718;
                if (as_bool(dense_X)) {
                    current_statement_begin__ = 718;
                    stan::math::assign(get_base1_lhs(alpha,1,"alpha",1), (get_base1(gamma,1,"gamma",1) - dot_product(xbar,beta)));
                } else {
                    current_statement_begin__ = 719;
                    stan::math::assign(get_base1_lhs(alpha,1,"alpha",1), get_base1(gamma,1,"gamma",1));
                }
            }
            {
            current_statement_begin__ = 722;
            validate_non_negative_index("nu", "N", N);
            vector_d nu(static_cast<Eigen::VectorXd::Index>(N));
            (void) nu;  // dummy to suppress unused var warning

            stan::math::initialize(nu, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(nu,DUMMY_VAR__);
            current_statement_begin__ = 723;
            validate_non_negative_index("eta", "N", N);
            vector_d eta(static_cast<Eigen::VectorXd::Index>(N));
            (void) eta;  // dummy to suppress unused var warning

            stan::math::initialize(eta, std::numeric_limits<double>::quiet_NaN());
            stan::math::fill(eta,DUMMY_VAR__);


            current_statement_begin__ = 724;
            if (as_bool(logical_gt(K,0))) {

                current_statement_begin__ = 725;
                if (as_bool(dense_X)) {
                    current_statement_begin__ = 725;
                    stan::math::assign(eta, multiply(get_base1(X,1,"X",1),beta));
                } else {
                    current_statement_begin__ = 726;
                    stan::math::assign(eta, csr_matrix_times_vector2(N,K,w_X,v_X,u_X,beta, pstream__));
                }
            } else {
                current_statement_begin__ = 728;
                stan::math::assign(eta, rep_vector(0.0,N));
            }
            current_statement_begin__ = 729;
            if (as_bool(logical_eq(has_offset,1))) {
                current_statement_begin__ = 729;
                stan::math::assign(eta, add(eta,offset));
            }
            current_statement_begin__ = 730;
            if (as_bool(K_smooth)) {
                current_statement_begin__ = 730;
                stan::math::assign(eta, add(eta,multiply(S,beta_smooth)));
            }
            current_statement_begin__ = 731;
            if (as_bool(logical_gt(t,0))) {

                current_statement_begin__ = 732;
                if (as_bool(special_case)) {
                    current_statement_begin__ = 732;
                    for (int i = 1; i <= t; ++i) {
                        current_statement_begin__ = 732;
                        stan::math::assign(eta, add(eta,stan::model::rvalue(b, stan::model::cons_list(stan::model::index_multi(get_base1(V,i,"V",1)), stan::model::nil_index_list()), "b")));
                    }
                } else {
                    current_statement_begin__ = 733;
                    stan::math::assign(eta, add(eta,csr_matrix_times_vector2(N,q,w,v,u,b, pstream__)));
                }
            }
            current_statement_begin__ = 735;
            if (as_bool(logical_eq(has_intercept,1))) {

                current_statement_begin__ = 736;
                if (as_bool(logical_eq(link,1))) {
                    current_statement_begin__ = 736;
                    stan::math::assign(eta, add(eta,get_base1(gamma,1,"gamma",1)));
                } else {
                    {
                    current_statement_begin__ = 738;
                    double shift(0.0);
                    (void) shift;  // dummy to suppress unused var warning

                    stan::math::initialize(shift, std::numeric_limits<double>::quiet_NaN());
                    stan::math::fill(shift,DUMMY_VAR__);
                    stan::math::assign(shift,min(eta));


                    current_statement_begin__ = 739;
                    stan::math::assign(eta, add(subtract(eta,shift),get_base1(gamma,1,"gamma",1)));
                    current_statement_begin__ = 740;
                    stan::math::assign(get_base1_lhs(alpha,1,"alpha",1), (get_base1(alpha,1,"alpha",1) - shift));
                    }
                }
            } else {

                current_statement_begin__ = 745;
                stan::math::assign(eta, add(eta,dot_product(xbar,beta)));
            }
            current_statement_begin__ = 748;
            if (as_bool(logical_eq(family,8))) {

                current_statement_begin__ = 749;
                if (as_bool(logical_eq(link,1))) {
                    current_statement_begin__ = 749;
                    stan::math::assign(eta, add(add(eta,log(aux)),log(get_base1(noise,1,"noise",1))));
                } else if (as_bool(logical_eq(link,2))) {
                    current_statement_begin__ = 750;
                    stan::math::assign(eta, elt_multiply(multiply(eta,aux),get_base1(noise,1,"noise",1)));
                } else {
                    current_statement_begin__ = 751;
                    stan::math::assign(eta, add(add(eta,sqrt(aux)),sqrt(get_base1(noise,1,"noise",1))));
                }
            }
            current_statement_begin__ = 753;
            stan::math::assign(nu, linkinv_count(eta,link, pstream__));
            current_statement_begin__ = 754;
            if (as_bool(logical_neq(family,7))) {
                current_statement_begin__ = 754;
                for (int n = 1; n <= N; ++n) {

                    current_statement_begin__ = 755;
                    if (as_bool(logical_lt(get_base1(nu,n,"nu",1),poisson_max))) {
                        current_statement_begin__ = 755;
                        stan::math::assign(mean_PPD, (mean_PPD + poisson_rng(get_base1(nu,n,"nu",1), base_rng__)));
                    } else {
                        current_statement_begin__ = 756;
                        stan::math::assign(mean_PPD, (mean_PPD + normal_rng(get_base1(nu,n,"nu",1),sqrt(get_base1(nu,n,"nu",1)), base_rng__)));
                    }
                }
            } else {
                current_statement_begin__ = 758;
                for (int n = 1; n <= N; ++n) {
                    {
                    current_statement_begin__ = 759;
                    double gamma_temp(0.0);
                    (void) gamma_temp;  // dummy to suppress unused var warning

                    stan::math::initialize(gamma_temp, std::numeric_limits<double>::quiet_NaN());
                    stan::math::fill(gamma_temp,DUMMY_VAR__);


                    current_statement_begin__ = 760;
                    if (as_bool(is_inf(aux))) {
                        current_statement_begin__ = 760;
                        stan::math::assign(gamma_temp, get_base1(nu,n,"nu",1));
                    } else {
                        current_statement_begin__ = 761;
                        stan::math::assign(gamma_temp, gamma_rng(aux,(aux / get_base1(nu,n,"nu",1)), base_rng__));
                    }
                    current_statement_begin__ = 762;
                    if (as_bool(logical_lt(gamma_temp,poisson_max))) {
                        current_statement_begin__ = 763;
                        stan::math::assign(mean_PPD, (mean_PPD + poisson_rng(gamma_temp, base_rng__)));
                    } else {
                        current_statement_begin__ = 764;
                        stan::math::assign(mean_PPD, (mean_PPD + normal_rng(gamma_temp,sqrt(gamma_temp), base_rng__)));
                    }
                    }
                }
            }
            current_statement_begin__ = 766;
            stan::math::assign(mean_PPD, (mean_PPD / N));
            }

            // validate generated quantities
            current_statement_begin__ = 715;
            current_statement_begin__ = 716;

            // write generated quantities
            for (int k_0__ = 0; k_0__ < has_intercept; ++k_0__) {
            vars__.push_back(alpha[k_0__]);
            }
        vars__.push_back(mean_PPD);

        } 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];
    }

    static std::string model_name() {
        return "model_count";
    }


    void constrained_param_names(std::vector<std::string>& param_names__,
                                 bool include_tparams__ = true,
                                 bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        for (int k_0__ = 1; k_0__ <= has_intercept; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gamma" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= (logical_eq(prior_dist,7) ? sum(num_normals) : K ); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_beta" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= K_smooth; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_beta_smooth" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "smooth_sd_raw" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= hs; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "global" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_1__ = 1; k_1__ <= K; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= hs; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "local" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_0__ = 1; k_0__ <= logical_gt(hs,0); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "caux" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_1__ = 1; k_1__ <= K; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))); ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "mix" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_0__ = 1; k_0__ <= logical_eq(prior_dist,6); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "one_over_lambda" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= q; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_b" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_z_T; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_T" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_rho; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "rho" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_concentration; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "zeta" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= t; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= logical_gt(family,6); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "aux_unscaled" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_1__ = 1; k_1__ <= N; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= logical_eq(family,8); ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "noise" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__ && !include_tparams__) return;
        param_name_stream__.str(std::string());
        param_name_stream__ << "aux";
        param_names__.push_back(param_name_stream__.str());
        for (int k_0__ = 1; k_0__ <= K; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "beta" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= K_smooth; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "beta_smooth" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "smooth_sd" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= q; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "b" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_theta_L; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "theta_L" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }

        if (!include_gqs__) return;
        for (int k_0__ = 1; k_0__ <= has_intercept; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "alpha" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        param_name_stream__.str(std::string());
        param_name_stream__ << "mean_PPD";
        param_names__.push_back(param_name_stream__.str());
    }


    void unconstrained_param_names(std::vector<std::string>& param_names__,
                                   bool include_tparams__ = true,
                                   bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        for (int k_0__ = 1; k_0__ <= has_intercept; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gamma" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= (logical_eq(prior_dist,7) ? sum(num_normals) : K ); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_beta" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= K_smooth; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_beta_smooth" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "smooth_sd_raw" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= hs; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "global" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_1__ = 1; k_1__ <= K; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= hs; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "local" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_0__ = 1; k_0__ <= logical_gt(hs,0); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "caux" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_1__ = 1; k_1__ <= K; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= (primitive_value(logical_eq(prior_dist,5)) || primitive_value(logical_eq(prior_dist,6))); ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "mix" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_0__ = 1; k_0__ <= logical_eq(prior_dist,6); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "one_over_lambda" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= q; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_b" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_z_T; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "z_T" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_rho; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "rho" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_concentration; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "zeta" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= t; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= logical_gt(family,6); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "aux_unscaled" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_1__ = 1; k_1__ <= N; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= logical_eq(family,8); ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "noise" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__ && !include_tparams__) return;
        param_name_stream__.str(std::string());
        param_name_stream__ << "aux";
        param_names__.push_back(param_name_stream__.str());
        for (int k_0__ = 1; k_0__ <= K; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "beta" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= K_smooth; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "beta_smooth" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= (logical_gt(K_smooth,0) ? get_base1(smooth_map,K_smooth,"smooth_map",1) : 0 ); ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "smooth_sd" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= q; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "b" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        for (int k_0__ = 1; k_0__ <= len_theta_L; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "theta_L" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }

        if (!include_gqs__) return;
        for (int k_0__ = 1; k_0__ <= has_intercept; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "alpha" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        param_name_stream__.str(std::string());
        param_name_stream__ << "mean_PPD";
        param_names__.push_back(param_name_stream__.str());
    }

}; // model

}

typedef model_count_namespace::model_count stan_model;