// Code generated by Stan version 2.15.0

#include <stan/model/model_header.hpp>

namespace model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c_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__;

class model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c : public prob_grad {
private:
    int R;
    int S;
    int A;
    int K;
    matrix_d B;
    vector_d std_schedule;
    vector_d agg_wts;
    vector_d death_frac;
    vector_d pihat1;
    vector_d pihat_all;
    vector_d a3;
    vector_d b3;
    matrix_d N;
    vector<vector<int> > D;
    matrix_d lambda_star;
public:
    model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        typedef boost::ecuyer1988 rng_t;
        rng_t base_rng(0);  // 0 seed default
        ctor_body(context__, base_rng, pstream__);
    }

    template <class RNG>
    model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c(stan::io::var_context& context__,
        RNG& base_rng__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, base_rng__, pstream__);
    }

    template <class RNG>
    void ctor_body(stan::io::var_context& context__,
                   RNG& base_rng__,
                   std::ostream* pstream__) {
        current_statement_begin__ = -1;

        static const char* function__ = "model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c_namespace::model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c";
        (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
        context__.validate_dims("data initialization", "R", "int", context__.to_vec());
        R = int(0);
        vals_i__ = context__.vals_i("R");
        pos__ = 0;
        R = vals_i__[pos__++];
        context__.validate_dims("data initialization", "S", "int", context__.to_vec());
        S = int(0);
        vals_i__ = context__.vals_i("S");
        pos__ = 0;
        S = vals_i__[pos__++];
        context__.validate_dims("data initialization", "A", "int", context__.to_vec());
        A = int(0);
        vals_i__ = context__.vals_i("A");
        pos__ = 0;
        A = vals_i__[pos__++];
        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__++];
        validate_non_negative_index("B", "A", A);
        validate_non_negative_index("B", "K", K);
        context__.validate_dims("data initialization", "B", "matrix_d", context__.to_vec(A,K));
        validate_non_negative_index("B", "A", A);
        validate_non_negative_index("B", "K", K);
        B = matrix_d(static_cast<Eigen::VectorXd::Index>(A),static_cast<Eigen::VectorXd::Index>(K));
        vals_r__ = context__.vals_r("B");
        pos__ = 0;
        size_t B_m_mat_lim__ = A;
        size_t B_n_mat_lim__ = K;
        for (size_t n_mat__ = 0; n_mat__ < B_n_mat_lim__; ++n_mat__) {
            for (size_t m_mat__ = 0; m_mat__ < B_m_mat_lim__; ++m_mat__) {
                B(m_mat__,n_mat__) = vals_r__[pos__++];
            }
        }
        validate_non_negative_index("std_schedule", "A", A);
        context__.validate_dims("data initialization", "std_schedule", "vector_d", context__.to_vec(A));
        validate_non_negative_index("std_schedule", "A", A);
        std_schedule = vector_d(static_cast<Eigen::VectorXd::Index>(A));
        vals_r__ = context__.vals_r("std_schedule");
        pos__ = 0;
        size_t std_schedule_i_vec_lim__ = A;
        for (size_t i_vec__ = 0; i_vec__ < std_schedule_i_vec_lim__; ++i_vec__) {
            std_schedule[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("agg_wts", "R", R);
        context__.validate_dims("data initialization", "agg_wts", "vector_d", context__.to_vec(R));
        validate_non_negative_index("agg_wts", "R", R);
        agg_wts = vector_d(static_cast<Eigen::VectorXd::Index>(R));
        vals_r__ = context__.vals_r("agg_wts");
        pos__ = 0;
        size_t agg_wts_i_vec_lim__ = R;
        for (size_t i_vec__ = 0; i_vec__ < agg_wts_i_vec_lim__; ++i_vec__) {
            agg_wts[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("death_frac", "3", 3);
        context__.validate_dims("data initialization", "death_frac", "vector_d", context__.to_vec(3));
        validate_non_negative_index("death_frac", "3", 3);
        death_frac = vector_d(static_cast<Eigen::VectorXd::Index>(3));
        vals_r__ = context__.vals_r("death_frac");
        pos__ = 0;
        size_t death_frac_i_vec_lim__ = 3;
        for (size_t i_vec__ = 0; i_vec__ < death_frac_i_vec_lim__; ++i_vec__) {
            death_frac[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("pihat1", "R", R);
        context__.validate_dims("data initialization", "pihat1", "vector_d", context__.to_vec(R));
        validate_non_negative_index("pihat1", "R", R);
        pihat1 = vector_d(static_cast<Eigen::VectorXd::Index>(R));
        vals_r__ = context__.vals_r("pihat1");
        pos__ = 0;
        size_t pihat1_i_vec_lim__ = R;
        for (size_t i_vec__ = 0; i_vec__ < pihat1_i_vec_lim__; ++i_vec__) {
            pihat1[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("pihat_all", "R", R);
        context__.validate_dims("data initialization", "pihat_all", "vector_d", context__.to_vec(R));
        validate_non_negative_index("pihat_all", "R", R);
        pihat_all = vector_d(static_cast<Eigen::VectorXd::Index>(R));
        vals_r__ = context__.vals_r("pihat_all");
        pos__ = 0;
        size_t pihat_all_i_vec_lim__ = R;
        for (size_t i_vec__ = 0; i_vec__ < pihat_all_i_vec_lim__; ++i_vec__) {
            pihat_all[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("a3", "S", S);
        context__.validate_dims("data initialization", "a3", "vector_d", context__.to_vec(S));
        validate_non_negative_index("a3", "S", S);
        a3 = vector_d(static_cast<Eigen::VectorXd::Index>(S));
        vals_r__ = context__.vals_r("a3");
        pos__ = 0;
        size_t a3_i_vec_lim__ = S;
        for (size_t i_vec__ = 0; i_vec__ < a3_i_vec_lim__; ++i_vec__) {
            a3[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("b3", "S", S);
        context__.validate_dims("data initialization", "b3", "vector_d", context__.to_vec(S));
        validate_non_negative_index("b3", "S", S);
        b3 = vector_d(static_cast<Eigen::VectorXd::Index>(S));
        vals_r__ = context__.vals_r("b3");
        pos__ = 0;
        size_t b3_i_vec_lim__ = S;
        for (size_t i_vec__ = 0; i_vec__ < b3_i_vec_lim__; ++i_vec__) {
            b3[i_vec__] = vals_r__[pos__++];
        }
        validate_non_negative_index("N", "A", A);
        validate_non_negative_index("N", "R", R);
        context__.validate_dims("data initialization", "N", "matrix_d", context__.to_vec(A,R));
        validate_non_negative_index("N", "A", A);
        validate_non_negative_index("N", "R", R);
        N = matrix_d(static_cast<Eigen::VectorXd::Index>(A),static_cast<Eigen::VectorXd::Index>(R));
        vals_r__ = context__.vals_r("N");
        pos__ = 0;
        size_t N_m_mat_lim__ = A;
        size_t N_n_mat_lim__ = R;
        for (size_t n_mat__ = 0; n_mat__ < N_n_mat_lim__; ++n_mat__) {
            for (size_t m_mat__ = 0; m_mat__ < N_m_mat_lim__; ++m_mat__) {
                N(m_mat__,n_mat__) = vals_r__[pos__++];
            }
        }
        validate_non_negative_index("D", "A", A);
        validate_non_negative_index("D", "R", R);
        context__.validate_dims("data initialization", "D", "int", context__.to_vec(A,R));
        validate_non_negative_index("D", "A", A);
        validate_non_negative_index("D", "R", R);
        D = std::vector<std::vector<int> >(A,std::vector<int>(R,int(0)));
        vals_i__ = context__.vals_i("D");
        pos__ = 0;
        size_t D_limit_1__ = R;
        for (size_t i_1__ = 0; i_1__ < D_limit_1__; ++i_1__) {
            size_t D_limit_0__ = A;
            for (size_t i_0__ = 0; i_0__ < D_limit_0__; ++i_0__) {
                D[i_0__][i_1__] = vals_i__[pos__++];
            }
        }

        // validate, data variables
        check_greater_or_equal(function__,"R",R,0);
        check_greater_or_equal(function__,"S",S,0);
        check_greater_or_equal(function__,"A",A,0);
        check_greater_or_equal(function__,"K",K,0);
        check_greater_or_equal(function__,"pihat1",pihat1,0);
        check_less_or_equal(function__,"pihat1",pihat1,1);
        check_greater_or_equal(function__,"pihat_all",pihat_all,0);
        check_less_or_equal(function__,"pihat_all",pihat_all,1);
        check_greater_or_equal(function__,"N",N,0);
        for (int k0__ = 0; k0__ < A; ++k0__) {
            for (int k1__ = 0; k1__ < R; ++k1__) {
                check_greater_or_equal(function__,"D[k0__][k1__]",D[k0__][k1__],0);
            }
        }
        // initialize data variables
        validate_non_negative_index("lambda_star", "A", A);
        validate_non_negative_index("lambda_star", "R", R);
        lambda_star = matrix_d(static_cast<Eigen::VectorXd::Index>(A),static_cast<Eigen::VectorXd::Index>(R));
        stan::math::fill(lambda_star,DUMMY_VAR__);

        try {
            current_statement_begin__ = 28;
            stan::math::assign(lambda_star, rep_matrix(std_schedule,R));
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(e,current_statement_begin__);
            // Next line prevents compiler griping about no return
            throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
        }

        // validate transformed data

        // validate, set parameter ranges
        num_params_r__ = 0U;
        param_ranges_i__.clear();
        validate_non_negative_index("alpha", "K", K);
        validate_non_negative_index("alpha", "R", R);
        num_params_r__ += K * R;
        validate_non_negative_index("eta", "4", 4);
        validate_non_negative_index("eta", "R", R);
        num_params_r__ += (4 - 1) * R;
        ++num_params_r__;
        ++num_params_r__;
    }

    ~model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c() { }


    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("alpha")))
            throw std::runtime_error("variable alpha missing");
        vals_r__ = context__.vals_r("alpha");
        pos__ = 0U;
        validate_non_negative_index("alpha", "K", K);
        validate_non_negative_index("alpha", "R", R);
        context__.validate_dims("initialization", "alpha", "matrix_d", context__.to_vec(K,R));
        // generate_declaration alpha
        matrix_d alpha(static_cast<Eigen::VectorXd::Index>(K),static_cast<Eigen::VectorXd::Index>(R));
        for (int j2__ = 0U; j2__ < R; ++j2__)
            for (int j1__ = 0U; j1__ < K; ++j1__)
                alpha(j1__,j2__) = vals_r__[pos__++];
        try {
            writer__.matrix_unconstrain(alpha);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable alpha: ") + e.what());
        }

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

        if (!(context__.contains_r("K1")))
            throw std::runtime_error("variable K1 missing");
        vals_r__ = context__.vals_r("K1");
        pos__ = 0U;
        context__.validate_dims("initialization", "K1", "double", context__.to_vec());
        // generate_declaration K1
        double K1(0);
        K1 = vals_r__[pos__++];
        try {
            writer__.scalar_lb_unconstrain(0,K1);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable K1: ") + e.what());
        }

        if (!(context__.contains_r("Kall")))
            throw std::runtime_error("variable Kall missing");
        vals_r__ = context__.vals_r("Kall");
        pos__ = 0U;
        context__.validate_dims("initialization", "Kall", "double", context__.to_vec());
        // generate_declaration Kall
        double Kall(0);
        Kall = vals_r__[pos__++];
        try {
            writer__.scalar_lb_unconstrain(0,Kall);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable Kall: ") + 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__;

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

        Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic>  alpha;
        (void) alpha;  // dummy to suppress unused var warning
        if (jacobian__)
            alpha = in__.matrix_constrain(K,R,lp__);
        else
            alpha = in__.matrix_constrain(K,R);

        vector<Eigen::Matrix<T__,Eigen::Dynamic,1> > eta;
        size_t dim_eta_0__ = R;
        eta.reserve(dim_eta_0__);
        for (size_t k_0__ = 0; k_0__ < dim_eta_0__; ++k_0__) {
            if (jacobian__)
                eta.push_back(in__.simplex_constrain(4,lp__));
            else
                eta.push_back(in__.simplex_constrain(4));
        }

        T__ K1;
        (void) K1;  // dummy to suppress unused var warning
        if (jacobian__)
            K1 = in__.scalar_lb_constrain(0,lp__);
        else
            K1 = in__.scalar_lb_constrain(0);

        T__ Kall;
        (void) Kall;  // dummy to suppress unused var warning
        if (jacobian__)
            Kall = in__.scalar_lb_constrain(0,lp__);
        else
            Kall = in__.scalar_lb_constrain(0);


        // transformed parameters
        validate_non_negative_index("pi", "R", R);
        validate_non_negative_index("pi", "3", 3);
        Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic>  pi(static_cast<Eigen::VectorXd::Index>(R),static_cast<Eigen::VectorXd::Index>(3));
        (void) pi;  // dummy to suppress unused var warning

        stan::math::initialize(pi, DUMMY_VAR__);
        stan::math::fill(pi,DUMMY_VAR__);
        validate_non_negative_index("piall", "R", R);
        Eigen::Matrix<T__,Eigen::Dynamic,1>  piall(static_cast<Eigen::VectorXd::Index>(R));
        (void) piall;  // dummy to suppress unused var warning

        stan::math::initialize(piall, DUMMY_VAR__);
        stan::math::fill(piall,DUMMY_VAR__);
        T__ agg_pi3;
        (void) agg_pi3;  // dummy to suppress unused var warning

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


        try {
            current_statement_begin__ = 45;
            for (int region = 1; region <= R; ++region) {

                current_statement_begin__ = 47;
                stan::math::assign(get_base1_lhs(pi,region,1,"pi",1), get_base1(get_base1(eta,region,"eta",1),1,"eta",2));
                current_statement_begin__ = 48;
                stan::math::assign(get_base1_lhs(pi,region,2,"pi",1), ((get_base1(get_base1(eta,region,"eta",1),1,"eta",2) + get_base1(get_base1(eta,region,"eta",1),2,"eta",2)) + get_base1(get_base1(eta,region,"eta",1),3,"eta",2)));
                current_statement_begin__ = 49;
                stan::math::assign(get_base1_lhs(pi,region,3,"pi",1), (get_base1(get_base1(eta,region,"eta",1),1,"eta",2) + get_base1(get_base1(eta,region,"eta",1),2,"eta",2)));
            }
            current_statement_begin__ = 53;
            stan::math::assign(piall, multiply(pi,death_frac));
            current_statement_begin__ = 55;
            stan::math::assign(agg_pi3, dot_product(stan::model::rvalue(pi, stan::model::cons_list(stan::model::index_omni(), stan::model::cons_list(stan::model::index_uni(3), stan::model::nil_index_list())), "pi"),agg_wts));
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(e,current_statement_begin__);
            // Next line prevents compiler griping about no return
            throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
        }

        // validate transformed parameters
        for (int i0__ = 0; i0__ < R; ++i0__) {
            for (int i1__ = 0; i1__ < 3; ++i1__) {
                if (stan::math::is_uninitialized(pi(i0__,i1__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: pi" << '[' << i0__ << ']' << '[' << i1__ << ']';
                    throw std::runtime_error(msg__.str());
                }
            }
        }
        for (int i0__ = 0; i0__ < R; ++i0__) {
            if (stan::math::is_uninitialized(piall(i0__))) {
                std::stringstream msg__;
                msg__ << "Undefined transformed parameter: piall" << '[' << i0__ << ']';
                throw std::runtime_error(msg__.str());
            }
        }
        if (stan::math::is_uninitialized(agg_pi3)) {
            std::stringstream msg__;
            msg__ << "Undefined transformed parameter: agg_pi3";
            throw std::runtime_error(msg__.str());
        }

        const char* function__ = "validate transformed params";
        (void) function__;  // dummy to suppress unused var warning
        check_greater_or_equal(function__,"pi",pi,0);
        check_less_or_equal(function__,"pi",pi,1);
        check_greater_or_equal(function__,"piall",piall,0);
        check_less_or_equal(function__,"piall",piall,1);
        check_greater_or_equal(function__,"agg_pi3",agg_pi3,0);
        check_less_or_equal(function__,"agg_pi3",agg_pi3,1);

        // model body
        try {
            {
                validate_non_negative_index("lambda", "A", A);
                validate_non_negative_index("lambda", "R", R);
                Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic>  lambda(static_cast<Eigen::VectorXd::Index>(A),static_cast<Eigen::VectorXd::Index>(R));
                (void) lambda;  // dummy to suppress unused var warning

                stan::math::initialize(lambda, DUMMY_VAR__);
                stan::math::fill(lambda,DUMMY_VAR__);
                validate_non_negative_index("alpha_diff", "(K - 1)", (K - 1));
                validate_non_negative_index("alpha_diff", "R", R);
                Eigen::Matrix<T__,Eigen::Dynamic,Eigen::Dynamic>  alpha_diff(static_cast<Eigen::VectorXd::Index>((K - 1)),static_cast<Eigen::VectorXd::Index>(R));
                (void) alpha_diff;  // dummy to suppress unused var warning

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


                current_statement_begin__ = 63;
                stan::math::assign(lambda, add(lambda_star,multiply(B,alpha)));
                current_statement_begin__ = 65;
                stan::math::assign(alpha_diff, subtract(stan::model::rvalue(alpha, stan::model::cons_list(stan::model::index_min_max(2, K), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "alpha"),stan::model::rvalue(alpha, stan::model::cons_list(stan::model::index_min_max(1, (K - 1)), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "alpha")));
                current_statement_begin__ = 70;
                for (int region = 1; region <= R; ++region) {

                    current_statement_begin__ = 71;
                    lp_accum__.add(poisson_log_log<propto__>(get_base1(get_base1(D,1,"D",1),region,"D",2), ((log(get_base1(N,1,region,"N",1)) + get_base1(lambda,1,region,"lambda",1)) + log(get_base1(pi,region,1,"pi",1)))));
                    current_statement_begin__ = 72;
                    lp_accum__.add(poisson_log_log<propto__>(stan::model::rvalue(D, stan::model::cons_list(stan::model::index_min_max(2, 30), stan::model::cons_list(stan::model::index_uni(region), stan::model::nil_index_list())), "D"), add(add(log(stan::model::rvalue(N, stan::model::cons_list(stan::model::index_min_max(2, 30), stan::model::cons_list(stan::model::index_uni(region), stan::model::nil_index_list())), "N")),stan::model::rvalue(lambda, stan::model::cons_list(stan::model::index_min_max(2, 30), stan::model::cons_list(stan::model::index_uni(region), stan::model::nil_index_list())), "lambda")),log(get_base1(pi,region,2,"pi",1)))));
                    current_statement_begin__ = 73;
                    lp_accum__.add(poisson_log_log<propto__>(stan::model::rvalue(D, stan::model::cons_list(stan::model::index_min_max(31, A), stan::model::cons_list(stan::model::index_uni(region), stan::model::nil_index_list())), "D"), add(add(log(stan::model::rvalue(N, stan::model::cons_list(stan::model::index_min_max(31, A), stan::model::cons_list(stan::model::index_uni(region), stan::model::nil_index_list())), "N")),stan::model::rvalue(lambda, stan::model::cons_list(stan::model::index_min_max(31, A), stan::model::cons_list(stan::model::index_uni(region), stan::model::nil_index_list())), "lambda")),log(get_base1(pi,region,3,"pi",1)))));
                }
                current_statement_begin__ = 77;
                lp_accum__.add(normal_log<propto__>(to_vector(alpha), 0, 4));
                current_statement_begin__ = 78;
                lp_accum__.add(normal_log<propto__>(to_vector(alpha_diff), 0, (1 / sqrt(2))));
                current_statement_begin__ = 80;
                lp_accum__.add(exponential_log<propto__>(K1, 0.050000000000000003));
                current_statement_begin__ = 81;
                lp_accum__.add(exponential_log<propto__>(Kall, 0.050000000000000003));
                current_statement_begin__ = 83;
                lp_accum__.add(beta_log<propto__>(stan::model::rvalue(pi, stan::model::cons_list(stan::model::index_omni(), stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list())), "pi"), multiply((5 + K1),pihat1), multiply((5 + K1),subtract(1,pihat1))));
                current_statement_begin__ = 85;
                lp_accum__.add(beta_log<propto__>(piall, multiply((5 + Kall),pihat_all), multiply((5 + Kall),subtract(1,pihat_all))));
                current_statement_begin__ = 87;
                if (as_bool(logical_gt(S,1))) {
                    current_statement_begin__ = 87;
                    lp_accum__.add(beta_log<propto__>(stan::model::rvalue(pi, stan::model::cons_list(stan::model::index_omni(), stan::model::cons_list(stan::model::index_uni(3), stan::model::nil_index_list())), "pi"), a3, b3));
                }
                current_statement_begin__ = 88;
                if (as_bool(logical_eq(S,1))) {
                    current_statement_begin__ = 88;
                    lp_accum__.add(beta_log<propto__>(agg_pi3, a3, b3));
                }
            }
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(e,current_statement_begin__);
            // 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("alpha");
        names__.push_back("eta");
        names__.push_back("K1");
        names__.push_back("Kall");
        names__.push_back("pi");
        names__.push_back("piall");
        names__.push_back("agg_pi3");
    }


    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(K);
        dims__.push_back(R);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(R);
        dims__.push_back(4);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(R);
        dims__.push_back(3);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(R);
        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__ = "model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c_namespace::write_array";
        (void) function__;  // dummy to suppress unused var warning
        // read-transform, write parameters
        matrix_d alpha = in__.matrix_constrain(K,R);
        vector<vector_d> eta;
        size_t dim_eta_0__ = R;
        for (size_t k_0__ = 0; k_0__ < dim_eta_0__; ++k_0__) {
            eta.push_back(in__.simplex_constrain(4));
        }
        double K1 = in__.scalar_lb_constrain(0);
        double Kall = in__.scalar_lb_constrain(0);
        for (int k_1__ = 0; k_1__ < R; ++k_1__) {
            for (int k_0__ = 0; k_0__ < K; ++k_0__) {
                vars__.push_back(alpha(k_0__, k_1__));
            }
        }
        for (int k_1__ = 0; k_1__ < 4; ++k_1__) {
            for (int k_0__ = 0; k_0__ < R; ++k_0__) {
                vars__.push_back(eta[k_0__][k_1__]);
            }
        }
        vars__.push_back(K1);
        vars__.push_back(Kall);

        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

        validate_non_negative_index("pi", "R", R);
        validate_non_negative_index("pi", "3", 3);
        matrix_d pi(static_cast<Eigen::VectorXd::Index>(R),static_cast<Eigen::VectorXd::Index>(3));
        (void) pi;  // dummy to suppress unused var warning

        stan::math::initialize(pi, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(pi,DUMMY_VAR__);
        validate_non_negative_index("piall", "R", R);
        vector_d piall(static_cast<Eigen::VectorXd::Index>(R));
        (void) piall;  // dummy to suppress unused var warning

        stan::math::initialize(piall, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(piall,DUMMY_VAR__);
        double agg_pi3(0.0);
        (void) agg_pi3;  // dummy to suppress unused var warning

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


        try {
            current_statement_begin__ = 45;
            for (int region = 1; region <= R; ++region) {

                current_statement_begin__ = 47;
                stan::math::assign(get_base1_lhs(pi,region,1,"pi",1), get_base1(get_base1(eta,region,"eta",1),1,"eta",2));
                current_statement_begin__ = 48;
                stan::math::assign(get_base1_lhs(pi,region,2,"pi",1), ((get_base1(get_base1(eta,region,"eta",1),1,"eta",2) + get_base1(get_base1(eta,region,"eta",1),2,"eta",2)) + get_base1(get_base1(eta,region,"eta",1),3,"eta",2)));
                current_statement_begin__ = 49;
                stan::math::assign(get_base1_lhs(pi,region,3,"pi",1), (get_base1(get_base1(eta,region,"eta",1),1,"eta",2) + get_base1(get_base1(eta,region,"eta",1),2,"eta",2)));
            }
            current_statement_begin__ = 53;
            stan::math::assign(piall, multiply(pi,death_frac));
            current_statement_begin__ = 55;
            stan::math::assign(agg_pi3, dot_product(stan::model::rvalue(pi, stan::model::cons_list(stan::model::index_omni(), stan::model::cons_list(stan::model::index_uni(3), stan::model::nil_index_list())), "pi"),agg_wts));
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(e,current_statement_begin__);
            // Next line prevents compiler griping about no return
            throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
        }

        // validate transformed parameters
        check_greater_or_equal(function__,"pi",pi,0);
        check_less_or_equal(function__,"pi",pi,1);
        check_greater_or_equal(function__,"piall",piall,0);
        check_less_or_equal(function__,"piall",piall,1);
        check_greater_or_equal(function__,"agg_pi3",agg_pi3,0);
        check_less_or_equal(function__,"agg_pi3",agg_pi3,1);

        // write transformed parameters
        for (int k_1__ = 0; k_1__ < 3; ++k_1__) {
            for (int k_0__ = 0; k_0__ < R; ++k_0__) {
                vars__.push_back(pi(k_0__, k_1__));
            }
        }
        for (int k_0__ = 0; k_0__ < R; ++k_0__) {
            vars__.push_back(piall[k_0__]);
        }
        vars__.push_back(agg_pi3);

        if (!include_gqs__) return;
        // declare and define generated quantities


        try {
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(e,current_statement_begin__);
            // Next line prevents compiler griping about no return
            throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
        }

        // validate generated quantities

        // write generated quantities
    }

    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 "model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c";
    }


    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_1__ = 1; k_1__ <= R; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= K; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "alpha" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_1__ = 1; k_1__ <= 4; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= R; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "eta" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        param_name_stream__.str(std::string());
        param_name_stream__ << "K1";
        param_names__.push_back(param_name_stream__.str());
        param_name_stream__.str(std::string());
        param_name_stream__ << "Kall";
        param_names__.push_back(param_name_stream__.str());

        if (!include_gqs__ && !include_tparams__) return;
        for (int k_1__ = 1; k_1__ <= 3; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= R; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "pi" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_0__ = 1; k_0__ <= R; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "piall" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        param_name_stream__.str(std::string());
        param_name_stream__ << "agg_pi3";
        param_names__.push_back(param_name_stream__.str());

        if (!include_gqs__) return;
    }


    void unconstrained_param_names(std::vector<std::string>& param_names__,
                                   bool include_tparams__ = true,
                                   bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        for (int k_1__ = 1; k_1__ <= R; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= K; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "alpha" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_1__ = 1; k_1__ <= (4 - 1); ++k_1__) {
            for (int k_0__ = 1; k_0__ <= R; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "eta" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        param_name_stream__.str(std::string());
        param_name_stream__ << "K1";
        param_names__.push_back(param_name_stream__.str());
        param_name_stream__.str(std::string());
        param_name_stream__ << "Kall";
        param_names__.push_back(param_name_stream__.str());

        if (!include_gqs__ && !include_tparams__) return;
        for (int k_1__ = 1; k_1__ <= 3; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= R; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "pi" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_0__ = 1; k_0__ <= R; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "piall" << '.' << k_0__;
            param_names__.push_back(param_name_stream__.str());
        }
        param_name_stream__.str(std::string());
        param_name_stream__ << "agg_pi3";
        param_names__.push_back(param_name_stream__.str());

        if (!include_gqs__) return;
    }

}; // model

}

typedef model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c_namespace::model6c44f047e07_dd37d1dbcaf77db616116e448a86dd9c stan_model;