// Code generated by Stan version 2.19.1

#include <stan/model/model_header.hpp>

namespace model3083780ecf59_multistage_mark_recapture_namespace {

using std::istream;
using std::string;
using std::stringstream;
using std::vector;
using stan::io::dump;
using stan::math::lgamma;
using stan::model::prob_grad;
using namespace stan::math;

static int current_statement_begin__;

stan::io::program_reader prog_reader__() {
    stan::io::program_reader reader;
    reader.add_event(0, 0, "start", "model3083780ecf59_multistage_mark_recapture");
    reader.add_event(267, 265, "end", "model3083780ecf59_multistage_mark_recapture");
    return reader;
}

template <typename T3__, typename T4__, typename T5__>
typename boost::math::tools::promote_args<T3__, T4__, T5__>::type
impute_missing(const int& k,
                   const int& T0,
                   const int& T1,
                   const std::vector<std::vector<T3__> >& p,
                   const std::vector<std::vector<T4__> >& s,
                   const std::vector<std::vector<T5__> >& f, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T3__, T4__, T5__>::type local_scalar_t__;
    typedef local_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        local_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__ = 4;
        local_scalar_t__ prob(DUMMY_VAR__);
        (void) prob;  // dummy to suppress unused var warning
        stan::math::initialize(prob, DUMMY_VAR__);
        stan::math::fill(prob, DUMMY_VAR__);

        current_statement_begin__ = 5;
        int T(0);
        (void) T;  // dummy to suppress unused var warning
        stan::math::fill(T, std::numeric_limits<int>::min());
        stan::math::assign(T,(T1 - T0));

        current_statement_begin__ = 6;
        validate_non_negative_index("A", "T", T);
        validate_non_negative_index("A", "T", T);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> A(T, T);
        stan::math::initialize(A, DUMMY_VAR__);
        stan::math::fill(A, DUMMY_VAR__);


        current_statement_begin__ = 7;
        for (int row = 1; row <= T; ++row) {

            current_statement_begin__ = 8;
            if (as_bool(logical_gt(row, 1))) {

                current_statement_begin__ = 9;
                stan::model::assign(A, 
                            stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list())), 
                            ((get_base1(get_base1(s, k, "s", 1), T0, "s", 2) * (1 - get_base1(get_base1(f, k, "f", 1), T0, "f", 2))) * (1 - get_base1(get_base1(p, 1, "p", 1), (T0 + 1), "p", 2))), 
                            "assigning variable A");
            } else {

                current_statement_begin__ = 13;
                stan::model::assign(A, 
                            stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list())), 
                            ((1 - get_base1(get_base1(s, k, "s", 1), T0, "s", 2)) * (1 - get_base1(get_base1(p, 2, "p", 1), (T0 + 1), "p", 2))), 
                            "assigning variable A");
            }
            current_statement_begin__ = 16;
            for (int col = 2; col <= T; ++col) {

                current_statement_begin__ = 17;
                if (as_bool(logical_lt(col, T))) {

                    current_statement_begin__ = 18;
                    if (as_bool(logical_gt(col, row))) {

                        current_statement_begin__ = 19;
                        stan::model::assign(A, 
                                    stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(col), stan::model::nil_index_list())), 
                                    (get_base1(A, row, (col - 1), "A", 1) * (1 - get_base1(get_base1(p, 2, "p", 1), (T0 + col), "p", 2))), 
                                    "assigning variable A");
                    } else if (as_bool(logical_eq(col, row))) {

                        current_statement_begin__ = 22;
                        stan::model::assign(A, 
                                    stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(col), stan::model::nil_index_list())), 
                                    ((get_base1(A, row, (col - 1), "A", 1) * (1 - get_base1(get_base1(s, 1, "s", 1), ((T0 - 1) + col), "s", 2))) * (1 - get_base1(get_base1(p, 2, "p", 1), (T0 + col), "p", 2))), 
                                    "assigning variable A");
                    } else {

                        current_statement_begin__ = 26;
                        stan::model::assign(A, 
                                    stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(col), stan::model::nil_index_list())), 
                                    (((get_base1(A, row, (col - 1), "A", 1) * get_base1(get_base1(s, 1, "s", 1), ((T0 - 1) + col), "s", 2)) * (1 - get_base1(get_base1(f, 1, "f", 1), ((T0 - 1) + col), "f", 2))) * (1 - get_base1(get_base1(p, 1, "p", 1), (T0 + col), "p", 2))), 
                                    "assigning variable A");
                    }
                } else {

                    current_statement_begin__ = 32;
                    if (as_bool(logical_lt(row, T))) {

                        current_statement_begin__ = 33;
                        stan::model::assign(A, 
                                    stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(T), stan::model::nil_index_list())), 
                                    (get_base1(A, row, (T - 1), "A", 1) * get_base1(get_base1(p, 2, "p", 1), T1, "p", 2)), 
                                    "assigning variable A");
                    } else {

                        current_statement_begin__ = 36;
                        stan::model::assign(A, 
                                    stan::model::cons_list(stan::model::index_uni(row), stan::model::cons_list(stan::model::index_uni(T), stan::model::nil_index_list())), 
                                    ((get_base1(A, row, (T - 1), "A", 1) * (1 - get_base1(get_base1(s, 1, "s", 1), (T1 - 1), "s", 2))) * get_base1(get_base1(p, 2, "p", 1), T1, "p", 2)), 
                                    "assigning variable A");
                    }
                }
            }
        }
        current_statement_begin__ = 44;
        stan::math::assign(prob, 0.0);
        current_statement_begin__ = 45;
        for (int i = 1; i <= T; ++i) {
            current_statement_begin__ = 46;
            stan::math::assign(prob, (prob + get_base1(A, i, T, "A", 1)));
        }
        current_statement_begin__ = 47;
        return stan::math::promote_scalar<fun_return_scalar_t__>(prob);
        }
    } 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 impute_missing_functor__ {
    template <typename T3__, typename T4__, typename T5__>
        typename boost::math::tools::promote_args<T3__, T4__, T5__>::type
    operator()(const int& k,
                   const int& T0,
                   const int& T1,
                   const std::vector<std::vector<T3__> >& p,
                   const std::vector<std::vector<T4__> >& s,
                   const std::vector<std::vector<T5__> >& f, std::ostream* pstream__) const {
        return impute_missing(k, T0, T1, p, s, f, pstream__);
    }
};

int
first_obs(const std::vector<int>& y_i, std::ostream* pstream__) {
    typedef double local_scalar_t__;
    typedef int fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        local_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__ = 51;
        for (int k = 1; k <= size(y_i); ++k) {

            current_statement_begin__ = 52;
            if (as_bool(get_base1(y_i, k, "y_i", 1))) {
                current_statement_begin__ = 53;
                return stan::math::promote_scalar<fun_return_scalar_t__>(k);
            }
        }
        current_statement_begin__ = 55;
        return stan::math::promote_scalar<fun_return_scalar_t__>(0);
    } 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 first_obs_functor__ {
            int
    operator()(const std::vector<int>& y_i, std::ostream* pstream__) const {
        return first_obs(y_i, pstream__);
    }
};

int
last_obs(const std::vector<int>& y_i,
             const int& i, std::ostream* pstream__) {
    typedef double local_scalar_t__;
    typedef int fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        local_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__ = 59;
        for (int k = 0; k <= (size(y_i) - 1); ++k) {

            current_statement_begin__ = 60;
            if (as_bool(get_base1(y_i, (size(y_i) - k), "y_i", 1))) {

                current_statement_begin__ = 61;
                return stan::math::promote_scalar<fun_return_scalar_t__>((size(y_i) - k));
            }
        }
        current_statement_begin__ = 64;
        return stan::math::promote_scalar<fun_return_scalar_t__>(-(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 last_obs_functor__ {
            int
    operator()(const std::vector<int>& y_i,
             const int& i, std::ostream* pstream__) const {
        return last_obs(y_i, i, pstream__);
    }
};

int
find_block_start(const std::vector<int>& y_i,
                     const int& start_ind,
                     const int& i, std::ostream* pstream__) {
    typedef double local_scalar_t__;
    typedef int fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        local_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__ = 68;
        int kback(0);
        (void) kback;  // dummy to suppress unused var warning
        stan::math::fill(kback, std::numeric_limits<int>::min());


        current_statement_begin__ = 69;
        for (int k = 1; k <= start_ind; ++k) {

            current_statement_begin__ = 70;
            stan::math::assign(kback, (start_ind - k));
            current_statement_begin__ = 71;
            if (as_bool(get_base1(y_i, kback, "y_i", 1))) {

                current_statement_begin__ = 72;
                return stan::math::promote_scalar<fun_return_scalar_t__>(kback);
            }
        }
        current_statement_begin__ = 75;
        return stan::math::promote_scalar<fun_return_scalar_t__>(0);
        }
    } 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 find_block_start_functor__ {
            int
    operator()(const std::vector<int>& y_i,
                     const int& start_ind,
                     const int& i, std::ostream* pstream__) const {
        return find_block_start(y_i, start_ind, i, pstream__);
    }
};

class model3083780ecf59_multistage_mark_recapture : public prob_grad {
private:
        int T_years;
        int N_plant;
        std::vector<std::vector<int> > y;
public:
    model3083780ecf59_multistage_mark_recapture(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, 0, pstream__);
    }

    model3083780ecf59_multistage_mark_recapture(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__) {
        typedef double local_scalar_t__;

        boost::ecuyer1988 base_rng__ =
          stan::services::util::create_rng(random_seed__, 0);
        (void) base_rng__;  // suppress unused var warning

        current_statement_begin__ = -1;

        static const char* function__ = "model3083780ecf59_multistage_mark_recapture_namespace::model3083780ecf59_multistage_mark_recapture";
        (void) function__;  // dummy to suppress unused var warning
        size_t pos__;
        (void) pos__;  // dummy to suppress unused var warning
        std::vector<int> vals_i__;
        std::vector<double> vals_r__;
        local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

        try {
            // initialize data block variables from context__
            current_statement_begin__ = 80;
            context__.validate_dims("data initialization", "T_years", "int", context__.to_vec());
            T_years = int(0);
            vals_i__ = context__.vals_i("T_years");
            pos__ = 0;
            T_years = vals_i__[pos__++];
            check_greater_or_equal(function__, "T_years", T_years, 0);

            current_statement_begin__ = 81;
            context__.validate_dims("data initialization", "N_plant", "int", context__.to_vec());
            N_plant = int(0);
            vals_i__ = context__.vals_i("N_plant");
            pos__ = 0;
            N_plant = vals_i__[pos__++];
            check_greater_or_equal(function__, "N_plant", N_plant, 0);

            current_statement_begin__ = 82;
            validate_non_negative_index("y", "N_plant", N_plant);
            validate_non_negative_index("y", "T_years", T_years);
            context__.validate_dims("data initialization", "y", "int", context__.to_vec(N_plant,T_years));
            y = std::vector<std::vector<int> >(N_plant, std::vector<int>(T_years, int(0)));
            vals_i__ = context__.vals_i("y");
            pos__ = 0;
            size_t y_k_0_max__ = N_plant;
            size_t y_k_1_max__ = T_years;
            for (size_t k_1__ = 0; k_1__ < y_k_1_max__; ++k_1__) {
                for (size_t k_0__ = 0; k_0__ < y_k_0_max__; ++k_0__) {
                    y[k_0__][k_1__] = vals_i__[pos__++];
                }
            }
            size_t y_i_0_max__ = N_plant;
            size_t y_i_1_max__ = T_years;
            for (size_t i_0__ = 0; i_0__ < y_i_0_max__; ++i_0__) {
                for (size_t i_1__ = 0; i_1__ < y_i_1_max__; ++i_1__) {
                    check_greater_or_equal(function__, "y[i_0__][i_1__]", y[i_0__][i_1__], 0);
                    check_less_or_equal(function__, "y[i_0__][i_1__]", y[i_0__][i_1__], 3);
                }
            }


            // initialize transformed data variables
            // execute transformed data statements

            // validate transformed data

            // validate, set parameter ranges
            num_params_r__ = 0U;
            param_ranges_i__.clear();
            current_statement_begin__ = 85;
            validate_non_negative_index("p", "2", 2);
            validate_non_negative_index("p", "T_years", T_years);
            num_params_r__ += ((1 * 2) * T_years);
            current_statement_begin__ = 88;
            validate_non_negative_index("s", "2", 2);
            validate_non_negative_index("s", "(T_years - 1)", (T_years - 1));
            num_params_r__ += ((1 * 2) * (T_years - 1));
            current_statement_begin__ = 91;
            validate_non_negative_index("f", "2", 2);
            validate_non_negative_index("f", "(T_years - 1)", (T_years - 1));
            num_params_r__ += ((1 * 2) * (T_years - 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 ***");
        }
    }

    ~model3083780ecf59_multistage_mark_recapture() { }


    void transform_inits(const stan::io::var_context& context__,
                         std::vector<int>& params_i__,
                         std::vector<double>& params_r__,
                         std::ostream* pstream__) const {
        typedef double local_scalar_t__;
        stan::io::writer<double> writer__(params_r__, params_i__);
        size_t pos__;
        (void) pos__; // dummy call to supress warning
        std::vector<double> vals_r__;
        std::vector<int> vals_i__;

        current_statement_begin__ = 85;
        if (!(context__.contains_r("p")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable p missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("p");
        pos__ = 0U;
        validate_non_negative_index("p", "2", 2);
        validate_non_negative_index("p", "T_years", T_years);
        context__.validate_dims("parameter initialization", "p", "double", context__.to_vec(2,T_years));
        std::vector<std::vector<double> > p(2, std::vector<double>(T_years, double(0)));
        size_t p_k_0_max__ = 2;
        size_t p_k_1_max__ = T_years;
        for (size_t k_1__ = 0; k_1__ < p_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < p_k_0_max__; ++k_0__) {
                p[k_0__][k_1__] = vals_r__[pos__++];
            }
        }
        size_t p_i_0_max__ = 2;
        size_t p_i_1_max__ = T_years;
        for (size_t i_0__ = 0; i_0__ < p_i_0_max__; ++i_0__) {
            for (size_t i_1__ = 0; i_1__ < p_i_1_max__; ++i_1__) {
                try {
                    writer__.scalar_lub_unconstrain(0, 1, p[i_0__][i_1__]);
                } catch (const std::exception& e) {
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable p: ") + e.what()), current_statement_begin__, prog_reader__());
                }
            }
        }

        current_statement_begin__ = 88;
        if (!(context__.contains_r("s")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable s missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("s");
        pos__ = 0U;
        validate_non_negative_index("s", "2", 2);
        validate_non_negative_index("s", "(T_years - 1)", (T_years - 1));
        context__.validate_dims("parameter initialization", "s", "double", context__.to_vec(2,(T_years - 1)));
        std::vector<std::vector<double> > s(2, std::vector<double>((T_years - 1), double(0)));
        size_t s_k_0_max__ = 2;
        size_t s_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < s_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < s_k_0_max__; ++k_0__) {
                s[k_0__][k_1__] = vals_r__[pos__++];
            }
        }
        size_t s_i_0_max__ = 2;
        size_t s_i_1_max__ = (T_years - 1);
        for (size_t i_0__ = 0; i_0__ < s_i_0_max__; ++i_0__) {
            for (size_t i_1__ = 0; i_1__ < s_i_1_max__; ++i_1__) {
                try {
                    writer__.scalar_lub_unconstrain(0, 1, s[i_0__][i_1__]);
                } catch (const std::exception& e) {
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable s: ") + e.what()), current_statement_begin__, prog_reader__());
                }
            }
        }

        current_statement_begin__ = 91;
        if (!(context__.contains_r("f")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable f missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("f");
        pos__ = 0U;
        validate_non_negative_index("f", "2", 2);
        validate_non_negative_index("f", "(T_years - 1)", (T_years - 1));
        context__.validate_dims("parameter initialization", "f", "double", context__.to_vec(2,(T_years - 1)));
        std::vector<std::vector<double> > f(2, std::vector<double>((T_years - 1), double(0)));
        size_t f_k_0_max__ = 2;
        size_t f_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < f_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < f_k_0_max__; ++k_0__) {
                f[k_0__][k_1__] = vals_r__[pos__++];
            }
        }
        size_t f_i_0_max__ = 2;
        size_t f_i_1_max__ = (T_years - 1);
        for (size_t i_0__ = 0; i_0__ < f_i_0_max__; ++i_0__) {
            for (size_t i_1__ = 0; i_1__ < f_i_1_max__; ++i_1__) {
                try {
                    writer__.scalar_lub_unconstrain(0, 1, f[i_0__][i_1__]);
                } catch (const std::exception& e) {
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable f: ") + e.what()), current_statement_begin__, prog_reader__());
                }
            }
        }

        params_r__ = writer__.data_r();
        params_i__ = writer__.data_i();
    }

    void transform_inits(const stan::io::var_context& context,
                         Eigen::Matrix<double, Eigen::Dynamic, 1>& params_r,
                         std::ostream* pstream__) const {
      std::vector<double> params_r_vec;
      std::vector<int> params_i_vec;
      transform_inits(context, params_i_vec, params_r_vec, pstream__);
      params_r.resize(params_r_vec.size());
      for (int i = 0; i < params_r.size(); ++i)
        params_r(i) = params_r_vec[i];
    }


    template <bool propto__, bool jacobian__, typename T__>
    T__ log_prob(std::vector<T__>& params_r__,
                 std::vector<int>& params_i__,
                 std::ostream* pstream__ = 0) const {

        typedef T__ local_scalar_t__;

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

        T__ lp__(0.0);
        stan::math::accumulator<T__> lp_accum__;
        try {
            stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);

            // model parameters
            current_statement_begin__ = 85;
            std::vector<std::vector<local_scalar_t__> > p;
            size_t p_d_0_max__ = 2;
            size_t p_d_1_max__ = T_years;
            p.resize(p_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < p_d_0_max__; ++d_0__) {
                p[d_0__].reserve(p_d_1_max__);
                for (size_t d_1__ = 0; d_1__ < p_d_1_max__; ++d_1__) {
                    if (jacobian__)
                        p[d_0__].push_back(in__.scalar_lub_constrain(0, 1, lp__));
                    else
                        p[d_0__].push_back(in__.scalar_lub_constrain(0, 1));
                }
            }

            current_statement_begin__ = 88;
            std::vector<std::vector<local_scalar_t__> > s;
            size_t s_d_0_max__ = 2;
            size_t s_d_1_max__ = (T_years - 1);
            s.resize(s_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < s_d_0_max__; ++d_0__) {
                s[d_0__].reserve(s_d_1_max__);
                for (size_t d_1__ = 0; d_1__ < s_d_1_max__; ++d_1__) {
                    if (jacobian__)
                        s[d_0__].push_back(in__.scalar_lub_constrain(0, 1, lp__));
                    else
                        s[d_0__].push_back(in__.scalar_lub_constrain(0, 1));
                }
            }

            current_statement_begin__ = 91;
            std::vector<std::vector<local_scalar_t__> > f;
            size_t f_d_0_max__ = 2;
            size_t f_d_1_max__ = (T_years - 1);
            f.resize(f_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < f_d_0_max__; ++d_0__) {
                f[d_0__].reserve(f_d_1_max__);
                for (size_t d_1__ = 0; d_1__ < f_d_1_max__; ++d_1__) {
                    if (jacobian__)
                        f[d_0__].push_back(in__.scalar_lub_constrain(0, 1, lp__));
                    else
                        f[d_0__].push_back(in__.scalar_lub_constrain(0, 1));
                }
            }

            // model body

            current_statement_begin__ = 102;
            for (int i = 1; i <= 2; ++i) {
                current_statement_begin__ = 103;
                for (int j = 1; j <= (T_years - 1); ++j) {

                    current_statement_begin__ = 104;
                    lp_accum__.add(beta_log<propto__>(get_base1(get_base1(s, i, "s", 1), j, "s", 2), 12, 2));
                    current_statement_begin__ = 105;
                    lp_accum__.add(beta_log<propto__>(get_base1(get_base1(f, i, "f", 1), j, "f", 2), 2, 4));
                    current_statement_begin__ = 106;
                    lp_accum__.add(beta_log<propto__>(get_base1(get_base1(p, i, "p", 1), j, "p", 2), 0.75, 0.8));
                }
            }
            current_statement_begin__ = 109;
            for (int i = 1; i <= N_plant; ++i) {
                {
                current_statement_begin__ = 110;
                int first_year(0);
                (void) first_year;  // dummy to suppress unused var warning
                stan::math::fill(first_year, std::numeric_limits<int>::min());

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


                current_statement_begin__ = 112;
                stan::math::assign(first_year, first_obs(get_base1(y, i, "y", 1), pstream__));
                current_statement_begin__ = 113;
                stan::math::assign(last_year, last_obs(get_base1(y, i, "y", 1), i, pstream__));
                current_statement_begin__ = 118;
                if (as_bool(logical_lt(get_base1(get_base1(y, i, "y", 1), last_year, "y", 2), 3))) {

                    current_statement_begin__ = 120;
                    for (int t = (first_year + 1); t <= last_year; ++t) {

                        current_statement_begin__ = 123;
                        if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), t, "y", 2), 1))) {

                            current_statement_begin__ = 125;
                            lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(p, 1, "p", 1), t, "p", 2)));
                            current_statement_begin__ = 127;
                            if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (t - 1), "y", 2), 2))) {

                                current_statement_begin__ = 129;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 2, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 131;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (t - 1), "s", 2)));
                            } else {

                                current_statement_begin__ = 135;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 1, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 137;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (t - 1), "s", 2)));
                            }
                        } else if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), t, "y", 2), 2))) {

                            current_statement_begin__ = 143;
                            if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (t - 1), "y", 2), 2))) {

                                current_statement_begin__ = 145;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(f, 2, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 147;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (t - 1), "s", 2)));
                            } else {

                                current_statement_begin__ = 151;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(f, 1, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 153;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (t - 1), "s", 2)));
                            }
                        } else {

                            current_statement_begin__ = 158;
                            lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(p, 1, "p", 1), t, "p", 2)));
                            current_statement_begin__ = 160;
                            if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (t - 1), "y", 2), 2))) {

                                current_statement_begin__ = 162;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 2, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 164;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (t - 1), "s", 2)));
                            } else {

                                current_statement_begin__ = 168;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 1, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 170;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (t - 1), "s", 2)));
                            }
                        }
                    }
                } else {
                    {
                    current_statement_begin__ = 179;
                    int block_start(0);
                    (void) block_start;  // dummy to suppress unused var warning
                    stan::math::fill(block_start, std::numeric_limits<int>::min());
                    stan::math::assign(block_start,find_block_start(get_base1(y, i, "y", 1), last_year, i, pstream__));


                    current_statement_begin__ = 185;
                    if (as_bool(get_base1(get_base1(y, i, "y", 1), (last_year - 1), "y", 2))) {

                        current_statement_begin__ = 187;
                        lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(p, 2, "p", 1), last_year, "p", 2)));
                        current_statement_begin__ = 189;
                        if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (last_year - 1), "y", 2), 2))) {

                            current_statement_begin__ = 191;
                            lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 2, "f", 1), (last_year - 1), "f", 2)));
                            current_statement_begin__ = 193;
                            lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (last_year - 1), "s", 2)));
                        } else {

                            current_statement_begin__ = 197;
                            lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 1, "f", 1), (last_year - 1), "f", 2)));
                            current_statement_begin__ = 199;
                            lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (last_year - 1), "s", 2)));
                        }
                    } else {
                        {
                        current_statement_begin__ = 202;
                        int k_ind(0);
                        (void) k_ind;  // dummy to suppress unused var warning
                        stan::math::fill(k_ind, std::numeric_limits<int>::min());
                        stan::math::assign(k_ind,get_base1(get_base1(y, i, "y", 1), block_start, "y", 2));

                        current_statement_begin__ = 203;
                        local_scalar_t__ prob_combined(DUMMY_VAR__);
                        (void) prob_combined;  // dummy to suppress unused var warning
                        stan::math::initialize(prob_combined, DUMMY_VAR__);
                        stan::math::fill(prob_combined, DUMMY_VAR__);
                        stan::math::assign(prob_combined,impute_missing(k_ind, block_start, last_year, p, s, f, pstream__));


                        current_statement_begin__ = 204;
                        if (pstream__) {
                            stan_print(pstream__,prob_combined);
                            *pstream__ << std::endl;
                        }
                        current_statement_begin__ = 206;
                        lp_accum__.add(bernoulli_log<propto__>(1, prob_combined));
                        }
                    }
                    current_statement_begin__ = 209;
                    for (int t = (first_year + 1); t <= block_start; ++t) {

                        current_statement_begin__ = 212;
                        if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), t, "y", 2), 1))) {

                            current_statement_begin__ = 214;
                            lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(p, 1, "p", 1), t, "p", 2)));
                            current_statement_begin__ = 216;
                            if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (t - 1), "y", 2), 2))) {

                                current_statement_begin__ = 218;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 2, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 220;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (t - 1), "s", 2)));
                            } else {

                                current_statement_begin__ = 224;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 1, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 226;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (t - 1), "s", 2)));
                            }
                        } else if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), t, "y", 2), 2))) {

                            current_statement_begin__ = 232;
                            if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (t - 1), "y", 2), 2))) {

                                current_statement_begin__ = 234;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(f, 2, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 236;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (t - 1), "s", 2)));
                            } else {

                                current_statement_begin__ = 240;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(f, 1, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 242;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (t - 1), "s", 2)));
                            }
                        } else {

                            current_statement_begin__ = 247;
                            lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(p, 1, "p", 1), t, "p", 2)));
                            current_statement_begin__ = 249;
                            if (as_bool(logical_eq(get_base1(get_base1(y, i, "y", 1), (t - 1), "y", 2), 2))) {

                                current_statement_begin__ = 251;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 2, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 253;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 2, "s", 1), (t - 1), "s", 2)));
                            } else {

                                current_statement_begin__ = 257;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f, 1, "f", 1), (t - 1), "f", 2)));
                                current_statement_begin__ = 259;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s, 1, "s", 1), (t - 1), "s", 2)));
                            }
                        }
                    }
                    }
                }
                }
            }

        } 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("p");
        names__.push_back("s");
        names__.push_back("f");
    }


    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(2);
        dims__.push_back(T_years);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(2);
        dims__.push_back((T_years - 1));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(2);
        dims__.push_back((T_years - 1));
        dimss__.push_back(dims__);
    }

    template <typename RNG>
    void write_array(RNG& base_rng__,
                     std::vector<double>& params_r__,
                     std::vector<int>& params_i__,
                     std::vector<double>& vars__,
                     bool include_tparams__ = true,
                     bool include_gqs__ = true,
                     std::ostream* pstream__ = 0) const {
        typedef double local_scalar_t__;

        vars__.resize(0);
        stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
        static const char* function__ = "model3083780ecf59_multistage_mark_recapture_namespace::write_array";
        (void) function__;  // dummy to suppress unused var warning

        // read-transform, write parameters
        std::vector<std::vector<double> > p;
        size_t p_d_0_max__ = 2;
        size_t p_d_1_max__ = T_years;
        p.resize(p_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < p_d_0_max__; ++d_0__) {
            p[d_0__].reserve(p_d_1_max__);
            for (size_t d_1__ = 0; d_1__ < p_d_1_max__; ++d_1__) {
                p[d_0__].push_back(in__.scalar_lub_constrain(0, 1));
            }
        }
        size_t p_k_0_max__ = 2;
        size_t p_k_1_max__ = T_years;
        for (size_t k_1__ = 0; k_1__ < p_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < p_k_0_max__; ++k_0__) {
                vars__.push_back(p[k_0__][k_1__]);
            }
        }

        std::vector<std::vector<double> > s;
        size_t s_d_0_max__ = 2;
        size_t s_d_1_max__ = (T_years - 1);
        s.resize(s_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < s_d_0_max__; ++d_0__) {
            s[d_0__].reserve(s_d_1_max__);
            for (size_t d_1__ = 0; d_1__ < s_d_1_max__; ++d_1__) {
                s[d_0__].push_back(in__.scalar_lub_constrain(0, 1));
            }
        }
        size_t s_k_0_max__ = 2;
        size_t s_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < s_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < s_k_0_max__; ++k_0__) {
                vars__.push_back(s[k_0__][k_1__]);
            }
        }

        std::vector<std::vector<double> > f;
        size_t f_d_0_max__ = 2;
        size_t f_d_1_max__ = (T_years - 1);
        f.resize(f_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < f_d_0_max__; ++d_0__) {
            f[d_0__].reserve(f_d_1_max__);
            for (size_t d_1__ = 0; d_1__ < f_d_1_max__; ++d_1__) {
                f[d_0__].push_back(in__.scalar_lub_constrain(0, 1));
            }
        }
        size_t f_k_0_max__ = 2;
        size_t f_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < f_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < f_k_0_max__; ++k_0__) {
                vars__.push_back(f[k_0__][k_1__]);
            }
        }

        double lp__ = 0.0;
        (void) lp__;  // dummy to suppress unused var warning
        stan::math::accumulator<double> lp_accum__;

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

        if (!include_tparams__ && !include_gqs__) return;

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

    template <typename RNG>
    void write_array(RNG& base_rng,
                     Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,
                     Eigen::Matrix<double,Eigen::Dynamic,1>& vars,
                     bool include_tparams = true,
                     bool include_gqs = true,
                     std::ostream* pstream = 0) const {
      std::vector<double> params_r_vec(params_r.size());
      for (int i = 0; i < params_r.size(); ++i)
        params_r_vec[i] = params_r(i);
      std::vector<double> vars_vec;
      std::vector<int> params_i_vec;
      write_array(base_rng, params_r_vec, params_i_vec, vars_vec, include_tparams, include_gqs, pstream);
      vars.resize(vars_vec.size());
      for (int i = 0; i < vars.size(); ++i)
        vars(i) = vars_vec[i];
    }

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


    void constrained_param_names(std::vector<std::string>& param_names__,
                                 bool include_tparams__ = true,
                                 bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        size_t p_k_0_max__ = 2;
        size_t p_k_1_max__ = T_years;
        for (size_t k_1__ = 0; k_1__ < p_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < p_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "p" << '.' << k_0__ + 1 << '.' << k_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t s_k_0_max__ = 2;
        size_t s_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < s_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < s_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "s" << '.' << k_0__ + 1 << '.' << k_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t f_k_0_max__ = 2;
        size_t f_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < f_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < f_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "f" << '.' << k_0__ + 1 << '.' << k_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
        }

        if (!include_gqs__) return;
    }


    void unconstrained_param_names(std::vector<std::string>& param_names__,
                                   bool include_tparams__ = true,
                                   bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        size_t p_k_0_max__ = 2;
        size_t p_k_1_max__ = T_years;
        for (size_t k_1__ = 0; k_1__ < p_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < p_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "p" << '.' << k_0__ + 1 << '.' << k_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t s_k_0_max__ = 2;
        size_t s_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < s_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < s_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "s" << '.' << k_0__ + 1 << '.' << k_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t f_k_0_max__ = 2;
        size_t f_k_1_max__ = (T_years - 1);
        for (size_t k_1__ = 0; k_1__ < f_k_1_max__; ++k_1__) {
            for (size_t k_0__ = 0; k_0__ < f_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "f" << '.' << k_0__ + 1 << '.' << k_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
        }

        if (!include_gqs__) return;
    }

}; // model

}  // namespace

typedef model3083780ecf59_multistage_mark_recapture_namespace::model3083780ecf59_multistage_mark_recapture stan_model;