// Code generated by Stan version 2.17.0

#include <stan/model/model_header.hpp>

namespace modeld3143b1c2538_mark_simple_ends_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", "modeld3143b1c2538_mark_simple_ends");
    reader.add_event(292, 292, "end", "modeld3143b1c2538_mark_simple_ends");
    return reader;
}

template <typename T3__, typename T4__, typename T5__>
typename boost::math::tools::promote_args<T3__, T4__, T5__>::type
impute_observed(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 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__ = 4;
        fun_scalar_t__ prob;
        (void) prob;  // dummy to suppress unused var warning

        stan::math::initialize(prob, std::numeric_limits<double>::quiet_NaN());
        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<fun_scalar_t__,Eigen::Dynamic,Eigen::Dynamic>  A(static_cast<Eigen::VectorXd::Index>(T),static_cast<Eigen::VectorXd::Index>(T));
        (void) A;  // dummy to suppress unused var warning

        stan::math::initialize(A, std::numeric_limits<double>::quiet_NaN());
        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::math::assign(get_base1_lhs(A,row,1,"A",1), ((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))));
            } else {

                current_statement_begin__ = 13;
                stan::math::assign(get_base1_lhs(A,row,1,"A",1), ((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))));
            }
            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::math::assign(get_base1_lhs(A,row,col,"A",1), (get_base1(A,row,(col - 1),"A",1) * (1 - get_base1(get_base1(p,2,"p",1),(T0 + col),"p",2))));
                    } else if (as_bool(logical_eq(col,row))) {

                        current_statement_begin__ = 22;
                        stan::math::assign(get_base1_lhs(A,row,col,"A",1), ((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))));
                    } else {

                        current_statement_begin__ = 26;
                        stan::math::assign(get_base1_lhs(A,row,col,"A",1), (((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))));
                    }
                } else {

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

                        current_statement_begin__ = 33;
                        stan::math::assign(get_base1_lhs(A,row,T,"A",1), (get_base1(A,row,(T - 1),"A",1) * get_base1(get_base1(p,2,"p",1),T1,"p",2)));
                    } else {

                        current_statement_begin__ = 36;
                        stan::math::assign(get_base1_lhs(A,row,T,"A",1), ((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)));
                    }
                }
            }
        }
        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_observed_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_observed(k, T0, T1, p, s, f, pstream__);
    }
};

template <typename T3__, typename T4__, typename T5__>
typename boost::math::tools::promote_args<T3__, T4__, T5__>::type
impute_unobserved(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 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__ = 51;
        fun_scalar_t__ prob;
        (void) prob;  // dummy to suppress unused var warning

        stan::math::initialize(prob, std::numeric_limits<double>::quiet_NaN());
        stan::math::fill(prob,DUMMY_VAR__);
        current_statement_begin__ = 52;
        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) + 1));
        current_statement_begin__ = 53;
        validate_non_negative_index("A", "T", T);
        validate_non_negative_index("A", "T", T);
        Eigen::Matrix<fun_scalar_t__,Eigen::Dynamic,Eigen::Dynamic>  A(static_cast<Eigen::VectorXd::Index>(T),static_cast<Eigen::VectorXd::Index>(T));
        (void) A;  // dummy to suppress unused var warning

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


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

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

                current_statement_begin__ = 56;
                stan::math::assign(get_base1_lhs(A,row,1,"A",1), ((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))));
            } else {

                current_statement_begin__ = 60;
                stan::math::assign(get_base1_lhs(A,row,1,"A",1), ((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))));
            }
            current_statement_begin__ = 63;
            for (int col = 2; col <= (T - 1); ++col) {

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

                    current_statement_begin__ = 65;
                    stan::math::assign(get_base1_lhs(A,row,col,"A",1), (get_base1(A,row,(col - 1),"A",1) * (1 - get_base1(get_base1(p,2,"p",1),(T0 + col),"p",2))));
                } else if (as_bool(logical_eq(col,row))) {

                    current_statement_begin__ = 68;
                    stan::math::assign(get_base1_lhs(A,row,col,"A",1), ((get_base1(A,row,(col - 1),"A",1) * (1 - get_base1(get_base1(s,1,"s",1),(T0 + (col - 1)),"s",2))) * (1 - get_base1(get_base1(p,2,"p",1),(T0 + col),"p",2))));
                } else {

                    current_statement_begin__ = 72;
                    stan::math::assign(get_base1_lhs(A,row,col,"A",1), (((get_base1(A,row,(col - 1),"A",1) * get_base1(get_base1(s,1,"s",1),(T0 + (col - 1)),"s",2)) * (1 - get_base1(get_base1(f,1,"f",1),(T0 + (col - 1)),"f",2))) * (1 - get_base1(get_base1(p,1,"p",1),(T0 + col),"p",2))));
                }
            }
            current_statement_begin__ = 78;
            if (as_bool(logical_lt(row,T))) {

                current_statement_begin__ = 79;
                stan::math::assign(get_base1_lhs(A,row,T,"A",1), (get_base1(A,row,(T - 1),"A",1) * (1 - get_base1(get_base1(p,2,"p",1),T1,"p",2))));
            } else {

                current_statement_begin__ = 82;
                stan::math::assign(get_base1_lhs(A,row,T,"A",1), ((get_base1(A,row,(T - 1),"A",1) * (1 - get_base1(get_base1(s,1,"s",1),(T1 - 1),"s",2))) * (1 - get_base1(get_base1(p,2,"p",1),T1,"p",2))));
            }
        }
        current_statement_begin__ = 88;
        stan::math::assign(prob, 0.0);
        current_statement_begin__ = 89;
        for (int i = 1; i <= T; ++i) {
            current_statement_begin__ = 90;
            stan::math::assign(prob, (prob + get_base1(A,i,T,"A",1)));
        }
        current_statement_begin__ = 91;
        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_unobserved_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_unobserved(k, T0, T1, p, s, f, pstream__);
    }
};

int
first_obs(const std::vector<int>& y_i, 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__ = 95;
        for (int k = 1; k <= size(y_i); ++k) {

            current_statement_begin__ = 96;
            if (as_bool(get_base1(y_i,k,"y_i",1))) {
                current_statement_begin__ = 96;
                return stan::math::promote_scalar<fun_return_scalar_t__>(k);
            }
        }
        current_statement_begin__ = 98;
        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 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__ = 102;
        for (int k = 1; k <= size(y_i); ++k) {

            current_statement_begin__ = 103;
            if (as_bool(logical_eq(get_base1(y_i,k,"y_i",1),3))) {
                current_statement_begin__ = 103;
                return stan::math::promote_scalar<fun_return_scalar_t__>(k);
            }
        }
        current_statement_begin__ = 105;
        return stan::math::promote_scalar<fun_return_scalar_t__>(size(y_i));
    } 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 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__ = 109;
        int k_back(0);
        (void) k_back;  // dummy to suppress unused var warning

        stan::math::fill(k_back, std::numeric_limits<int>::min());


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

            current_statement_begin__ = 111;
            stan::math::assign(k_back, (start_ind - k));
            current_statement_begin__ = 112;
            if (as_bool(get_base1(y_i,k_back,"y_i",1))) {

                current_statement_begin__ = 113;
                return stan::math::promote_scalar<fun_return_scalar_t__>(k_back);
            }
        }
        current_statement_begin__ = 116;
        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 modeld3143b1c2538_mark_simple_ends : public prob_grad {
private:
    int T_years;
    int N_plant;
    vector<vector<int> > y;
public:
    modeld3143b1c2538_mark_simple_ends(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, 0, pstream__);
    }

    modeld3143b1c2538_mark_simple_ends(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__ = "modeld3143b1c2538_mark_simple_ends_namespace::modeld3143b1c2538_mark_simple_ends";
        (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__ = 121;
            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__++];
            current_statement_begin__ = 122;
            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__++];
            current_statement_begin__ = 123;
            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));
            validate_non_negative_index("y", "N_plant", N_plant);
            validate_non_negative_index("y", "T_years", 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_limit_1__ = T_years;
            for (size_t i_1__ = 0; i_1__ < y_limit_1__; ++i_1__) {
                size_t y_limit_0__ = N_plant;
                for (size_t i_0__ = 0; i_0__ < y_limit_0__; ++i_0__) {
                    y[i_0__][i_1__] = vals_i__[pos__++];
                }
            }

            // validate, data variables
            current_statement_begin__ = 121;
            check_greater_or_equal(function__,"T_years",T_years,0);
            current_statement_begin__ = 122;
            check_greater_or_equal(function__,"N_plant",N_plant,0);
            current_statement_begin__ = 123;
            for (int k0__ = 0; k0__ < N_plant; ++k0__) {
                for (int k1__ = 0; k1__ < T_years; ++k1__) {
                    check_greater_or_equal(function__,"y[k0__][k1__]",y[k0__][k1__],0);
                    check_less_or_equal(function__,"y[k0__][k1__]",y[k0__][k1__],3);
                }
            }
            // initialize data variables


            // validate transformed data

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

    ~modeld3143b1c2538_mark_simple_ends() { }


    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("p")))
            throw std::runtime_error("variable p missing");
        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("initialization", "p", "double", context__.to_vec(2,T_years));
        std::vector<std::vector<double> > p(2,std::vector<double>(T_years,double(0)));
        for (int i1__ = 0U; i1__ < T_years; ++i1__)
            for (int i0__ = 0U; i0__ < 2; ++i0__)
                p[i0__][i1__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < 2; ++i0__)
            for (int i1__ = 0U; i1__ < T_years; ++i1__)
                try {
            writer__.scalar_lub_unconstrain(0,1,p[i0__][i1__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable p: ") + e.what());
        }

        if (!(context__.contains_r("s")))
            throw std::runtime_error("variable s missing");
        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("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)));
        for (int i1__ = 0U; i1__ < (T_years - 1); ++i1__)
            for (int i0__ = 0U; i0__ < 2; ++i0__)
                s[i0__][i1__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < 2; ++i0__)
            for (int i1__ = 0U; i1__ < (T_years - 1); ++i1__)
                try {
            writer__.scalar_lub_unconstrain(0,1,s[i0__][i1__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable s: ") + e.what());
        }

        if (!(context__.contains_r("f")))
            throw std::runtime_error("variable f missing");
        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("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)));
        for (int i1__ = 0U; i1__ < (T_years - 1); ++i1__)
            for (int i0__ = 0U; i0__ < 2; ++i0__)
                f[i0__][i1__] = vals_r__[pos__++];
        for (int i0__ = 0U; i0__ < 2; ++i0__)
            for (int i1__ = 0U; i1__ < (T_years - 1); ++i1__)
                try {
            writer__.scalar_lub_unconstrain(0,1,f[i0__][i1__]);
        } catch (const std::exception& e) { 
            throw std::runtime_error(std::string("Error transforming variable f: ") + 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<vector<T__> > p;
            size_t dim_p_0__ = 2;
            p.resize(dim_p_0__);
            for (size_t k_0__ = 0; k_0__ < dim_p_0__; ++k_0__) {
                size_t dim_p_1__ = T_years;
                p[k_0__].reserve(dim_p_1__);
                for (size_t k_1__ = 0; k_1__ < dim_p_1__; ++k_1__) {
                    if (jacobian__)
                        p[k_0__].push_back(in__.scalar_lub_constrain(0,1,lp__));
                    else
                        p[k_0__].push_back(in__.scalar_lub_constrain(0,1));
                }
            }

            vector<vector<T__> > s;
            size_t dim_s_0__ = 2;
            s.resize(dim_s_0__);
            for (size_t k_0__ = 0; k_0__ < dim_s_0__; ++k_0__) {
                size_t dim_s_1__ = (T_years - 1);
                s[k_0__].reserve(dim_s_1__);
                for (size_t k_1__ = 0; k_1__ < dim_s_1__; ++k_1__) {
                    if (jacobian__)
                        s[k_0__].push_back(in__.scalar_lub_constrain(0,1,lp__));
                    else
                        s[k_0__].push_back(in__.scalar_lub_constrain(0,1));
                }
            }

            vector<vector<T__> > f;
            size_t dim_f_0__ = 2;
            f.resize(dim_f_0__);
            for (size_t k_0__ = 0; k_0__ < dim_f_0__; ++k_0__) {
                size_t dim_f_1__ = (T_years - 1);
                f[k_0__].reserve(dim_f_1__);
                for (size_t k_1__ = 0; k_1__ < dim_f_1__; ++k_1__) {
                    if (jacobian__)
                        f[k_0__].push_back(in__.scalar_lub_constrain(0,1,lp__));
                    else
                        f[k_0__].push_back(in__.scalar_lub_constrain(0,1));
                }
            }


            // transformed parameters



            // validate transformed parameters

            const char* function__ = "validate transformed params";
            (void) function__;  // dummy to suppress unused var warning

            // model body

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

                    current_statement_begin__ = 145;
                    lp_accum__.add(beta_log<propto__>(get_base1(get_base1(s,i,"s",1),j,"s",2), 12, 2));
                    current_statement_begin__ = 146;
                    lp_accum__.add(beta_log<propto__>(get_base1(get_base1(f,i,"f",1),j,"f",2), 2, 4));
                    current_statement_begin__ = 147;
                    lp_accum__.add(beta_log<propto__>(get_base1(get_base1(p,i,"p",1),j,"p",2), 0.75, 0.80000000000000004));
                }
            }
            current_statement_begin__ = 150;
            for (int i = 1; i <= N_plant; ++i) {
                {
                current_statement_begin__ = 151;
                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__ = 152;
                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__ = 153;
                stan::math::assign(first_year, first_obs(get_base1(y,i,"y",1), pstream__));
                current_statement_begin__ = 154;
                stan::math::assign(last_year, last_obs(get_base1(y,i,"y",1),i, pstream__));
                current_statement_begin__ = 158;
                if (as_bool((primitive_value(logical_eq(get_base1(get_base1(y,i,"y",1),last_year,"y",2),1)) || primitive_value(logical_eq(get_base1(get_base1(y,i,"y",1),last_year,"y",2),2))))) {

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

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

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

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

                                current_statement_begin__ = 175;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f,1,"f",1),(t - 1),"f",2)));
                                current_statement_begin__ = 177;
                                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__ = 183;
                            if (as_bool(logical_eq(get_base1(get_base1(y,i,"y",1),(t - 1),"y",2),2))) {

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

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

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

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

                                current_statement_begin__ = 210;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f,1,"f",1),(t - 1),"f",2)));
                                current_statement_begin__ = 212;
                                lp_accum__.add(bernoulli_log<propto__>(1, get_base1(get_base1(s,1,"s",1),(t - 1),"s",2)));
                            }
                        }
                    }
                } else {
                    {
                    current_statement_begin__ = 221;
                    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__ = 223;
                    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__ = 224;
                    T__ prob_combined;
                    (void) prob_combined;  // dummy to suppress unused var warning

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


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

                        current_statement_begin__ = 227;
                        stan::math::assign(prob_combined, impute_observed(k_ind,block_start,last_year,p,s,f, pstream__));
                    } else {

                        current_statement_begin__ = 231;
                        stan::math::assign(prob_combined, impute_unobserved(k_ind,block_start,last_year,p,s,f, pstream__));
                    }
                    current_statement_begin__ = 234;
                    lp_accum__.add(bernoulli_log<propto__>(1, prob_combined));
                    current_statement_begin__ = 236;
                    for (int t = (first_year + 1); t <= block_start; ++t) {

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

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

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

                                current_statement_begin__ = 251;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f,1,"f",1),(t - 1),"f",2)));
                                current_statement_begin__ = 253;
                                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__ = 259;
                            if (as_bool(logical_eq(get_base1(get_base1(y,i,"y",1),(t - 1),"y",2),2))) {

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

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

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

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

                                current_statement_begin__ = 284;
                                lp_accum__.add(bernoulli_log<propto__>(0, get_base1(get_base1(f,1,"f",1),(t - 1),"f",2)));
                                current_statement_begin__ = 286;
                                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 {
        vars__.resize(0);
        stan::io::reader<double> in__(params_r__,params_i__);
        static const char* function__ = "modeld3143b1c2538_mark_simple_ends_namespace::write_array";
        (void) function__;  // dummy to suppress unused var warning
        // read-transform, write parameters
        vector<vector<double> > p;
        size_t dim_p_0__ = 2;
        p.resize(dim_p_0__);
        for (size_t k_0__ = 0; k_0__ < dim_p_0__; ++k_0__) {
            size_t dim_p_1__ = T_years;
            for (size_t k_1__ = 0; k_1__ < dim_p_1__; ++k_1__) {
                p[k_0__].push_back(in__.scalar_lub_constrain(0,1));
            }
        }
        vector<vector<double> > s;
        size_t dim_s_0__ = 2;
        s.resize(dim_s_0__);
        for (size_t k_0__ = 0; k_0__ < dim_s_0__; ++k_0__) {
            size_t dim_s_1__ = (T_years - 1);
            for (size_t k_1__ = 0; k_1__ < dim_s_1__; ++k_1__) {
                s[k_0__].push_back(in__.scalar_lub_constrain(0,1));
            }
        }
        vector<vector<double> > f;
        size_t dim_f_0__ = 2;
        f.resize(dim_f_0__);
        for (size_t k_0__ = 0; k_0__ < dim_f_0__; ++k_0__) {
            size_t dim_f_1__ = (T_years - 1);
            for (size_t k_1__ = 0; k_1__ < dim_f_1__; ++k_1__) {
                f[k_0__].push_back(in__.scalar_lub_constrain(0,1));
            }
        }
            for (int k_1__ = 0; k_1__ < T_years; ++k_1__) {
                for (int k_0__ = 0; k_0__ < 2; ++k_0__) {
                vars__.push_back(p[k_0__][k_1__]);
                }
            }
            for (int k_1__ = 0; k_1__ < (T_years - 1); ++k_1__) {
                for (int k_0__ = 0; k_0__ < 2; ++k_0__) {
                vars__.push_back(s[k_0__][k_1__]);
                }
            }
            for (int k_1__ = 0; k_1__ < (T_years - 1); ++k_1__) {
                for (int k_0__ = 0; k_0__ < 2; ++k_0__) {
                vars__.push_back(f[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 {



            // validate transformed parameters

            // write transformed parameters

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



            // validate generated quantities

            // write generated quantities
        } 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 "modeld3143b1c2538_mark_simple_ends";
    }


    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__ <= T_years; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= 2; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "p" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_1__ = 1; k_1__ <= (T_years - 1); ++k_1__) {
            for (int k_0__ = 1; k_0__ <= 2; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "s" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_1__ = 1; k_1__ <= (T_years - 1); ++k_1__) {
            for (int k_0__ = 1; k_0__ <= 2; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "f" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__ && !include_tparams__) return;

        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__ <= T_years; ++k_1__) {
            for (int k_0__ = 1; k_0__ <= 2; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "p" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_1__ = 1; k_1__ <= (T_years - 1); ++k_1__) {
            for (int k_0__ = 1; k_0__ <= 2; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "s" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        for (int k_1__ = 1; k_1__ <= (T_years - 1); ++k_1__) {
            for (int k_0__ = 1; k_0__ <= 2; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "f" << '.' << k_0__ << '.' << k_1__;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (!include_gqs__) return;
    }

}; // model

}

typedef modeld3143b1c2538_mark_simple_ends_namespace::modeld3143b1c2538_mark_simple_ends stan_model;