// Code generated by Stan version 2.21.0

#include <stan/model/model_header.hpp>

namespace model_estimate_infections_namespace {

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

static int current_statement_begin__;

stan::io::program_reader prog_reader__() {
    stan::io::program_reader reader;
    reader.add_event(0, 0, "start", "model_estimate_infections");
    reader.add_event(712, 710, "end", "model_estimate_infections");
    return reader;
}

template <typename T1__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__>::type, Eigen::Dynamic, 1>
discretised_gamma_pmf(const std::vector<int>& y,
                          const T1__& mu,
                          const T2__& sigma,
                          const int& max_val, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__, T2__>::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__ = 5;
        int n(0);
        (void) n;  // dummy to suppress unused var warning
        stan::math::fill(n, std::numeric_limits<int>::min());
        stan::math::assign(n,num_elements(y));

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

        current_statement_begin__ = 7;
        local_scalar_t__ trunc_pmf(DUMMY_VAR__);
        (void) trunc_pmf;  // dummy to suppress unused var warning
        stan::math::initialize(trunc_pmf, DUMMY_VAR__);
        stan::math::fill(trunc_pmf, DUMMY_VAR__);

        current_statement_begin__ = 9;
        local_scalar_t__ small(DUMMY_VAR__);
        (void) small;  // dummy to suppress unused var warning
        stan::math::initialize(small, DUMMY_VAR__);
        stan::math::fill(small, DUMMY_VAR__);
        stan::math::assign(small,1e-5);

        current_statement_begin__ = 10;
        local_scalar_t__ large(DUMMY_VAR__);
        (void) large;  // dummy to suppress unused var warning
        stan::math::initialize(large, DUMMY_VAR__);
        stan::math::fill(large, DUMMY_VAR__);
        stan::math::assign(large,1e8);

        current_statement_begin__ = 11;
        local_scalar_t__ c_sigma(DUMMY_VAR__);
        (void) c_sigma;  // dummy to suppress unused var warning
        stan::math::initialize(c_sigma, DUMMY_VAR__);
        stan::math::fill(c_sigma, DUMMY_VAR__);
        stan::math::assign(c_sigma,stan::math::fmax(small, sigma));

        current_statement_begin__ = 12;
        local_scalar_t__ c_mu(DUMMY_VAR__);
        (void) c_mu;  // dummy to suppress unused var warning
        stan::math::initialize(c_mu, DUMMY_VAR__);
        stan::math::fill(c_mu, DUMMY_VAR__);
        stan::math::assign(c_mu,stan::math::fmax(small, mu));

        current_statement_begin__ = 13;
        local_scalar_t__ alpha(DUMMY_VAR__);
        (void) alpha;  // dummy to suppress unused var warning
        stan::math::initialize(alpha, DUMMY_VAR__);
        stan::math::fill(alpha, DUMMY_VAR__);
        stan::math::assign(alpha,pow((c_mu / c_sigma), 2));

        current_statement_begin__ = 14;
        local_scalar_t__ beta(DUMMY_VAR__);
        (void) beta;  // dummy to suppress unused var warning
        stan::math::initialize(beta, DUMMY_VAR__);
        stan::math::fill(beta, DUMMY_VAR__);
        stan::math::assign(beta,(c_mu / pow(c_sigma, 2)));


        current_statement_begin__ = 16;
        stan::math::assign(alpha, stan::math::fmax(small, alpha));
        current_statement_begin__ = 17;
        stan::math::assign(alpha, stan::math::fmin(large, alpha));
        current_statement_begin__ = 18;
        stan::math::assign(beta, stan::math::fmax(small, beta));
        current_statement_begin__ = 19;
        stan::math::assign(beta, stan::math::fmin(large, beta));
        current_statement_begin__ = 21;
        stan::math::assign(trunc_pmf, (gamma_cdf((max_val + 1), alpha, beta) - gamma_cdf(1, alpha, beta)));
        current_statement_begin__ = 22;
        for (int i = 1; i <= n; ++i) {

            current_statement_begin__ = 23;
            stan::model::assign(pmf, 
                        stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                        ((gamma_cdf((get_base1(y, i, "y", 1) + 1), alpha, beta) - gamma_cdf(get_base1(y, i, "y", 1), alpha, beta)) / trunc_pmf), 
                        "assigning variable pmf");
        }
        current_statement_begin__ = 26;
        return stan::math::promote_scalar<fun_return_scalar_t__>(pmf);
        }
    } 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 discretised_gamma_pmf_functor__ {
    template <typename T1__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__>::type, Eigen::Dynamic, 1>
    operator()(const std::vector<int>& y,
                          const T1__& mu,
                          const T2__& sigma,
                          const int& max_val, std::ostream* pstream__) const {
        return discretised_gamma_pmf(y, mu, sigma, max_val, pstream__);
    }
};

template <typename T1__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__>::type, Eigen::Dynamic, 1>
discretised_lognormal_pmf(const std::vector<int>& y,
                              const T1__& mu,
                              const T2__& sigma,
                              const int& max_val, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__, T2__>::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__ = 31;
        int n(0);
        (void) n;  // dummy to suppress unused var warning
        stan::math::fill(n, std::numeric_limits<int>::min());
        stan::math::assign(n,num_elements(y));

        current_statement_begin__ = 32;
        validate_non_negative_index("pmf", "n", n);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> pmf(n);
        stan::math::initialize(pmf, DUMMY_VAR__);
        stan::math::fill(pmf, DUMMY_VAR__);

        current_statement_begin__ = 33;
        local_scalar_t__ small(DUMMY_VAR__);
        (void) small;  // dummy to suppress unused var warning
        stan::math::initialize(small, DUMMY_VAR__);
        stan::math::fill(small, DUMMY_VAR__);
        stan::math::assign(small,1e-5);

        current_statement_begin__ = 34;
        local_scalar_t__ c_sigma(DUMMY_VAR__);
        (void) c_sigma;  // dummy to suppress unused var warning
        stan::math::initialize(c_sigma, DUMMY_VAR__);
        stan::math::fill(c_sigma, DUMMY_VAR__);
        stan::math::assign(c_sigma,stan::math::fmax(small, sigma));

        current_statement_begin__ = 35;
        local_scalar_t__ c_mu(DUMMY_VAR__);
        (void) c_mu;  // dummy to suppress unused var warning
        stan::math::initialize(c_mu, DUMMY_VAR__);
        stan::math::fill(c_mu, DUMMY_VAR__);
        stan::math::assign(c_mu,stan::math::fmax(small, mu));

        current_statement_begin__ = 36;
        validate_non_negative_index("adj_y", "n", n);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> adj_y(n);
        stan::math::initialize(adj_y, DUMMY_VAR__);
        stan::math::fill(adj_y, DUMMY_VAR__);
        stan::math::assign(adj_y,add(to_vector(y), small));

        current_statement_begin__ = 37;
        validate_non_negative_index("upper_y", "n", n);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> upper_y(n);
        stan::math::initialize(upper_y, DUMMY_VAR__);
        stan::math::fill(upper_y, DUMMY_VAR__);
        stan::math::assign(upper_y,divide(subtract(stan::math::log(add(adj_y, 1)), c_mu), c_sigma));

        current_statement_begin__ = 38;
        validate_non_negative_index("lower_y", "n", n);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> lower_y(n);
        stan::math::initialize(lower_y, DUMMY_VAR__);
        stan::math::fill(lower_y, DUMMY_VAR__);
        stan::math::assign(lower_y,divide(subtract(stan::math::log(adj_y), c_mu), c_sigma));

        current_statement_begin__ = 39;
        local_scalar_t__ max_cdf(DUMMY_VAR__);
        (void) max_cdf;  // dummy to suppress unused var warning
        stan::math::initialize(max_cdf, DUMMY_VAR__);
        stan::math::fill(max_cdf, DUMMY_VAR__);
        stan::math::assign(max_cdf,normal_cdf(((stan::math::log((max_val + small)) - c_mu) / c_sigma), 0.0, 1.0));

        current_statement_begin__ = 40;
        local_scalar_t__ min_cdf(DUMMY_VAR__);
        (void) min_cdf;  // dummy to suppress unused var warning
        stan::math::initialize(min_cdf, DUMMY_VAR__);
        stan::math::fill(min_cdf, DUMMY_VAR__);
        stan::math::assign(min_cdf,normal_cdf(((stan::math::log(small) - c_mu) / c_sigma), 0.0, 1.0));

        current_statement_begin__ = 41;
        local_scalar_t__ trunc_cdf(DUMMY_VAR__);
        (void) trunc_cdf;  // dummy to suppress unused var warning
        stan::math::initialize(trunc_cdf, DUMMY_VAR__);
        stan::math::fill(trunc_cdf, DUMMY_VAR__);
        stan::math::assign(trunc_cdf,(max_cdf - min_cdf));


        current_statement_begin__ = 42;
        for (int i = 1; i <= n; ++i) {

            current_statement_begin__ = 43;
            stan::model::assign(pmf, 
                        stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                        ((normal_cdf(get_base1(upper_y, i, "upper_y", 1), 0.0, 1.0) - normal_cdf(get_base1(lower_y, i, "lower_y", 1), 0.0, 1.0)) / trunc_cdf), 
                        "assigning variable pmf");
        }
        current_statement_begin__ = 46;
        return stan::math::promote_scalar<fun_return_scalar_t__>(pmf);
        }
    } 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 discretised_lognormal_pmf_functor__ {
    template <typename T1__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__>::type, Eigen::Dynamic, 1>
    operator()(const std::vector<int>& y,
                              const T1__& mu,
                              const T2__& sigma,
                              const int& max_val, std::ostream* pstream__) const {
        return discretised_lognormal_pmf(y, mu, sigma, max_val, pstream__);
    }
};

template <typename T0__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__>::type, Eigen::Dynamic, 1>
reverse_mf(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& pmf,
               const int& max_pmf, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__>::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__ = 51;
        validate_non_negative_index("rev_pmf", "max_pmf", max_pmf);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> rev_pmf(max_pmf);
        stan::math::initialize(rev_pmf, DUMMY_VAR__);
        stan::math::fill(rev_pmf, DUMMY_VAR__);


        current_statement_begin__ = 52;
        for (int d = 1; d <= max_pmf; ++d) {

            current_statement_begin__ = 53;
            stan::model::assign(rev_pmf, 
                        stan::model::cons_list(stan::model::index_uni(d), stan::model::nil_index_list()), 
                        get_base1(pmf, ((max_pmf - d) + 1), "pmf", 1), 
                        "assigning variable rev_pmf");
        }
        current_statement_begin__ = 55;
        return stan::math::promote_scalar<fun_return_scalar_t__>(rev_pmf);
        }
    } 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 reverse_mf_functor__ {
    template <typename T0__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& pmf,
               const int& max_pmf, std::ostream* pstream__) const {
        return reverse_mf(pmf, max_pmf, pstream__);
    }
};

template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
convolve(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& cases,
             const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& rev_pmf, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::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__ = 60;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(cases));

        current_statement_begin__ = 61;
        int max_pmf(0);
        (void) max_pmf;  // dummy to suppress unused var warning
        stan::math::fill(max_pmf, std::numeric_limits<int>::min());
        stan::math::assign(max_pmf,num_elements(rev_pmf));

        current_statement_begin__ = 62;
        validate_non_negative_index("conv_cases", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> conv_cases(t);
        stan::math::initialize(conv_cases, DUMMY_VAR__);
        stan::math::fill(conv_cases, DUMMY_VAR__);
        stan::math::assign(conv_cases,rep_vector(1e-5, t));


        current_statement_begin__ = 63;
        for (int s = 1; s <= t; ++s) {

            current_statement_begin__ = 64;
            stan::model::assign(conv_cases, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        (stan::model::rvalue(conv_cases, stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), "conv_cases") + dot_product(stan::model::rvalue(cases, stan::model::cons_list(stan::model::index_min_max(std::max(1, ((s - max_pmf) + 1)), s), stan::model::nil_index_list()), "cases"), tail(rev_pmf, std::min(max_pmf, s)))), 
                        "assigning variable conv_cases");
        }
        current_statement_begin__ = 67;
        return stan::math::promote_scalar<fun_return_scalar_t__>(conv_cases);
        }
    } 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 convolve_functor__ {
    template <typename T0__, typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& cases,
             const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& rev_pmf, std::ostream* pstream__) const {
        return convolve(cases, rev_pmf, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__>::type, Eigen::Dynamic, 1>
convolve_to_report(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& infections,
                       const std::vector<T1__>& delay_mean,
                       const std::vector<T2__>& delay_sd,
                       const std::vector<int>& max_delay,
                       const int& seeding_time, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__>::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__ = 77;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(infections));

        current_statement_begin__ = 78;
        validate_non_negative_index("reports", "(t - seeding_time)", (t - seeding_time));
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> reports((t - seeding_time));
        stan::math::initialize(reports, DUMMY_VAR__);
        stan::math::fill(reports, DUMMY_VAR__);

        current_statement_begin__ = 79;
        validate_non_negative_index("unobs_reports", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> unobs_reports(t);
        stan::math::initialize(unobs_reports, DUMMY_VAR__);
        stan::math::fill(unobs_reports, DUMMY_VAR__);
        stan::math::assign(unobs_reports,infections);

        current_statement_begin__ = 80;
        int delays(0);
        (void) delays;  // dummy to suppress unused var warning
        stan::math::fill(delays, std::numeric_limits<int>::min());
        stan::math::assign(delays,num_elements(delay_mean));


        current_statement_begin__ = 81;
        if (as_bool(delays)) {

            current_statement_begin__ = 82;
            for (int s = 1; s <= delays; ++s) {
                {
                current_statement_begin__ = 83;
                validate_non_negative_index("pmf", "get_base1(max_delay, s, \"max_delay\", 1)", get_base1(max_delay, s, "max_delay", 1));
                Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> pmf(get_base1(max_delay, s, "max_delay", 1));
                stan::math::initialize(pmf, DUMMY_VAR__);
                stan::math::fill(pmf, DUMMY_VAR__);
                stan::math::assign(pmf,rep_vector(1e-5, get_base1(max_delay, s, "max_delay", 1)));

                current_statement_begin__ = 84;
                validate_non_negative_index("delay_indexes", "get_base1(max_delay, s, \"max_delay\", 1)", get_base1(max_delay, s, "max_delay", 1));
                std::vector<int  > delay_indexes(get_base1(max_delay, s, "max_delay", 1), int(0));
                stan::math::fill(delay_indexes, std::numeric_limits<int>::min());


                current_statement_begin__ = 85;
                for (int i = 1; i <= get_base1(max_delay, s, "max_delay", 1); ++i) {

                    current_statement_begin__ = 86;
                    stan::model::assign(delay_indexes, 
                                stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                                (get_base1(max_delay, s, "max_delay", 1) - i), 
                                "assigning variable delay_indexes");
                }
                current_statement_begin__ = 88;
                stan::math::assign(pmf, add(pmf, discretised_lognormal_pmf(delay_indexes, get_base1(delay_mean, s, "delay_mean", 1), get_base1(delay_sd, s, "delay_sd", 1), get_base1(max_delay, s, "max_delay", 1), pstream__)));
                current_statement_begin__ = 90;
                stan::math::assign(unobs_reports, convolve(unobs_reports, pmf, pstream__));
                }
            }
            current_statement_begin__ = 92;
            stan::math::assign(reports, stan::model::rvalue(unobs_reports, stan::model::cons_list(stan::model::index_min_max((seeding_time + 1), t), stan::model::nil_index_list()), "unobs_reports"));
        } else {

            current_statement_begin__ = 94;
            stan::math::assign(reports, stan::model::rvalue(infections, stan::model::cons_list(stan::model::index_min_max((seeding_time + 1), t), stan::model::nil_index_list()), "infections"));
        }
        current_statement_begin__ = 96;
        return stan::math::promote_scalar<fun_return_scalar_t__>(reports);
        }
    } 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 convolve_to_report_functor__ {
    template <typename T0__, typename T1__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& infections,
                       const std::vector<T1__>& delay_mean,
                       const std::vector<T2__>& delay_sd,
                       const std::vector<int>& max_delay,
                       const int& seeding_time, std::ostream* pstream__) const {
        return convolve_to_report(infections, delay_mean, delay_sd, max_delay, seeding_time, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T_lp__, typename T_lp_accum__>
void
delays_lp(const std::vector<T0__>& delay_mean,
              const std::vector<T1__>& delay_mean_mean,
              const std::vector<T2__>& delay_mean_sd,
              const std::vector<T3__>& delay_sd,
              const std::vector<T4__>& delay_sd_mean,
              const std::vector<T5__>& delay_sd_sd,
              const int& weight, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T_lp__>::type>::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__ = 101;
        int delays(0);
        (void) delays;  // dummy to suppress unused var warning
        stan::math::fill(delays, std::numeric_limits<int>::min());
        stan::math::assign(delays,num_elements(delay_mean));


        current_statement_begin__ = 102;
        if (as_bool(delays)) {

            current_statement_begin__ = 103;
            for (int s = 1; s <= delays; ++s) {

                current_statement_begin__ = 104;
                lp_accum__.add((normal_log(get_base1(delay_mean, s, "delay_mean", 1), get_base1(delay_mean_mean, s, "delay_mean_mean", 1), get_base1(delay_mean_sd, s, "delay_mean_sd", 1)) * weight));
                current_statement_begin__ = 105;
                lp_accum__.add((normal_log(get_base1(delay_sd, s, "delay_sd", 1), get_base1(delay_sd_mean, s, "delay_sd_mean", 1), get_base1(delay_sd_sd, s, "delay_sd_sd", 1)) * weight));
            }
        }
        }
    } 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 delays_lp_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const std::vector<T0__>& delay_mean,
              const std::vector<T1__>& delay_mean_mean,
              const std::vector<T2__>& delay_mean_sd,
              const std::vector<T3__>& delay_sd,
              const std::vector<T4__>& delay_sd_mean,
              const std::vector<T5__>& delay_sd_sd,
              const int& weight, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return delays_lp(delay_mean, delay_mean_mean, delay_mean_sd, delay_sd, delay_sd_mean, delay_sd_sd, weight, lp__, lp_accum__, pstream__);
    }
};

template <typename T0__>
typename boost::math::tools::promote_args<T0__>::type
lambda(const T0__& L,
           const int& m, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__>::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__ = 113;
        local_scalar_t__ lam(DUMMY_VAR__);
        (void) lam;  // dummy to suppress unused var warning
        stan::math::initialize(lam, DUMMY_VAR__);
        stan::math::fill(lam, DUMMY_VAR__);


        current_statement_begin__ = 114;
        stan::math::assign(lam, pow(((m * stan::math::pi()) / (2 * L)), 2));
        current_statement_begin__ = 115;
        return stan::math::promote_scalar<fun_return_scalar_t__>(lam);
        }
    } 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 lambda_functor__ {
    template <typename T0__>
        typename boost::math::tools::promote_args<T0__>::type
    operator()(const T0__& L,
           const int& m, std::ostream* pstream__) const {
        return lambda(L, m, pstream__);
    }
};

template <typename T0__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__>::type, Eigen::Dynamic, 1>
phi(const T0__& L,
        const int& m,
        const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& x, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T2__>::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__ = 120;
        validate_non_negative_index("fi", "rows(x)", rows(x));
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> fi(rows(x));
        stan::math::initialize(fi, DUMMY_VAR__);
        stan::math::fill(fi, DUMMY_VAR__);


        current_statement_begin__ = 121;
        stan::math::assign(fi, multiply((1 / stan::math::sqrt(L)), stan::math::sin(multiply(((m * stan::math::pi()) / (2 * L)), add(x, L)))));
        current_statement_begin__ = 122;
        return stan::math::promote_scalar<fun_return_scalar_t__>(fi);
        }
    } 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 phi_functor__ {
    template <typename T0__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__>::type, Eigen::Dynamic, 1>
    operator()(const T0__& L,
        const int& m,
        const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& x, std::ostream* pstream__) const {
        return phi(L, m, x, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__>
typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
spd_se(const T0__& alpha,
           const T1__& rho,
           const T2__& w, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__>::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__ = 126;
        local_scalar_t__ S(DUMMY_VAR__);
        (void) S;  // dummy to suppress unused var warning
        stan::math::initialize(S, DUMMY_VAR__);
        stan::math::fill(S, DUMMY_VAR__);


        current_statement_begin__ = 127;
        stan::math::assign(S, (((pow(alpha, 2) * stan::math::sqrt((2 * stan::math::pi()))) * rho) * stan::math::exp(((-(0.5) * pow(rho, 2)) * pow(w, 2)))));
        current_statement_begin__ = 128;
        return stan::math::promote_scalar<fun_return_scalar_t__>(S);
        }
    } catch (const std::exception& e) {
        stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
        // Next line prevents compiler griping about no return
        throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
    }
}


struct spd_se_functor__ {
    template <typename T0__, typename T1__, typename T2__>
        typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
    operator()(const T0__& alpha,
           const T1__& rho,
           const T2__& w, std::ostream* pstream__) const {
        return spd_se(alpha, rho, w, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__>
typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
spd_matern(const T0__& alpha,
               const T1__& rho,
               const T2__& w, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__>::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__ = 132;
        local_scalar_t__ S(DUMMY_VAR__);
        (void) S;  // dummy to suppress unused var warning
        stan::math::initialize(S, DUMMY_VAR__);
        stan::math::fill(S, DUMMY_VAR__);


        current_statement_begin__ = 133;
        stan::math::assign(S, ((((4 * pow(alpha, 2)) * pow((stan::math::sqrt(3) / rho), 3)) * 1) / pow((pow((stan::math::sqrt(3) / rho), 2) + pow(w, 2)), 2)));
        current_statement_begin__ = 134;
        return stan::math::promote_scalar<fun_return_scalar_t__>(S);
        }
    } catch (const std::exception& e) {
        stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
        // Next line prevents compiler griping about no return
        throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
    }
}


struct spd_matern_functor__ {
    template <typename T0__, typename T1__, typename T2__>
        typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
    operator()(const T0__& alpha,
               const T1__& rho,
               const T2__& w, std::ostream* pstream__) const {
        return spd_matern(alpha, rho, w, pstream__);
    }
};

int
setup_noise(const int& ot_h,
                const int& t,
                const int& horizon,
                const int& estimate_r,
                const int& stationary,
                const int& future_fixed,
                const int& fixed_from, 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__ = 139;
        int noise_time(0);
        (void) noise_time;  // dummy to suppress unused var warning
        stan::math::fill(noise_time, std::numeric_limits<int>::min());
        stan::math::assign(noise_time,(logical_gt(estimate_r, 0) ? (logical_gt(stationary, 0) ? ot_h : (ot_h - 1) ) : t ));

        current_statement_begin__ = 140;
        int noise_terms(0);
        (void) noise_terms;  // dummy to suppress unused var warning
        stan::math::fill(noise_terms, std::numeric_limits<int>::min());
        stan::math::assign(noise_terms,(logical_gt(future_fixed, 0) ? ((noise_time - horizon) + fixed_from) : noise_time ));


        current_statement_begin__ = 141;
        return stan::math::promote_scalar<fun_return_scalar_t__>(noise_terms);
        }
    } 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 setup_noise_functor__ {
            int
    operator()(const int& ot_h,
                const int& t,
                const int& horizon,
                const int& estimate_r,
                const int& stationary,
                const int& future_fixed,
                const int& fixed_from, std::ostream* pstream__) const {
        return setup_noise(ot_h, t, horizon, estimate_r, stationary, future_fixed, fixed_from, pstream__);
    }
};

template <typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__>::type, Eigen::Dynamic, Eigen::Dynamic>
setup_gp(const int& M,
             const T1__& L,
             const int& dimension, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__>::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__ = 145;
        validate_non_negative_index("time", "dimension", dimension);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> time(dimension);
        stan::math::initialize(time, DUMMY_VAR__);
        stan::math::fill(time, DUMMY_VAR__);

        current_statement_begin__ = 146;
        validate_non_negative_index("PHI", "dimension", dimension);
        validate_non_negative_index("PHI", "M", M);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> PHI(dimension, M);
        stan::math::initialize(PHI, DUMMY_VAR__);
        stan::math::fill(PHI, DUMMY_VAR__);

        current_statement_begin__ = 147;
        local_scalar_t__ half_dim(DUMMY_VAR__);
        (void) half_dim;  // dummy to suppress unused var warning
        stan::math::initialize(half_dim, DUMMY_VAR__);
        stan::math::fill(half_dim, DUMMY_VAR__);
        stan::math::assign(half_dim,(dimension / 2.0));


        current_statement_begin__ = 148;
        for (int s = 1; s <= dimension; ++s) {

            current_statement_begin__ = 149;
            stan::model::assign(time, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        ((s - half_dim) / half_dim), 
                        "assigning variable time");
        }
        current_statement_begin__ = 151;
        for (int m = 1; m <= M; ++m) {

            current_statement_begin__ = 152;
            stan::model::assign(PHI, 
                        stan::model::cons_list(stan::model::index_omni(), stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list())), 
                        phi(L, m, time, pstream__), 
                        "assigning variable PHI");
        }
        current_statement_begin__ = 154;
        return stan::math::promote_scalar<fun_return_scalar_t__>(PHI);
        }
    } 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 setup_gp_functor__ {
    template <typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__>::type, Eigen::Dynamic, Eigen::Dynamic>
    operator()(const int& M,
             const T1__& L,
             const int& dimension, std::ostream* pstream__) const {
        return setup_gp(M, L, dimension, pstream__);
    }
};

template <typename T0__, typename T2__, typename T3__, typename T4__, typename T5__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__, T3__, T4__, typename boost::math::tools::promote_args<T5__>::type>::type, Eigen::Dynamic, 1>
update_gp(const Eigen::Matrix<T0__, Eigen::Dynamic, Eigen::Dynamic>& PHI,
              const int& M,
              const T2__& L,
              const T3__& alpha,
              const T4__& rho,
              const Eigen::Matrix<T5__, Eigen::Dynamic, 1>& eta,
              const int& type, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T2__, T3__, T4__, typename boost::math::tools::promote_args<T5__>::type>::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__ = 159;
        validate_non_negative_index("diagSPD", "M", M);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> diagSPD(M);
        stan::math::initialize(diagSPD, DUMMY_VAR__);
        stan::math::fill(diagSPD, DUMMY_VAR__);

        current_statement_begin__ = 160;
        validate_non_negative_index("SPD_eta", "M", M);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> SPD_eta(M);
        stan::math::initialize(SPD_eta, DUMMY_VAR__);
        stan::math::fill(SPD_eta, DUMMY_VAR__);

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

        current_statement_begin__ = 162;
        validate_non_negative_index("noise", "noise_terms", noise_terms);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> noise(noise_terms);
        stan::math::initialize(noise, DUMMY_VAR__);
        stan::math::fill(noise, DUMMY_VAR__);
        stan::math::assign(noise,rep_vector(1e-6, noise_terms));

        current_statement_begin__ = 163;
        local_scalar_t__ unit_rho(DUMMY_VAR__);
        (void) unit_rho;  // dummy to suppress unused var warning
        stan::math::initialize(unit_rho, DUMMY_VAR__);
        stan::math::fill(unit_rho, DUMMY_VAR__);
        stan::math::assign(unit_rho,(rho / noise_terms));


        current_statement_begin__ = 165;
        if (as_bool(logical_eq(type, 0))) {

            current_statement_begin__ = 166;
            for (int m = 1; m <= M; ++m) {

                current_statement_begin__ = 167;
                stan::model::assign(diagSPD, 
                            stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list()), 
                            stan::math::sqrt(spd_se(alpha, unit_rho, stan::math::sqrt(lambda(L, m, pstream__)), pstream__)), 
                            "assigning variable diagSPD");
            }
        } else if (as_bool(logical_eq(type, 1))) {

            current_statement_begin__ = 170;
            for (int m = 1; m <= M; ++m) {

                current_statement_begin__ = 171;
                stan::model::assign(diagSPD, 
                            stan::model::cons_list(stan::model::index_uni(m), stan::model::nil_index_list()), 
                            stan::math::sqrt(spd_matern(alpha, unit_rho, stan::math::sqrt(lambda(L, m, pstream__)), pstream__)), 
                            "assigning variable diagSPD");
            }
        }
        current_statement_begin__ = 174;
        stan::math::assign(SPD_eta, elt_multiply(diagSPD, eta));
        current_statement_begin__ = 175;
        stan::math::assign(noise, add(noise, multiply(stan::model::rvalue(PHI, stan::model::cons_list(stan::model::index_omni(), stan::model::cons_list(stan::model::index_omni(), stan::model::nil_index_list())), "PHI"), SPD_eta)));
        current_statement_begin__ = 176;
        return stan::math::promote_scalar<fun_return_scalar_t__>(noise);
        }
    } 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 update_gp_functor__ {
    template <typename T0__, typename T2__, typename T3__, typename T4__, typename T5__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__, T3__, T4__, typename boost::math::tools::promote_args<T5__>::type>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, Eigen::Dynamic>& PHI,
              const int& M,
              const T2__& L,
              const T3__& alpha,
              const T4__& rho,
              const Eigen::Matrix<T5__, Eigen::Dynamic, 1>& eta,
              const int& type, std::ostream* pstream__) const {
        return update_gp(PHI, M, L, alpha, rho, eta, type, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T6__, typename T7__, typename T_lp__, typename T_lp_accum__>
void
gaussian_process_lp(const T0__& rho,
                        const T1__& alpha,
                        const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& eta,
                        const T3__& ls_meanlog,
                        const T4__& ls_sdlog,
                        const T5__& ls_min,
                        const T6__& ls_max,
                        const T7__& alpha_sd, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T6__, T7__, T_lp__>::type>::type 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__ = 182;
        lp_accum__.add(lognormal_log<propto__>(rho, ls_meanlog, ls_sdlog));
        if (rho < ls_min) lp_accum__.add(-std::numeric_limits<double>::infinity());
        else if (rho > ls_max) lp_accum__.add(-std::numeric_limits<double>::infinity());
        else lp_accum__.add(-log_diff_exp(lognormal_cdf_log(ls_max, ls_meanlog, ls_sdlog), lognormal_cdf_log(ls_min, ls_meanlog, ls_sdlog)));
        current_statement_begin__ = 183;
        lp_accum__.add(normal_log<propto__>(alpha, 0, alpha_sd));
        current_statement_begin__ = 184;
        lp_accum__.add(std_normal_log<propto__>(eta));
    } catch (const std::exception& e) {
        stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
        // Next line prevents compiler griping about no return
        throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
    }
}


struct gaussian_process_lp_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T6__, typename T7__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const T0__& rho,
                        const T1__& alpha,
                        const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& eta,
                        const T3__& ls_meanlog,
                        const T4__& ls_sdlog,
                        const T5__& ls_min,
                        const T6__& ls_max,
                        const T7__& alpha_sd, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return gaussian_process_lp(rho, alpha, eta, ls_meanlog, ls_sdlog, ls_min, ls_max, alpha_sd, lp__, lp_accum__, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T4__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__, T4__>::type, Eigen::Dynamic, 1>
update_Rt(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& input_R,
              const T1__& log_R,
              const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& noise,
              const std::vector<int>& bps,
              const std::vector<T4__>& bp_effects,
              const int& stationary, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T4__>::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__ = 191;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(input_R));

        current_statement_begin__ = 192;
        int bp_n(0);
        (void) bp_n;  // dummy to suppress unused var warning
        stan::math::fill(bp_n, std::numeric_limits<int>::min());
        stan::math::assign(bp_n,num_elements(bp_effects));

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

        current_statement_begin__ = 194;
        int gp_n(0);
        (void) gp_n;  // dummy to suppress unused var warning
        stan::math::fill(gp_n, std::numeric_limits<int>::min());
        stan::math::assign(gp_n,num_elements(noise));

        current_statement_begin__ = 196;
        validate_non_negative_index("bp", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> bp(t);
        stan::math::initialize(bp, DUMMY_VAR__);
        stan::math::fill(bp, DUMMY_VAR__);
        stan::math::assign(bp,rep_vector(0, t));

        current_statement_begin__ = 197;
        validate_non_negative_index("gp", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> gp(t);
        stan::math::initialize(gp, DUMMY_VAR__);
        stan::math::fill(gp, DUMMY_VAR__);
        stan::math::assign(gp,rep_vector(0, t));

        current_statement_begin__ = 198;
        validate_non_negative_index("R", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> R(t);
        stan::math::initialize(R, DUMMY_VAR__);
        stan::math::fill(R, DUMMY_VAR__);


        current_statement_begin__ = 200;
        if (as_bool(bp_n)) {

            current_statement_begin__ = 201;
            for (int s = 1; s <= t; ++s) {

                current_statement_begin__ = 202;
                if (as_bool(get_base1(bps, s, "bps", 1))) {

                    current_statement_begin__ = 203;
                    stan::math::assign(bp_c, (bp_c + get_base1(bps, s, "bps", 1)));
                    current_statement_begin__ = 204;
                    stan::model::assign(bp, 
                                stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                                get_base1(bp_effects, bp_c, "bp_effects", 1), 
                                "assigning variable bp");
                }
            }
            current_statement_begin__ = 207;
            stan::math::assign(bp, cumulative_sum(bp));
        }
        current_statement_begin__ = 210;
        if (as_bool(gp_n)) {

            current_statement_begin__ = 211;
            if (as_bool(stationary)) {

                current_statement_begin__ = 212;
                stan::model::assign(gp, 
                            stan::model::cons_list(stan::model::index_min_max(1, gp_n), stan::model::nil_index_list()), 
                            noise, 
                            "assigning variable gp");
                current_statement_begin__ = 214;
                if (as_bool(logical_gt(t, gp_n))) {

                    current_statement_begin__ = 215;
                    stan::model::assign(gp, 
                                stan::model::cons_list(stan::model::index_min_max((gp_n + 1), t), stan::model::nil_index_list()), 
                                rep_vector(get_base1(noise, gp_n, "noise", 1), (t - gp_n)), 
                                "assigning variable gp");
                }
            } else {

                current_statement_begin__ = 218;
                stan::model::assign(gp, 
                            stan::model::cons_list(stan::model::index_min_max(2, (gp_n + 1)), stan::model::nil_index_list()), 
                            noise, 
                            "assigning variable gp");
                current_statement_begin__ = 219;
                stan::math::assign(gp, cumulative_sum(gp));
            }
        }
        current_statement_begin__ = 223;
        stan::math::assign(R, add(add(rep_vector(log_R, t), bp), gp));
        current_statement_begin__ = 224;
        stan::math::assign(R, stan::math::exp(R));
        current_statement_begin__ = 225;
        return stan::math::promote_scalar<fun_return_scalar_t__>(R);
        }
    } 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 update_Rt_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T4__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__, T4__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& input_R,
              const T1__& log_R,
              const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& noise,
              const std::vector<int>& bps,
              const std::vector<T4__>& bp_effects,
              const int& stationary, std::ostream* pstream__) const {
        return update_Rt(input_R, log_R, noise, bps, bp_effects, stationary, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T7__, typename T8__, typename T9__, typename T10__, typename T_lp__, typename T_lp_accum__>
void
rt_lp(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& log_R,
          const std::vector<T1__>& initial_infections,
          const std::vector<T2__>& initial_growth,
          const std::vector<T3__>& bp_effects,
          const std::vector<T4__>& bp_sd,
          const int& bp_n,
          const int& seeding_time,
          const T7__& r_logmean,
          const T8__& r_logsd,
          const T9__& prior_infections,
          const T10__& prior_growth, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T7__, T8__, T9__, typename boost::math::tools::promote_args<T10__, T_lp__>::type>::type>::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__ = 233;
        lp_accum__.add(normal_log<propto__>(log_R, r_logmean, r_logsd));
        current_statement_begin__ = 235;
        if (as_bool(logical_gt(bp_n, 0))) {

            current_statement_begin__ = 236;
            lp_accum__.add(normal_log<propto__>(get_base1(bp_sd, 1, "bp_sd", 1), 0, 0.1));
            if (get_base1(bp_sd, 1, "bp_sd", 1) < 0) lp_accum__.add(-std::numeric_limits<double>::infinity());
            else lp_accum__.add(-normal_ccdf_log(0, 0, 0.1));
            current_statement_begin__ = 237;
            lp_accum__.add(normal_log<propto__>(bp_effects, 0, get_base1(bp_sd, 1, "bp_sd", 1)));
        }
        current_statement_begin__ = 240;
        lp_accum__.add(normal_log<propto__>(initial_infections, prior_infections, 0.2));
        current_statement_begin__ = 241;
        if (as_bool(logical_gt(seeding_time, 1))) {

            current_statement_begin__ = 242;
            lp_accum__.add(normal_log<propto__>(initial_growth, prior_growth, 0.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 ***");
    }
}


struct rt_lp_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T7__, typename T8__, typename T9__, typename T10__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& log_R,
          const std::vector<T1__>& initial_infections,
          const std::vector<T2__>& initial_growth,
          const std::vector<T3__>& bp_effects,
          const std::vector<T4__>& bp_sd,
          const int& bp_n,
          const int& seeding_time,
          const T7__& r_logmean,
          const T8__& r_logsd,
          const T9__& prior_infections,
          const T10__& prior_growth, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return rt_lp(log_R, initial_infections, initial_growth, bp_effects, bp_sd, bp_n, seeding_time, r_logmean, r_logsd, prior_infections, prior_growth, lp__, lp_accum__, pstream__);
    }
};

template <typename T0__, typename T1__>
typename boost::math::tools::promote_args<T0__, T1__>::type
update_infectiousness(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& infections,
                          const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& gt_pmf,
                          const int& seeding_time,
                          const int& max_gt,
                          const int& index, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::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__ = 250;
        int inf_start(0);
        (void) inf_start;  // dummy to suppress unused var warning
        stan::math::fill(inf_start, std::numeric_limits<int>::min());
        stan::math::assign(inf_start,std::max(1, ((index + seeding_time) - max_gt)));

        current_statement_begin__ = 251;
        int inf_end(0);
        (void) inf_end;  // dummy to suppress unused var warning
        stan::math::fill(inf_end, std::numeric_limits<int>::min());
        stan::math::assign(inf_end,((index + seeding_time) - 1));

        current_statement_begin__ = 252;
        int pmf_accessed(0);
        (void) pmf_accessed;  // dummy to suppress unused var warning
        stan::math::fill(pmf_accessed, std::numeric_limits<int>::min());
        stan::math::assign(pmf_accessed,std::min(max_gt, ((index + seeding_time) - 1)));

        current_statement_begin__ = 253;
        local_scalar_t__ new_inf(DUMMY_VAR__);
        (void) new_inf;  // dummy to suppress unused var warning
        stan::math::initialize(new_inf, DUMMY_VAR__);
        stan::math::fill(new_inf, DUMMY_VAR__);
        stan::math::assign(new_inf,dot_product(stan::model::rvalue(infections, stan::model::cons_list(stan::model::index_min_max(inf_start, inf_end), stan::model::nil_index_list()), "infections"), tail(gt_pmf, pmf_accessed)));


        current_statement_begin__ = 254;
        return stan::math::promote_scalar<fun_return_scalar_t__>(new_inf);
        }
    } 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 update_infectiousness_functor__ {
    template <typename T0__, typename T1__>
        typename boost::math::tools::promote_args<T0__, T1__>::type
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& infections,
                          const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& gt_pmf,
                          const int& seeding_time,
                          const int& max_gt,
                          const int& index, std::ostream* pstream__) const {
        return update_infectiousness(infections, gt_pmf, seeding_time, max_gt, index, pstream__);
    }
};

template <typename T0__, typename T2__, typename T3__, typename T5__, typename T6__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__, T3__, T5__, typename boost::math::tools::promote_args<T6__>::type>::type, Eigen::Dynamic, 1>
generate_infections(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& oR,
                        const int& uot,
                        const std::vector<T2__>& gt_mean,
                        const std::vector<T3__>& gt_sd,
                        const int& max_gt,
                        const std::vector<T5__>& initial_infections,
                        const std::vector<T6__>& initial_growth,
                        const int& pop,
                        const int& ht, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T2__, T3__, T5__, typename boost::math::tools::promote_args<T6__>::type>::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__ = 262;
        int ot(0);
        (void) ot;  // dummy to suppress unused var warning
        stan::math::fill(ot, std::numeric_limits<int>::min());
        stan::math::assign(ot,num_elements(oR));

        current_statement_begin__ = 263;
        int nht(0);
        (void) nht;  // dummy to suppress unused var warning
        stan::math::fill(nht, std::numeric_limits<int>::min());
        stan::math::assign(nht,(ot - ht));

        current_statement_begin__ = 264;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,(ot + uot));

        current_statement_begin__ = 265;
        validate_non_negative_index("R", "ot", ot);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> R(ot);
        stan::math::initialize(R, DUMMY_VAR__);
        stan::math::fill(R, DUMMY_VAR__);
        stan::math::assign(R,oR);

        current_statement_begin__ = 266;
        local_scalar_t__ exp_adj_Rt(DUMMY_VAR__);
        (void) exp_adj_Rt;  // dummy to suppress unused var warning
        stan::math::initialize(exp_adj_Rt, DUMMY_VAR__);
        stan::math::fill(exp_adj_Rt, DUMMY_VAR__);

        current_statement_begin__ = 267;
        validate_non_negative_index("infections", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> infections(t);
        stan::math::initialize(infections, DUMMY_VAR__);
        stan::math::fill(infections, DUMMY_VAR__);
        stan::math::assign(infections,rep_vector(1e-5, t));

        current_statement_begin__ = 268;
        validate_non_negative_index("cum_infections", "ot", ot);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> cum_infections(ot);
        stan::math::initialize(cum_infections, DUMMY_VAR__);
        stan::math::fill(cum_infections, DUMMY_VAR__);
        stan::math::assign(cum_infections,rep_vector(0, ot));

        current_statement_begin__ = 269;
        validate_non_negative_index("infectiousness", "ot", ot);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> infectiousness(ot);
        stan::math::initialize(infectiousness, DUMMY_VAR__);
        stan::math::fill(infectiousness, DUMMY_VAR__);
        stan::math::assign(infectiousness,rep_vector(1e-5, ot));

        current_statement_begin__ = 271;
        validate_non_negative_index("gt_pmf", "max_gt", max_gt);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> gt_pmf(max_gt);
        stan::math::initialize(gt_pmf, DUMMY_VAR__);
        stan::math::fill(gt_pmf, DUMMY_VAR__);
        stan::math::assign(gt_pmf,rep_vector(1e-5, max_gt));

        current_statement_begin__ = 272;
        validate_non_negative_index("gt_indexes", "max_gt", max_gt);
        std::vector<int  > gt_indexes(max_gt, int(0));
        stan::math::fill(gt_indexes, std::numeric_limits<int>::min());


        current_statement_begin__ = 273;
        for (int i = 1; i <= max_gt; ++i) {

            current_statement_begin__ = 274;
            stan::model::assign(gt_indexes, 
                        stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                        ((max_gt - i) + 1), 
                        "assigning variable gt_indexes");
        }
        current_statement_begin__ = 276;
        stan::math::assign(gt_pmf, add(gt_pmf, discretised_gamma_pmf(gt_indexes, get_base1(gt_mean, 1, "gt_mean", 1), get_base1(gt_sd, 1, "gt_sd", 1), max_gt, pstream__)));
        current_statement_begin__ = 278;
        stan::model::assign(infections, 
                    stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()), 
                    stan::math::exp(get_base1(initial_infections, 1, "initial_infections", 1)), 
                    "assigning variable infections");
        current_statement_begin__ = 279;
        if (as_bool(logical_gt(uot, 1))) {

            current_statement_begin__ = 280;
            for (int s = 2; s <= uot; ++s) {

                current_statement_begin__ = 281;
                stan::model::assign(infections, 
                            stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                            stan::math::exp((get_base1(initial_infections, 1, "initial_infections", 1) + (get_base1(initial_growth, 1, "initial_growth", 1) * (s - 1)))), 
                            "assigning variable infections");
            }
        }
        current_statement_begin__ = 285;
        if (as_bool(pop)) {

            current_statement_begin__ = 286;
            stan::model::assign(cum_infections, 
                        stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()), 
                        sum(stan::model::rvalue(infections, stan::model::cons_list(stan::model::index_min_max(1, uot), stan::model::nil_index_list()), "infections")), 
                        "assigning variable cum_infections");
        }
        current_statement_begin__ = 289;
        for (int s = 1; s <= ot; ++s) {

            current_statement_begin__ = 290;
            stan::model::assign(infectiousness, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        (stan::model::rvalue(infectiousness, stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), "infectiousness") + update_infectiousness(infections, gt_pmf, uot, max_gt, s, pstream__)), 
                        "assigning variable infectiousness");
            current_statement_begin__ = 291;
            if (as_bool((primitive_value(pop) && primitive_value(logical_gt(s, nht))))) {

                current_statement_begin__ = 292;
                stan::math::assign(exp_adj_Rt, stan::math::exp(((-(get_base1(R, s, "R", 1)) * get_base1(infectiousness, s, "infectiousness", 1)) / (pop - get_base1(cum_infections, nht, "cum_infections", 1)))));
                current_statement_begin__ = 293;
                stan::math::assign(exp_adj_Rt, (logical_gt(exp_adj_Rt, 1) ? stan::math::promote_scalar<local_scalar_t__>(1) : stan::math::promote_scalar<local_scalar_t__>(exp_adj_Rt) ));
                current_statement_begin__ = 294;
                stan::model::assign(infections, 
                            stan::model::cons_list(stan::model::index_uni((s + uot)), stan::model::nil_index_list()), 
                            ((pop - get_base1(cum_infections, s, "cum_infections", 1)) * (1 - exp_adj_Rt)), 
                            "assigning variable infections");
            } else {

                current_statement_begin__ = 296;
                stan::model::assign(infections, 
                            stan::model::cons_list(stan::model::index_uni((s + uot)), stan::model::nil_index_list()), 
                            (stan::model::rvalue(infections, stan::model::cons_list(stan::model::index_uni((s + uot)), stan::model::nil_index_list()), "infections") + (get_base1(R, s, "R", 1) * get_base1(infectiousness, s, "infectiousness", 1))), 
                            "assigning variable infections");
            }
            current_statement_begin__ = 298;
            if (as_bool((primitive_value(pop) && primitive_value(logical_lt(s, ot))))) {

                current_statement_begin__ = 299;
                stan::model::assign(cum_infections, 
                            stan::model::cons_list(stan::model::index_uni((s + 1)), stan::model::nil_index_list()), 
                            (get_base1(cum_infections, s, "cum_infections", 1) + get_base1(infections, (s + uot), "infections", 1)), 
                            "assigning variable cum_infections");
            }
        }
        current_statement_begin__ = 302;
        return stan::math::promote_scalar<fun_return_scalar_t__>(infections);
        }
    } 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 generate_infections_functor__ {
    template <typename T0__, typename T2__, typename T3__, typename T5__, typename T6__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__, T3__, T5__, typename boost::math::tools::promote_args<T6__>::type>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& oR,
                        const int& uot,
                        const std::vector<T2__>& gt_mean,
                        const std::vector<T3__>& gt_sd,
                        const int& max_gt,
                        const std::vector<T5__>& initial_infections,
                        const std::vector<T6__>& initial_growth,
                        const int& pop,
                        const int& ht, std::ostream* pstream__) const {
        return generate_infections(oR, uot, gt_mean, gt_sd, max_gt, initial_infections, initial_growth, pop, ht, pstream__);
    }
};

template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
deconvolve_infections(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& shifted_cases,
                          const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& noise,
                          const int& fixed,
                          const int& prior, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::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__ = 307;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(shifted_cases));

        current_statement_begin__ = 308;
        validate_non_negative_index("infections", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> infections(t);
        stan::math::initialize(infections, DUMMY_VAR__);
        stan::math::fill(infections, DUMMY_VAR__);
        stan::math::assign(infections,rep_vector(1e-5, t));


        current_statement_begin__ = 309;
        if (as_bool(logical_negation(fixed))) {
            {
            current_statement_begin__ = 310;
            validate_non_negative_index("exp_noise", "t", t);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> exp_noise(t);
            stan::math::initialize(exp_noise, DUMMY_VAR__);
            stan::math::fill(exp_noise, DUMMY_VAR__);
            stan::math::assign(exp_noise,stan::math::exp(noise));


            current_statement_begin__ = 311;
            if (as_bool(logical_eq(prior, 1))) {

                current_statement_begin__ = 312;
                stan::math::assign(infections, add(infections, elt_multiply(shifted_cases, exp_noise)));
            } else if (as_bool(logical_eq(prior, 0))) {

                current_statement_begin__ = 314;
                stan::math::assign(infections, add(infections, exp_noise));
            } else if (as_bool(logical_eq(prior, 2))) {

                current_statement_begin__ = 316;
                stan::model::assign(infections, 
                            stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()), 
                            (get_base1(infections, 1, "infections", 1) + (get_base1(shifted_cases, 1, "shifted_cases", 1) * get_base1(exp_noise, 1, "exp_noise", 1))), 
                            "assigning variable infections");
                current_statement_begin__ = 317;
                for (int i = 2; i <= t; ++i) {

                    current_statement_begin__ = 318;
                    stan::model::assign(infections, 
                                stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                                (get_base1(infections, (i - 1), "infections", 1) * get_base1(exp_noise, i, "exp_noise", 1)), 
                                "assigning variable infections");
                }
            }
            }
        } else {

            current_statement_begin__ = 322;
            stan::math::assign(infections, add(infections, shifted_cases));
        }
        current_statement_begin__ = 324;
        return stan::math::promote_scalar<fun_return_scalar_t__>(infections);
        }
    } 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 deconvolve_infections_functor__ {
    template <typename T0__, typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& shifted_cases,
                          const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& noise,
                          const int& fixed,
                          const int& prior, std::ostream* pstream__) const {
        return deconvolve_infections(shifted_cases, noise, fixed, prior, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T_lp__, typename T_lp_accum__>
void
generation_time_lp(const std::vector<T0__>& gt_mean,
                       const T1__& gt_mean_mean,
                       const T2__& gt_mean_sd,
                       const std::vector<T3__>& gt_sd,
                       const T4__& gt_sd_mean,
                       const T5__& gt_sd_sd,
                       const int& weight, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T_lp__>::type>::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__ = 329;
        lp_accum__.add((normal_log(get_base1(gt_mean, 1, "gt_mean", 1), gt_mean_mean, gt_mean_sd) * weight));
        current_statement_begin__ = 330;
        lp_accum__.add((normal_log(get_base1(gt_sd, 1, "gt_sd", 1), gt_sd_mean, gt_sd_sd) * weight));
    } 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 generation_time_lp_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const std::vector<T0__>& gt_mean,
                       const T1__& gt_mean_mean,
                       const T2__& gt_mean_sd,
                       const std::vector<T3__>& gt_sd,
                       const T4__& gt_sd_mean,
                       const T5__& gt_sd_sd,
                       const int& weight, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return generation_time_lp(gt_mean, gt_mean_mean, gt_mean_sd, gt_sd, gt_sd_mean, gt_sd_sd, weight, lp__, lp_accum__, pstream__);
    }
};

template <typename T0__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__>::type, Eigen::Dynamic, 1>
day_of_week_effect(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
                       const std::vector<int>& day_of_week,
                       const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& effect, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T2__>::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__ = 336;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(reports));

        current_statement_begin__ = 337;
        int wl(0);
        (void) wl;  // dummy to suppress unused var warning
        stan::math::fill(wl, std::numeric_limits<int>::min());
        stan::math::assign(wl,num_elements(effect));

        current_statement_begin__ = 339;
        validate_non_negative_index("scaled_effect", "wl", wl);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> scaled_effect(wl);
        stan::math::initialize(scaled_effect, DUMMY_VAR__);
        stan::math::fill(scaled_effect, DUMMY_VAR__);
        stan::math::assign(scaled_effect,multiply(wl, effect));

        current_statement_begin__ = 340;
        validate_non_negative_index("scaled_reports", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> scaled_reports(t);
        stan::math::initialize(scaled_reports, DUMMY_VAR__);
        stan::math::fill(scaled_reports, DUMMY_VAR__);


        current_statement_begin__ = 341;
        for (int s = 1; s <= t; ++s) {

            current_statement_begin__ = 343;
            stan::model::assign(scaled_reports, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        (get_base1(reports, s, "reports", 1) * get_base1(scaled_effect, get_base1(day_of_week, s, "day_of_week", 1), "scaled_effect", 1)), 
                        "assigning variable scaled_reports");
        }
        current_statement_begin__ = 345;
        return stan::math::promote_scalar<fun_return_scalar_t__>(scaled_reports);
        }
    } 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 day_of_week_effect_functor__ {
    template <typename T0__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
                       const std::vector<int>& day_of_week,
                       const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& effect, std::ostream* pstream__) const {
        return day_of_week_effect(reports, day_of_week, effect, pstream__);
    }
};

template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
scale_obs(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
              const T1__& frac_obs, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::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__ = 350;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(reports));

        current_statement_begin__ = 351;
        validate_non_negative_index("scaled_reports", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> scaled_reports(t);
        stan::math::initialize(scaled_reports, DUMMY_VAR__);
        stan::math::fill(scaled_reports, DUMMY_VAR__);


        current_statement_begin__ = 352;
        stan::math::assign(scaled_reports, multiply(reports, frac_obs));
        current_statement_begin__ = 353;
        return stan::math::promote_scalar<fun_return_scalar_t__>(scaled_reports);
        }
    } 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 scale_obs_functor__ {
    template <typename T0__, typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
              const T1__& frac_obs, std::ostream* pstream__) const {
        return scale_obs(reports, frac_obs, pstream__);
    }
};

template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
truncation_cmf(const T0__& trunc_mean,
                   const T1__& trunc_sd,
                   const int& trunc_max, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::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__ = 357;
        validate_non_negative_index("trunc_indexes", "trunc_max", trunc_max);
        std::vector<int  > trunc_indexes(trunc_max, int(0));
        stan::math::fill(trunc_indexes, std::numeric_limits<int>::min());

        current_statement_begin__ = 358;
        validate_non_negative_index("cmf", "trunc_max", trunc_max);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> cmf(trunc_max);
        stan::math::initialize(cmf, DUMMY_VAR__);
        stan::math::fill(cmf, DUMMY_VAR__);


        current_statement_begin__ = 359;
        for (int i = 1; i <= trunc_max; ++i) {

            current_statement_begin__ = 360;
            stan::model::assign(trunc_indexes, 
                        stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                        (i - 1), 
                        "assigning variable trunc_indexes");
        }
        current_statement_begin__ = 362;
        stan::math::assign(cmf, discretised_lognormal_pmf(trunc_indexes, trunc_mean, trunc_sd, trunc_max, pstream__));
        current_statement_begin__ = 363;
        stan::model::assign(cmf, 
                    stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()), 
                    (get_base1(cmf, 1, "cmf", 1) + 1e-8), 
                    "assigning variable cmf");
        current_statement_begin__ = 364;
        stan::math::assign(cmf, cumulative_sum(cmf));
        current_statement_begin__ = 365;
        stan::math::assign(cmf, reverse_mf(cmf, trunc_max, pstream__));
        current_statement_begin__ = 366;
        return stan::math::promote_scalar<fun_return_scalar_t__>(cmf);
        }
    } 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 truncation_cmf_functor__ {
    template <typename T0__, typename T1__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, 1>
    operator()(const T0__& trunc_mean,
                   const T1__& trunc_sd,
                   const int& trunc_max, std::ostream* pstream__) const {
        return truncation_cmf(trunc_mean, trunc_sd, trunc_max, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__>::type, Eigen::Dynamic, 1>
truncate(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
             const std::vector<T1__>& truncation_mean,
             const std::vector<T2__>& truncation_sd,
             const std::vector<int>& truncation_max,
             const int& reconstruct, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__>::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__ = 371;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(reports));

        current_statement_begin__ = 372;
        int truncation(0);
        (void) truncation;  // dummy to suppress unused var warning
        stan::math::fill(truncation, std::numeric_limits<int>::min());
        stan::math::assign(truncation,num_elements(truncation_mean));

        current_statement_begin__ = 373;
        validate_non_negative_index("trunc_reports", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> trunc_reports(t);
        stan::math::initialize(trunc_reports, DUMMY_VAR__);
        stan::math::fill(trunc_reports, DUMMY_VAR__);
        stan::math::assign(trunc_reports,reports);


        current_statement_begin__ = 374;
        if (as_bool(truncation)) {
            {
            current_statement_begin__ = 376;
            int trunc_max(0);
            (void) trunc_max;  // dummy to suppress unused var warning
            stan::math::fill(trunc_max, std::numeric_limits<int>::min());
            stan::math::assign(trunc_max,(logical_gt(get_base1(truncation_max, 1, "truncation_max", 1), t) ? t : get_base1(truncation_max, 1, "truncation_max", 1) ));

            current_statement_begin__ = 377;
            validate_non_negative_index("trunc_indexes", "trunc_max", trunc_max);
            std::vector<int  > trunc_indexes(trunc_max, int(0));
            stan::math::fill(trunc_indexes, std::numeric_limits<int>::min());

            current_statement_begin__ = 378;
            validate_non_negative_index("cmf", "trunc_max", trunc_max);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> cmf(trunc_max);
            stan::math::initialize(cmf, DUMMY_VAR__);
            stan::math::fill(cmf, DUMMY_VAR__);

            current_statement_begin__ = 379;
            int first_t(0);
            (void) first_t;  // dummy to suppress unused var warning
            stan::math::fill(first_t, std::numeric_limits<int>::min());
            stan::math::assign(first_t,((t - trunc_max) + 1));


            current_statement_begin__ = 380;
            stan::math::assign(cmf, truncation_cmf(get_base1(truncation_mean, 1, "truncation_mean", 1), get_base1(truncation_sd, 1, "truncation_sd", 1), trunc_max, pstream__));
            current_statement_begin__ = 382;
            if (as_bool(reconstruct)) {

                current_statement_begin__ = 383;
                stan::model::assign(trunc_reports, 
                            stan::model::cons_list(stan::model::index_min_max(first_t, t), stan::model::nil_index_list()), 
                            stan::model::deep_copy(elt_divide(stan::model::rvalue(trunc_reports, stan::model::cons_list(stan::model::index_min_max(first_t, t), stan::model::nil_index_list()), "trunc_reports"), cmf)), 
                            "assigning variable trunc_reports");
            } else {

                current_statement_begin__ = 385;
                stan::model::assign(trunc_reports, 
                            stan::model::cons_list(stan::model::index_min_max(first_t, t), stan::model::nil_index_list()), 
                            stan::model::deep_copy(elt_multiply(stan::model::rvalue(trunc_reports, stan::model::cons_list(stan::model::index_min_max(first_t, t), stan::model::nil_index_list()), "trunc_reports"), cmf)), 
                            "assigning variable trunc_reports");
            }
            }
        }
        current_statement_begin__ = 388;
        return stan::math::promote_scalar<fun_return_scalar_t__>(trunc_reports);
        }
    } 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 truncate_functor__ {
    template <typename T0__, typename T1__, typename T2__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__, T2__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
             const std::vector<T1__>& truncation_mean,
             const std::vector<T2__>& truncation_sd,
             const std::vector<int>& truncation_max,
             const int& reconstruct, std::ostream* pstream__) const {
        return truncate(reports, truncation_mean, truncation_sd, truncation_max, reconstruct, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T_lp__, typename T_lp_accum__>
void
truncation_lp(const std::vector<T0__>& truncation_mean,
                  const std::vector<T1__>& truncation_sd,
                  const std::vector<T2__>& trunc_mean_mean,
                  const std::vector<T3__>& trunc_mean_sd,
                  const std::vector<T4__>& trunc_sd_mean,
                  const std::vector<T5__>& trunc_sd_sd, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T_lp__>::type>::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__ = 394;
        int truncation(0);
        (void) truncation;  // dummy to suppress unused var warning
        stan::math::fill(truncation, std::numeric_limits<int>::min());
        stan::math::assign(truncation,num_elements(truncation_mean));


        current_statement_begin__ = 395;
        if (as_bool(truncation)) {

            current_statement_begin__ = 396;
            lp_accum__.add(normal_log<propto__>(truncation_mean, trunc_mean_mean, trunc_mean_sd));
            current_statement_begin__ = 397;
            lp_accum__.add(normal_log<propto__>(truncation_sd, trunc_sd_mean, trunc_sd_sd));
        }
        }
    } 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 truncation_lp_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const std::vector<T0__>& truncation_mean,
                  const std::vector<T1__>& truncation_sd,
                  const std::vector<T2__>& trunc_mean_mean,
                  const std::vector<T3__>& trunc_mean_sd,
                  const std::vector<T4__>& trunc_sd_mean,
                  const std::vector<T5__>& trunc_sd_sd, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return truncation_lp(truncation_mean, truncation_sd, trunc_mean_mean, trunc_mean_sd, trunc_sd_mean, trunc_sd_sd, lp__, lp_accum__, pstream__);
    }
};

template <typename T1__, typename T2__, typename T3__, typename T4__, typename T6__, typename T_lp__, typename T_lp_accum__>
void
report_lp(const std::vector<int>& cases,
              const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& reports,
              const std::vector<T2__>& rep_phi,
              const T3__& phi_mean,
              const T4__& phi_sd,
              const int& model_type,
              const T6__& weight, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__, T2__, T3__, T4__, typename boost::math::tools::promote_args<T6__, T_lp__>::type>::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__ = 404;
        local_scalar_t__ sqrt_phi(DUMMY_VAR__);
        (void) sqrt_phi;  // dummy to suppress unused var warning
        stan::math::initialize(sqrt_phi, DUMMY_VAR__);
        stan::math::fill(sqrt_phi, DUMMY_VAR__);
        stan::math::assign(sqrt_phi,1e5);


        current_statement_begin__ = 405;
        if (as_bool(model_type)) {

            current_statement_begin__ = 407;
            lp_accum__.add(normal_log<propto__>(get_base1(rep_phi, model_type, "rep_phi", 1), phi_mean, phi_sd));
            if (get_base1(rep_phi, model_type, "rep_phi", 1) < 0) lp_accum__.add(-std::numeric_limits<double>::infinity());
            else lp_accum__.add(-normal_ccdf_log(0, phi_mean, phi_sd));
            current_statement_begin__ = 408;
            stan::math::assign(sqrt_phi, (1 / stan::math::sqrt(get_base1(rep_phi, model_type, "rep_phi", 1))));
        }
        current_statement_begin__ = 412;
        if (as_bool(logical_gt(sqrt_phi, 1e4))) {

            current_statement_begin__ = 413;
            if (as_bool(logical_eq(weight, 1))) {

                current_statement_begin__ = 414;
                lp_accum__.add(poisson_log<propto__>(cases, reports));
            } else {

                current_statement_begin__ = 416;
                lp_accum__.add((poisson_log(cases, reports) * weight));
            }
        } else {

            current_statement_begin__ = 419;
            if (as_bool(logical_eq(weight, 1))) {

                current_statement_begin__ = 420;
                lp_accum__.add(neg_binomial_2_log<propto__>(cases, reports, sqrt_phi));
            } else {

                current_statement_begin__ = 422;
                lp_accum__.add(neg_binomial_2_log(cases, reports, sqrt_phi));
            }
        }
        }
    } 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 report_lp_functor__ {
    template <typename T1__, typename T2__, typename T3__, typename T4__, typename T6__, typename T_lp__, typename T_lp_accum__>
        void
    operator()(const std::vector<int>& cases,
              const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& reports,
              const std::vector<T2__>& rep_phi,
              const T3__& phi_mean,
              const T4__& phi_sd,
              const int& model_type,
              const T6__& weight, T_lp__& lp__, T_lp_accum__& lp_accum__, std::ostream* pstream__) const {
        return report_lp(cases, reports, rep_phi, phi_mean, phi_sd, model_type, weight, lp__, lp_accum__, pstream__);
    }
};

template <typename T1__, typename T2__, typename T4__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__, T4__>::type, Eigen::Dynamic, 1>
report_log_lik(const std::vector<int>& cases,
                   const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& reports,
                   const std::vector<T2__>& rep_phi,
                   const int& model_type,
                   const T4__& weight, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__, T2__, T4__>::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__ = 430;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(reports));

        current_statement_begin__ = 431;
        validate_non_negative_index("log_lik", "t", t);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> log_lik(t);
        stan::math::initialize(log_lik, DUMMY_VAR__);
        stan::math::fill(log_lik, DUMMY_VAR__);

        current_statement_begin__ = 432;
        local_scalar_t__ sqrt_phi(DUMMY_VAR__);
        (void) sqrt_phi;  // dummy to suppress unused var warning
        stan::math::initialize(sqrt_phi, DUMMY_VAR__);
        stan::math::fill(sqrt_phi, DUMMY_VAR__);
        stan::math::assign(sqrt_phi,1e5);


        current_statement_begin__ = 433;
        if (as_bool(model_type)) {

            current_statement_begin__ = 435;
            stan::math::assign(sqrt_phi, (1 / stan::math::sqrt(get_base1(rep_phi, model_type, "rep_phi", 1))));
        }
        current_statement_begin__ = 439;
        if (as_bool(logical_gt(sqrt_phi, 1e4))) {

            current_statement_begin__ = 440;
            for (int i = 1; i <= t; ++i) {

                current_statement_begin__ = 441;
                stan::model::assign(log_lik, 
                            stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                            (poisson_log(get_base1(cases, i, "cases", 1), get_base1(reports, i, "reports", 1)) * weight), 
                            "assigning variable log_lik");
            }
        } else {

            current_statement_begin__ = 444;
            for (int i = 1; i <= t; ++i) {

                current_statement_begin__ = 445;
                stan::model::assign(log_lik, 
                            stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                            (neg_binomial_2_log(get_base1(cases, i, "cases", 1), get_base1(reports, i, "reports", 1), sqrt_phi) * weight), 
                            "assigning variable log_lik");
            }
        }
        current_statement_begin__ = 448;
        return stan::math::promote_scalar<fun_return_scalar_t__>(log_lik);
        }
    } 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 report_log_lik_functor__ {
    template <typename T1__, typename T2__, typename T4__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__, T4__>::type, Eigen::Dynamic, 1>
    operator()(const std::vector<int>& cases,
                   const Eigen::Matrix<T1__, Eigen::Dynamic, 1>& reports,
                   const std::vector<T2__>& rep_phi,
                   const int& model_type,
                   const T4__& weight, std::ostream* pstream__) const {
        return report_log_lik(cases, reports, rep_phi, model_type, weight, pstream__);
    }
};

template <typename T0__, typename T1__, class RNG>
std::vector<int>
report_rng(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
               const std::vector<T1__>& rep_phi,
               const int& model_type, RNG& base_rng__, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::type 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__ = 452;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(reports));

        current_statement_begin__ = 453;
        validate_non_negative_index("sampled_reports", "t", t);
        std::vector<int  > sampled_reports(t, int(0));
        stan::math::fill(sampled_reports, std::numeric_limits<int>::min());

        current_statement_begin__ = 454;
        local_scalar_t__ sqrt_phi(DUMMY_VAR__);
        (void) sqrt_phi;  // dummy to suppress unused var warning
        stan::math::initialize(sqrt_phi, DUMMY_VAR__);
        stan::math::fill(sqrt_phi, DUMMY_VAR__);
        stan::math::assign(sqrt_phi,1e5);


        current_statement_begin__ = 455;
        if (as_bool(model_type)) {

            current_statement_begin__ = 456;
            stan::math::assign(sqrt_phi, (1 / stan::math::sqrt(get_base1(rep_phi, model_type, "rep_phi", 1))));
        }
        current_statement_begin__ = 459;
        for (int s = 1; s <= t; ++s) {

            current_statement_begin__ = 461;
            if (as_bool(logical_gt(sqrt_phi, 1e4))) {

                current_statement_begin__ = 462;
                stan::model::assign(sampled_reports, 
                            stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                            poisson_rng((logical_gt(get_base1(reports, s, "reports", 1), 1e8) ? 1e8 : get_base1(reports, s, "reports", 1) ), base_rng__), 
                            "assigning variable sampled_reports");
            } else {

                current_statement_begin__ = 464;
                stan::model::assign(sampled_reports, 
                            stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                            neg_binomial_2_rng((logical_gt(get_base1(reports, s, "reports", 1), 1e8) ? 1e8 : get_base1(reports, s, "reports", 1) ), sqrt_phi, base_rng__), 
                            "assigning variable sampled_reports");
            }
        }
        current_statement_begin__ = 467;
        return stan::math::promote_scalar<fun_return_scalar_t__>(sampled_reports);
        }
    } 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 report_rng_functor__ {
    template <typename T0__, typename T1__, class RNG>
        std::vector<int>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& reports,
               const std::vector<T1__>& rep_phi,
               const int& model_type, RNG& base_rng__, std::ostream* pstream__) const {
        return report_rng(reports, rep_phi, model_type, base_rng__, pstream__);
    }
};

template <typename T0__, typename T2__, typename T3__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__, T3__>::type, Eigen::Dynamic, 1>
calculate_Rt(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& infections,
                 const int& seeding_time,
                 const T2__& gt_mean,
                 const T3__& gt_sd,
                 const int& max_gt,
                 const int& smooth, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T2__, T3__>::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__ = 474;
        validate_non_negative_index("gt_pmf", "max_gt", max_gt);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> gt_pmf(max_gt);
        stan::math::initialize(gt_pmf, DUMMY_VAR__);
        stan::math::fill(gt_pmf, DUMMY_VAR__);

        current_statement_begin__ = 475;
        validate_non_negative_index("gt_indexes", "max_gt", max_gt);
        std::vector<int  > gt_indexes(max_gt, int(0));
        stan::math::fill(gt_indexes, std::numeric_limits<int>::min());

        current_statement_begin__ = 476;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(infections));

        current_statement_begin__ = 477;
        int ot(0);
        (void) ot;  // dummy to suppress unused var warning
        stan::math::fill(ot, std::numeric_limits<int>::min());
        stan::math::assign(ot,(t - seeding_time));

        current_statement_begin__ = 478;
        validate_non_negative_index("R", "ot", ot);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> R(ot);
        stan::math::initialize(R, DUMMY_VAR__);
        stan::math::fill(R, DUMMY_VAR__);

        current_statement_begin__ = 479;
        validate_non_negative_index("sR", "ot", ot);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> sR(ot);
        stan::math::initialize(sR, DUMMY_VAR__);
        stan::math::fill(sR, DUMMY_VAR__);

        current_statement_begin__ = 480;
        validate_non_negative_index("infectiousness", "ot", ot);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> infectiousness(ot);
        stan::math::initialize(infectiousness, DUMMY_VAR__);
        stan::math::fill(infectiousness, DUMMY_VAR__);
        stan::math::assign(infectiousness,rep_vector(1e-5, ot));


        current_statement_begin__ = 482;
        for (int i = 1; i <= max_gt; ++i) {

            current_statement_begin__ = 483;
            stan::model::assign(gt_indexes, 
                        stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                        ((max_gt - i) + 1), 
                        "assigning variable gt_indexes");
        }
        current_statement_begin__ = 485;
        stan::math::assign(gt_pmf, discretised_gamma_pmf(gt_indexes, gt_mean, gt_sd, max_gt, pstream__));
        current_statement_begin__ = 487;
        for (int s = 1; s <= ot; ++s) {

            current_statement_begin__ = 488;
            stan::model::assign(infectiousness, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        (stan::model::rvalue(infectiousness, stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), "infectiousness") + update_infectiousness(infections, gt_pmf, seeding_time, max_gt, s, pstream__)), 
                        "assigning variable infectiousness");
            current_statement_begin__ = 489;
            stan::model::assign(R, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        (get_base1(infections, (s + seeding_time), "infections", 1) / get_base1(infectiousness, s, "infectiousness", 1)), 
                        "assigning variable R");
        }
        current_statement_begin__ = 491;
        if (as_bool(smooth)) {

            current_statement_begin__ = 492;
            for (int s = 1; s <= ot; ++s) {
                {
                current_statement_begin__ = 493;
                local_scalar_t__ window(DUMMY_VAR__);
                (void) window;  // dummy to suppress unused var warning
                stan::math::initialize(window, DUMMY_VAR__);
                stan::math::fill(window, DUMMY_VAR__);
                stan::math::assign(window,0);


                current_statement_begin__ = 494;
                stan::model::assign(sR, 
                            stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                            0, 
                            "assigning variable sR");
                current_statement_begin__ = 495;
                for (int i = std::max(1, (s - smooth)); i <= std::min(ot, (s + smooth)); ++i) {

                    current_statement_begin__ = 496;
                    stan::model::assign(sR, 
                                stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                                (stan::model::rvalue(sR, stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), "sR") + get_base1(R, i, "R", 1)), 
                                "assigning variable sR");
                    current_statement_begin__ = 497;
                    stan::math::assign(window, (window + 1));
                }
                current_statement_begin__ = 499;
                stan::model::assign(sR, 
                            stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                            (get_base1(sR, s, "sR", 1) / window), 
                            "assigning variable sR");
                }
            }
        } else {

            current_statement_begin__ = 502;
            stan::math::assign(sR, R);
        }
        current_statement_begin__ = 504;
        return stan::math::promote_scalar<fun_return_scalar_t__>(sR);
        }
    } 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 calculate_Rt_functor__ {
    template <typename T0__, typename T2__, typename T3__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T2__, T3__>::type, Eigen::Dynamic, 1>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& infections,
                 const int& seeding_time,
                 const T2__& gt_mean,
                 const T3__& gt_sd,
                 const int& max_gt,
                 const int& smooth, std::ostream* pstream__) const {
        return calculate_Rt(infections, seeding_time, gt_mean, gt_sd, max_gt, smooth, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__>
std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__>::type>
R_to_growth(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& R,
                const T1__& gt_mean,
                const T2__& gt_sd, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__>::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__ = 508;
        local_scalar_t__ k(DUMMY_VAR__);
        (void) k;  // dummy to suppress unused var warning
        stan::math::initialize(k, DUMMY_VAR__);
        stan::math::fill(k, DUMMY_VAR__);
        stan::math::assign(k,pow((gt_sd / gt_mean), 2));

        current_statement_begin__ = 509;
        int t(0);
        (void) t;  // dummy to suppress unused var warning
        stan::math::fill(t, std::numeric_limits<int>::min());
        stan::math::assign(t,num_elements(R));

        current_statement_begin__ = 510;
        validate_non_negative_index("r", "t", t);
        std::vector<local_scalar_t__  > r(t, local_scalar_t__(DUMMY_VAR__));
        stan::math::initialize(r, DUMMY_VAR__);
        stan::math::fill(r, DUMMY_VAR__);


        current_statement_begin__ = 511;
        for (int s = 1; s <= t; ++s) {

            current_statement_begin__ = 512;
            stan::model::assign(r, 
                        stan::model::cons_list(stan::model::index_uni(s), stan::model::nil_index_list()), 
                        ((pow(get_base1(R, s, "R", 1), k) - 1) / (k * gt_mean)), 
                        "assigning variable r");
        }
        current_statement_begin__ = 514;
        return stan::math::promote_scalar<fun_return_scalar_t__>(r);
        }
    } 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 R_to_growth_functor__ {
    template <typename T0__, typename T1__, typename T2__>
        std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__>::type>
    operator()(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& R,
                const T1__& gt_mean,
                const T2__& gt_sd, std::ostream* pstream__) const {
        return R_to_growth(R, gt_mean, gt_sd, pstream__);
    }
};

class model_estimate_infections
  : public stan::model::model_base_crtp<model_estimate_infections> {
private:
        int t;
        int seeding_time;
        int horizon;
        int future_time;
        std::vector<int> cases;
        vector_d shifted_cases;
        int delays;
        std::vector<double> delay_mean_sd;
        std::vector<double> delay_mean_mean;
        std::vector<double> delay_sd_mean;
        std::vector<double> delay_sd_sd;
        std::vector<int> max_delay;
        double L;
        int M;
        double ls_meanlog;
        double ls_sdlog;
        double ls_min;
        double ls_max;
        double alpha_sd;
        int gp_type;
        int stationary;
        int fixed;
        double gt_mean_sd;
        double gt_mean_mean;
        double gt_sd_mean;
        double gt_sd_sd;
        int max_gt;
        int estimate_r;
        double prior_infections;
        double prior_growth;
        double r_mean;
        double r_sd;
        int bp_n;
        std::vector<int> breakpoints;
        int future_fixed;
        int fixed_from;
        int pop;
        int backcalc_prior;
        int rt_half_window;
        std::vector<int> day_of_week;
        int model_type;
        double phi_mean;
        double phi_sd;
        int week_effect;
        int truncation;
        std::vector<double> trunc_mean_mean;
        std::vector<double> trunc_mean_sd;
        std::vector<double> trunc_sd_mean;
        std::vector<double> trunc_sd_sd;
        std::vector<int> max_truncation;
        int obs_scale;
        double obs_scale_mean;
        double obs_scale_sd;
        double obs_weight;
        int likelihood;
        int return_likelihood;
        int ot;
        int ot_h;
        int noise_terms;
        matrix_d PHI;
        double r_logmean;
        double r_logsd;
public:
    model_estimate_infections(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : model_base_crtp(0) {
        ctor_body(context__, 0, pstream__);
    }

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

    void ctor_body(stan::io::var_context& context__,
                   unsigned int random_seed__,
                   std::ostream* pstream__) {
        typedef double local_scalar_t__;

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

        current_statement_begin__ = -1;

        static const char* function__ = "model_estimate_infections_namespace::model_estimate_infections";
        (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__ = 521;
            context__.validate_dims("data initialization", "t", "int", context__.to_vec());
            t = int(0);
            vals_i__ = context__.vals_i("t");
            pos__ = 0;
            t = vals_i__[pos__++];

            current_statement_begin__ = 522;
            context__.validate_dims("data initialization", "seeding_time", "int", context__.to_vec());
            seeding_time = int(0);
            vals_i__ = context__.vals_i("seeding_time");
            pos__ = 0;
            seeding_time = vals_i__[pos__++];

            current_statement_begin__ = 523;
            context__.validate_dims("data initialization", "horizon", "int", context__.to_vec());
            horizon = int(0);
            vals_i__ = context__.vals_i("horizon");
            pos__ = 0;
            horizon = vals_i__[pos__++];

            current_statement_begin__ = 524;
            context__.validate_dims("data initialization", "future_time", "int", context__.to_vec());
            future_time = int(0);
            vals_i__ = context__.vals_i("future_time");
            pos__ = 0;
            future_time = vals_i__[pos__++];

            current_statement_begin__ = 525;
            validate_non_negative_index("cases", "((t - horizon) - seeding_time)", ((t - horizon) - seeding_time));
            context__.validate_dims("data initialization", "cases", "int", context__.to_vec(((t - horizon) - seeding_time)));
            cases = std::vector<int>(((t - horizon) - seeding_time), int(0));
            vals_i__ = context__.vals_i("cases");
            pos__ = 0;
            size_t cases_k_0_max__ = ((t - horizon) - seeding_time);
            for (size_t k_0__ = 0; k_0__ < cases_k_0_max__; ++k_0__) {
                cases[k_0__] = vals_i__[pos__++];
            }
            size_t cases_i_0_max__ = ((t - horizon) - seeding_time);
            for (size_t i_0__ = 0; i_0__ < cases_i_0_max__; ++i_0__) {
                check_greater_or_equal(function__, "cases[i_0__]", cases[i_0__], 0);
            }

            current_statement_begin__ = 526;
            validate_non_negative_index("shifted_cases", "t", t);
            context__.validate_dims("data initialization", "shifted_cases", "vector_d", context__.to_vec(t));
            shifted_cases = Eigen::Matrix<double, Eigen::Dynamic, 1>(t);
            vals_r__ = context__.vals_r("shifted_cases");
            pos__ = 0;
            size_t shifted_cases_j_1_max__ = t;
            for (size_t j_1__ = 0; j_1__ < shifted_cases_j_1_max__; ++j_1__) {
                shifted_cases(j_1__) = vals_r__[pos__++];
            }
            check_greater_or_equal(function__, "shifted_cases", shifted_cases, 0);

            current_statement_begin__ = 528;
            context__.validate_dims("data initialization", "delays", "int", context__.to_vec());
            delays = int(0);
            vals_i__ = context__.vals_i("delays");
            pos__ = 0;
            delays = vals_i__[pos__++];

            current_statement_begin__ = 529;
            validate_non_negative_index("delay_mean_sd", "delays", delays);
            context__.validate_dims("data initialization", "delay_mean_sd", "double", context__.to_vec(delays));
            delay_mean_sd = std::vector<double>(delays, double(0));
            vals_r__ = context__.vals_r("delay_mean_sd");
            pos__ = 0;
            size_t delay_mean_sd_k_0_max__ = delays;
            for (size_t k_0__ = 0; k_0__ < delay_mean_sd_k_0_max__; ++k_0__) {
                delay_mean_sd[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 530;
            validate_non_negative_index("delay_mean_mean", "delays", delays);
            context__.validate_dims("data initialization", "delay_mean_mean", "double", context__.to_vec(delays));
            delay_mean_mean = std::vector<double>(delays, double(0));
            vals_r__ = context__.vals_r("delay_mean_mean");
            pos__ = 0;
            size_t delay_mean_mean_k_0_max__ = delays;
            for (size_t k_0__ = 0; k_0__ < delay_mean_mean_k_0_max__; ++k_0__) {
                delay_mean_mean[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 531;
            validate_non_negative_index("delay_sd_mean", "delays", delays);
            context__.validate_dims("data initialization", "delay_sd_mean", "double", context__.to_vec(delays));
            delay_sd_mean = std::vector<double>(delays, double(0));
            vals_r__ = context__.vals_r("delay_sd_mean");
            pos__ = 0;
            size_t delay_sd_mean_k_0_max__ = delays;
            for (size_t k_0__ = 0; k_0__ < delay_sd_mean_k_0_max__; ++k_0__) {
                delay_sd_mean[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 532;
            validate_non_negative_index("delay_sd_sd", "delays", delays);
            context__.validate_dims("data initialization", "delay_sd_sd", "double", context__.to_vec(delays));
            delay_sd_sd = std::vector<double>(delays, double(0));
            vals_r__ = context__.vals_r("delay_sd_sd");
            pos__ = 0;
            size_t delay_sd_sd_k_0_max__ = delays;
            for (size_t k_0__ = 0; k_0__ < delay_sd_sd_k_0_max__; ++k_0__) {
                delay_sd_sd[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 533;
            validate_non_negative_index("max_delay", "delays", delays);
            context__.validate_dims("data initialization", "max_delay", "int", context__.to_vec(delays));
            max_delay = std::vector<int>(delays, int(0));
            vals_i__ = context__.vals_i("max_delay");
            pos__ = 0;
            size_t max_delay_k_0_max__ = delays;
            for (size_t k_0__ = 0; k_0__ < max_delay_k_0_max__; ++k_0__) {
                max_delay[k_0__] = vals_i__[pos__++];
            }

            current_statement_begin__ = 535;
            context__.validate_dims("data initialization", "L", "double", context__.to_vec());
            L = double(0);
            vals_r__ = context__.vals_r("L");
            pos__ = 0;
            L = vals_r__[pos__++];

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

            current_statement_begin__ = 537;
            context__.validate_dims("data initialization", "ls_meanlog", "double", context__.to_vec());
            ls_meanlog = double(0);
            vals_r__ = context__.vals_r("ls_meanlog");
            pos__ = 0;
            ls_meanlog = vals_r__[pos__++];

            current_statement_begin__ = 538;
            context__.validate_dims("data initialization", "ls_sdlog", "double", context__.to_vec());
            ls_sdlog = double(0);
            vals_r__ = context__.vals_r("ls_sdlog");
            pos__ = 0;
            ls_sdlog = vals_r__[pos__++];

            current_statement_begin__ = 539;
            context__.validate_dims("data initialization", "ls_min", "double", context__.to_vec());
            ls_min = double(0);
            vals_r__ = context__.vals_r("ls_min");
            pos__ = 0;
            ls_min = vals_r__[pos__++];
            check_greater_or_equal(function__, "ls_min", ls_min, 0);

            current_statement_begin__ = 540;
            context__.validate_dims("data initialization", "ls_max", "double", context__.to_vec());
            ls_max = double(0);
            vals_r__ = context__.vals_r("ls_max");
            pos__ = 0;
            ls_max = vals_r__[pos__++];
            check_greater_or_equal(function__, "ls_max", ls_max, 0);

            current_statement_begin__ = 541;
            context__.validate_dims("data initialization", "alpha_sd", "double", context__.to_vec());
            alpha_sd = double(0);
            vals_r__ = context__.vals_r("alpha_sd");
            pos__ = 0;
            alpha_sd = vals_r__[pos__++];

            current_statement_begin__ = 542;
            context__.validate_dims("data initialization", "gp_type", "int", context__.to_vec());
            gp_type = int(0);
            vals_i__ = context__.vals_i("gp_type");
            pos__ = 0;
            gp_type = vals_i__[pos__++];

            current_statement_begin__ = 543;
            context__.validate_dims("data initialization", "stationary", "int", context__.to_vec());
            stationary = int(0);
            vals_i__ = context__.vals_i("stationary");
            pos__ = 0;
            stationary = vals_i__[pos__++];

            current_statement_begin__ = 544;
            context__.validate_dims("data initialization", "fixed", "int", context__.to_vec());
            fixed = int(0);
            vals_i__ = context__.vals_i("fixed");
            pos__ = 0;
            fixed = vals_i__[pos__++];

            current_statement_begin__ = 546;
            context__.validate_dims("data initialization", "gt_mean_sd", "double", context__.to_vec());
            gt_mean_sd = double(0);
            vals_r__ = context__.vals_r("gt_mean_sd");
            pos__ = 0;
            gt_mean_sd = vals_r__[pos__++];

            current_statement_begin__ = 547;
            context__.validate_dims("data initialization", "gt_mean_mean", "double", context__.to_vec());
            gt_mean_mean = double(0);
            vals_r__ = context__.vals_r("gt_mean_mean");
            pos__ = 0;
            gt_mean_mean = vals_r__[pos__++];

            current_statement_begin__ = 548;
            context__.validate_dims("data initialization", "gt_sd_mean", "double", context__.to_vec());
            gt_sd_mean = double(0);
            vals_r__ = context__.vals_r("gt_sd_mean");
            pos__ = 0;
            gt_sd_mean = vals_r__[pos__++];

            current_statement_begin__ = 549;
            context__.validate_dims("data initialization", "gt_sd_sd", "double", context__.to_vec());
            gt_sd_sd = double(0);
            vals_r__ = context__.vals_r("gt_sd_sd");
            pos__ = 0;
            gt_sd_sd = vals_r__[pos__++];

            current_statement_begin__ = 550;
            context__.validate_dims("data initialization", "max_gt", "int", context__.to_vec());
            max_gt = int(0);
            vals_i__ = context__.vals_i("max_gt");
            pos__ = 0;
            max_gt = vals_i__[pos__++];

            current_statement_begin__ = 552;
            context__.validate_dims("data initialization", "estimate_r", "int", context__.to_vec());
            estimate_r = int(0);
            vals_i__ = context__.vals_i("estimate_r");
            pos__ = 0;
            estimate_r = vals_i__[pos__++];

            current_statement_begin__ = 553;
            context__.validate_dims("data initialization", "prior_infections", "double", context__.to_vec());
            prior_infections = double(0);
            vals_r__ = context__.vals_r("prior_infections");
            pos__ = 0;
            prior_infections = vals_r__[pos__++];

            current_statement_begin__ = 554;
            context__.validate_dims("data initialization", "prior_growth", "double", context__.to_vec());
            prior_growth = double(0);
            vals_r__ = context__.vals_r("prior_growth");
            pos__ = 0;
            prior_growth = vals_r__[pos__++];

            current_statement_begin__ = 555;
            context__.validate_dims("data initialization", "r_mean", "double", context__.to_vec());
            r_mean = double(0);
            vals_r__ = context__.vals_r("r_mean");
            pos__ = 0;
            r_mean = vals_r__[pos__++];
            check_greater_or_equal(function__, "r_mean", r_mean, 0);

            current_statement_begin__ = 556;
            context__.validate_dims("data initialization", "r_sd", "double", context__.to_vec());
            r_sd = double(0);
            vals_r__ = context__.vals_r("r_sd");
            pos__ = 0;
            r_sd = vals_r__[pos__++];
            check_greater_or_equal(function__, "r_sd", r_sd, 0);

            current_statement_begin__ = 557;
            context__.validate_dims("data initialization", "bp_n", "int", context__.to_vec());
            bp_n = int(0);
            vals_i__ = context__.vals_i("bp_n");
            pos__ = 0;
            bp_n = vals_i__[pos__++];

            current_statement_begin__ = 558;
            validate_non_negative_index("breakpoints", "(t - seeding_time)", (t - seeding_time));
            context__.validate_dims("data initialization", "breakpoints", "int", context__.to_vec((t - seeding_time)));
            breakpoints = std::vector<int>((t - seeding_time), int(0));
            vals_i__ = context__.vals_i("breakpoints");
            pos__ = 0;
            size_t breakpoints_k_0_max__ = (t - seeding_time);
            for (size_t k_0__ = 0; k_0__ < breakpoints_k_0_max__; ++k_0__) {
                breakpoints[k_0__] = vals_i__[pos__++];
            }

            current_statement_begin__ = 559;
            context__.validate_dims("data initialization", "future_fixed", "int", context__.to_vec());
            future_fixed = int(0);
            vals_i__ = context__.vals_i("future_fixed");
            pos__ = 0;
            future_fixed = vals_i__[pos__++];

            current_statement_begin__ = 560;
            context__.validate_dims("data initialization", "fixed_from", "int", context__.to_vec());
            fixed_from = int(0);
            vals_i__ = context__.vals_i("fixed_from");
            pos__ = 0;
            fixed_from = vals_i__[pos__++];

            current_statement_begin__ = 561;
            context__.validate_dims("data initialization", "pop", "int", context__.to_vec());
            pop = int(0);
            vals_i__ = context__.vals_i("pop");
            pos__ = 0;
            pop = vals_i__[pos__++];

            current_statement_begin__ = 563;
            context__.validate_dims("data initialization", "backcalc_prior", "int", context__.to_vec());
            backcalc_prior = int(0);
            vals_i__ = context__.vals_i("backcalc_prior");
            pos__ = 0;
            backcalc_prior = vals_i__[pos__++];

            current_statement_begin__ = 564;
            context__.validate_dims("data initialization", "rt_half_window", "int", context__.to_vec());
            rt_half_window = int(0);
            vals_i__ = context__.vals_i("rt_half_window");
            pos__ = 0;
            rt_half_window = vals_i__[pos__++];

            current_statement_begin__ = 566;
            validate_non_negative_index("day_of_week", "(t - seeding_time)", (t - seeding_time));
            context__.validate_dims("data initialization", "day_of_week", "int", context__.to_vec((t - seeding_time)));
            day_of_week = std::vector<int>((t - seeding_time), int(0));
            vals_i__ = context__.vals_i("day_of_week");
            pos__ = 0;
            size_t day_of_week_k_0_max__ = (t - seeding_time);
            for (size_t k_0__ = 0; k_0__ < day_of_week_k_0_max__; ++k_0__) {
                day_of_week[k_0__] = vals_i__[pos__++];
            }

            current_statement_begin__ = 567;
            context__.validate_dims("data initialization", "model_type", "int", context__.to_vec());
            model_type = int(0);
            vals_i__ = context__.vals_i("model_type");
            pos__ = 0;
            model_type = vals_i__[pos__++];

            current_statement_begin__ = 568;
            context__.validate_dims("data initialization", "phi_mean", "double", context__.to_vec());
            phi_mean = double(0);
            vals_r__ = context__.vals_r("phi_mean");
            pos__ = 0;
            phi_mean = vals_r__[pos__++];

            current_statement_begin__ = 569;
            context__.validate_dims("data initialization", "phi_sd", "double", context__.to_vec());
            phi_sd = double(0);
            vals_r__ = context__.vals_r("phi_sd");
            pos__ = 0;
            phi_sd = vals_r__[pos__++];

            current_statement_begin__ = 570;
            context__.validate_dims("data initialization", "week_effect", "int", context__.to_vec());
            week_effect = int(0);
            vals_i__ = context__.vals_i("week_effect");
            pos__ = 0;
            week_effect = vals_i__[pos__++];

            current_statement_begin__ = 571;
            context__.validate_dims("data initialization", "truncation", "int", context__.to_vec());
            truncation = int(0);
            vals_i__ = context__.vals_i("truncation");
            pos__ = 0;
            truncation = vals_i__[pos__++];

            current_statement_begin__ = 572;
            validate_non_negative_index("trunc_mean_mean", "truncation", truncation);
            context__.validate_dims("data initialization", "trunc_mean_mean", "double", context__.to_vec(truncation));
            trunc_mean_mean = std::vector<double>(truncation, double(0));
            vals_r__ = context__.vals_r("trunc_mean_mean");
            pos__ = 0;
            size_t trunc_mean_mean_k_0_max__ = truncation;
            for (size_t k_0__ = 0; k_0__ < trunc_mean_mean_k_0_max__; ++k_0__) {
                trunc_mean_mean[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 573;
            validate_non_negative_index("trunc_mean_sd", "truncation", truncation);
            context__.validate_dims("data initialization", "trunc_mean_sd", "double", context__.to_vec(truncation));
            trunc_mean_sd = std::vector<double>(truncation, double(0));
            vals_r__ = context__.vals_r("trunc_mean_sd");
            pos__ = 0;
            size_t trunc_mean_sd_k_0_max__ = truncation;
            for (size_t k_0__ = 0; k_0__ < trunc_mean_sd_k_0_max__; ++k_0__) {
                trunc_mean_sd[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 574;
            validate_non_negative_index("trunc_sd_mean", "truncation", truncation);
            context__.validate_dims("data initialization", "trunc_sd_mean", "double", context__.to_vec(truncation));
            trunc_sd_mean = std::vector<double>(truncation, double(0));
            vals_r__ = context__.vals_r("trunc_sd_mean");
            pos__ = 0;
            size_t trunc_sd_mean_k_0_max__ = truncation;
            for (size_t k_0__ = 0; k_0__ < trunc_sd_mean_k_0_max__; ++k_0__) {
                trunc_sd_mean[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 575;
            validate_non_negative_index("trunc_sd_sd", "truncation", truncation);
            context__.validate_dims("data initialization", "trunc_sd_sd", "double", context__.to_vec(truncation));
            trunc_sd_sd = std::vector<double>(truncation, double(0));
            vals_r__ = context__.vals_r("trunc_sd_sd");
            pos__ = 0;
            size_t trunc_sd_sd_k_0_max__ = truncation;
            for (size_t k_0__ = 0; k_0__ < trunc_sd_sd_k_0_max__; ++k_0__) {
                trunc_sd_sd[k_0__] = vals_r__[pos__++];
            }

            current_statement_begin__ = 576;
            validate_non_negative_index("max_truncation", "truncation", truncation);
            context__.validate_dims("data initialization", "max_truncation", "int", context__.to_vec(truncation));
            max_truncation = std::vector<int>(truncation, int(0));
            vals_i__ = context__.vals_i("max_truncation");
            pos__ = 0;
            size_t max_truncation_k_0_max__ = truncation;
            for (size_t k_0__ = 0; k_0__ < max_truncation_k_0_max__; ++k_0__) {
                max_truncation[k_0__] = vals_i__[pos__++];
            }

            current_statement_begin__ = 577;
            context__.validate_dims("data initialization", "obs_scale", "int", context__.to_vec());
            obs_scale = int(0);
            vals_i__ = context__.vals_i("obs_scale");
            pos__ = 0;
            obs_scale = vals_i__[pos__++];

            current_statement_begin__ = 578;
            context__.validate_dims("data initialization", "obs_scale_mean", "double", context__.to_vec());
            obs_scale_mean = double(0);
            vals_r__ = context__.vals_r("obs_scale_mean");
            pos__ = 0;
            obs_scale_mean = vals_r__[pos__++];

            current_statement_begin__ = 579;
            context__.validate_dims("data initialization", "obs_scale_sd", "double", context__.to_vec());
            obs_scale_sd = double(0);
            vals_r__ = context__.vals_r("obs_scale_sd");
            pos__ = 0;
            obs_scale_sd = vals_r__[pos__++];

            current_statement_begin__ = 580;
            context__.validate_dims("data initialization", "obs_weight", "double", context__.to_vec());
            obs_weight = double(0);
            vals_r__ = context__.vals_r("obs_weight");
            pos__ = 0;
            obs_weight = vals_r__[pos__++];

            current_statement_begin__ = 581;
            context__.validate_dims("data initialization", "likelihood", "int", context__.to_vec());
            likelihood = int(0);
            vals_i__ = context__.vals_i("likelihood");
            pos__ = 0;
            likelihood = vals_i__[pos__++];

            current_statement_begin__ = 582;
            context__.validate_dims("data initialization", "return_likelihood", "int", context__.to_vec());
            return_likelihood = int(0);
            vals_i__ = context__.vals_i("return_likelihood");
            pos__ = 0;
            return_likelihood = vals_i__[pos__++];


            // initialize transformed data variables
            current_statement_begin__ = 587;
            ot = int(0);
            stan::math::fill(ot, std::numeric_limits<int>::min());
            stan::math::assign(ot,((t - seeding_time) - horizon));

            current_statement_begin__ = 588;
            ot_h = int(0);
            stan::math::fill(ot_h, std::numeric_limits<int>::min());
            stan::math::assign(ot_h,(ot + horizon));

            current_statement_begin__ = 590;
            noise_terms = int(0);
            stan::math::fill(noise_terms, std::numeric_limits<int>::min());
            stan::math::assign(noise_terms,setup_noise(ot_h, t, horizon, estimate_r, stationary, future_fixed, fixed_from, pstream__));

            current_statement_begin__ = 591;
            validate_non_negative_index("PHI", "noise_terms", noise_terms);
            validate_non_negative_index("PHI", "M", M);
            PHI = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(noise_terms, M);
            stan::math::fill(PHI, DUMMY_VAR__);
            stan::math::assign(PHI,setup_gp(M, L, noise_terms, pstream__));

            current_statement_begin__ = 593;
            r_logmean = double(0);
            stan::math::fill(r_logmean, DUMMY_VAR__);
            stan::math::assign(r_logmean,stan::math::log((pow(r_mean, 2) / stan::math::sqrt((pow(r_sd, 2) + pow(r_mean, 2))))));

            current_statement_begin__ = 594;
            r_logsd = double(0);
            stan::math::fill(r_logsd, DUMMY_VAR__);
            stan::math::assign(r_logsd,stan::math::sqrt(stan::math::log((1 + (pow(r_sd, 2) / pow(r_mean, 2))))));

            // execute transformed data statements

            // validate transformed data

            // validate, set parameter ranges
            num_params_r__ = 0U;
            param_ranges_i__.clear();
            current_statement_begin__ = 599;
            validate_non_negative_index("rho", "(fixed ? 0 : 1 )", (fixed ? 0 : 1 ));
            num_params_r__ += (1 * (fixed ? 0 : 1 ));
            current_statement_begin__ = 600;
            validate_non_negative_index("alpha", "(fixed ? 0 : 1 )", (fixed ? 0 : 1 ));
            num_params_r__ += (1 * (fixed ? 0 : 1 ));
            current_statement_begin__ = 601;
            validate_non_negative_index("eta", "(fixed ? 0 : M )", (fixed ? 0 : M ));
            num_params_r__ += (fixed ? 0 : M );
            current_statement_begin__ = 603;
            validate_non_negative_index("log_R", "estimate_r", estimate_r);
            num_params_r__ += estimate_r;
            current_statement_begin__ = 604;
            validate_non_negative_index("initial_infections", "estimate_r", estimate_r);
            num_params_r__ += (1 * estimate_r);
            current_statement_begin__ = 605;
            validate_non_negative_index("initial_growth", "((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 )", ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 ));
            num_params_r__ += (1 * ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 ));
            current_statement_begin__ = 606;
            validate_non_negative_index("gt_mean", "estimate_r", estimate_r);
            num_params_r__ += (1 * estimate_r);
            current_statement_begin__ = 607;
            validate_non_negative_index("gt_sd", "estimate_r", estimate_r);
            num_params_r__ += (1 * estimate_r);
            current_statement_begin__ = 608;
            validate_non_negative_index("bp_sd", "(logical_gt(bp_n, 0) ? 1 : 0 )", (logical_gt(bp_n, 0) ? 1 : 0 ));
            num_params_r__ += (1 * (logical_gt(bp_n, 0) ? 1 : 0 ));
            current_statement_begin__ = 609;
            validate_non_negative_index("bp_effects", "bp_n", bp_n);
            num_params_r__ += (1 * bp_n);
            current_statement_begin__ = 611;
            validate_non_negative_index("delay_mean", "delays", delays);
            num_params_r__ += (1 * delays);
            current_statement_begin__ = 612;
            validate_non_negative_index("delay_sd", "delays", delays);
            num_params_r__ += (1 * delays);
            current_statement_begin__ = 613;
            validate_non_negative_index("day_of_week_simplex", "week_effect", week_effect);
            num_params_r__ += (week_effect - 1);
            current_statement_begin__ = 614;
            validate_non_negative_index("frac_obs", "obs_scale", obs_scale);
            num_params_r__ += (1 * obs_scale);
            current_statement_begin__ = 615;
            validate_non_negative_index("truncation_mean", "truncation", truncation);
            num_params_r__ += (1 * truncation);
            current_statement_begin__ = 616;
            validate_non_negative_index("truncation_sd", "truncation", truncation);
            num_params_r__ += (1 * truncation);
            current_statement_begin__ = 617;
            validate_non_negative_index("rep_phi", "model_type", model_type);
            num_params_r__ += (1 * model_type);
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(e, current_statement_begin__, prog_reader__());
            // Next line prevents compiler griping about no return
            throw std::runtime_error("*** IF YOU SEE THIS, PLEASE REPORT A BUG ***");
        }
    }

    ~model_estimate_infections() { }


    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__ = 599;
        if (!(context__.contains_r("rho")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable rho missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("rho");
        pos__ = 0U;
        validate_non_negative_index("rho", "(fixed ? 0 : 1 )", (fixed ? 0 : 1 ));
        context__.validate_dims("parameter initialization", "rho", "double", context__.to_vec((fixed ? 0 : 1 )));
        std::vector<double> rho((fixed ? 0 : 1 ), double(0));
        size_t rho_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < rho_k_0_max__; ++k_0__) {
            rho[k_0__] = vals_r__[pos__++];
        }
        size_t rho_i_0_max__ = (fixed ? 0 : 1 );
        for (size_t i_0__ = 0; i_0__ < rho_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lub_unconstrain(ls_min, ls_max, rho[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable rho: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 600;
        if (!(context__.contains_r("alpha")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable alpha missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("alpha");
        pos__ = 0U;
        validate_non_negative_index("alpha", "(fixed ? 0 : 1 )", (fixed ? 0 : 1 ));
        context__.validate_dims("parameter initialization", "alpha", "double", context__.to_vec((fixed ? 0 : 1 )));
        std::vector<double> alpha((fixed ? 0 : 1 ), double(0));
        size_t alpha_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < alpha_k_0_max__; ++k_0__) {
            alpha[k_0__] = vals_r__[pos__++];
        }
        size_t alpha_i_0_max__ = (fixed ? 0 : 1 );
        for (size_t i_0__ = 0; i_0__ < alpha_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, alpha[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable alpha: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 601;
        if (!(context__.contains_r("eta")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable eta missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("eta");
        pos__ = 0U;
        validate_non_negative_index("eta", "(fixed ? 0 : M )", (fixed ? 0 : M ));
        context__.validate_dims("parameter initialization", "eta", "vector_d", context__.to_vec((fixed ? 0 : M )));
        Eigen::Matrix<double, Eigen::Dynamic, 1> eta((fixed ? 0 : M ));
        size_t eta_j_1_max__ = (fixed ? 0 : M );
        for (size_t j_1__ = 0; j_1__ < eta_j_1_max__; ++j_1__) {
            eta(j_1__) = vals_r__[pos__++];
        }
        try {
            writer__.vector_unconstrain(eta);
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable eta: ") + e.what()), current_statement_begin__, prog_reader__());
        }

        current_statement_begin__ = 603;
        if (!(context__.contains_r("log_R")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable log_R missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("log_R");
        pos__ = 0U;
        validate_non_negative_index("log_R", "estimate_r", estimate_r);
        context__.validate_dims("parameter initialization", "log_R", "vector_d", context__.to_vec(estimate_r));
        Eigen::Matrix<double, Eigen::Dynamic, 1> log_R(estimate_r);
        size_t log_R_j_1_max__ = estimate_r;
        for (size_t j_1__ = 0; j_1__ < log_R_j_1_max__; ++j_1__) {
            log_R(j_1__) = vals_r__[pos__++];
        }
        try {
            writer__.vector_unconstrain(log_R);
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable log_R: ") + e.what()), current_statement_begin__, prog_reader__());
        }

        current_statement_begin__ = 604;
        if (!(context__.contains_r("initial_infections")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable initial_infections missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("initial_infections");
        pos__ = 0U;
        validate_non_negative_index("initial_infections", "estimate_r", estimate_r);
        context__.validate_dims("parameter initialization", "initial_infections", "double", context__.to_vec(estimate_r));
        std::vector<double> initial_infections(estimate_r, double(0));
        size_t initial_infections_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < initial_infections_k_0_max__; ++k_0__) {
            initial_infections[k_0__] = vals_r__[pos__++];
        }
        size_t initial_infections_i_0_max__ = estimate_r;
        for (size_t i_0__ = 0; i_0__ < initial_infections_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_unconstrain(initial_infections[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable initial_infections: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 605;
        if (!(context__.contains_r("initial_growth")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable initial_growth missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("initial_growth");
        pos__ = 0U;
        validate_non_negative_index("initial_growth", "((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 )", ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 ));
        context__.validate_dims("parameter initialization", "initial_growth", "double", context__.to_vec(((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 )));
        std::vector<double> initial_growth(((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 ), double(0));
        size_t initial_growth_k_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < initial_growth_k_0_max__; ++k_0__) {
            initial_growth[k_0__] = vals_r__[pos__++];
        }
        size_t initial_growth_i_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
        for (size_t i_0__ = 0; i_0__ < initial_growth_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_unconstrain(initial_growth[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable initial_growth: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 606;
        if (!(context__.contains_r("gt_mean")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable gt_mean missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("gt_mean");
        pos__ = 0U;
        validate_non_negative_index("gt_mean", "estimate_r", estimate_r);
        context__.validate_dims("parameter initialization", "gt_mean", "double", context__.to_vec(estimate_r));
        std::vector<double> gt_mean(estimate_r, double(0));
        size_t gt_mean_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_mean_k_0_max__; ++k_0__) {
            gt_mean[k_0__] = vals_r__[pos__++];
        }
        size_t gt_mean_i_0_max__ = estimate_r;
        for (size_t i_0__ = 0; i_0__ < gt_mean_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lub_unconstrain(0, max_gt, gt_mean[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable gt_mean: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 607;
        if (!(context__.contains_r("gt_sd")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable gt_sd missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("gt_sd");
        pos__ = 0U;
        validate_non_negative_index("gt_sd", "estimate_r", estimate_r);
        context__.validate_dims("parameter initialization", "gt_sd", "double", context__.to_vec(estimate_r));
        std::vector<double> gt_sd(estimate_r, double(0));
        size_t gt_sd_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_sd_k_0_max__; ++k_0__) {
            gt_sd[k_0__] = vals_r__[pos__++];
        }
        size_t gt_sd_i_0_max__ = estimate_r;
        for (size_t i_0__ = 0; i_0__ < gt_sd_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, gt_sd[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable gt_sd: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 608;
        if (!(context__.contains_r("bp_sd")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable bp_sd missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("bp_sd");
        pos__ = 0U;
        validate_non_negative_index("bp_sd", "(logical_gt(bp_n, 0) ? 1 : 0 )", (logical_gt(bp_n, 0) ? 1 : 0 ));
        context__.validate_dims("parameter initialization", "bp_sd", "double", context__.to_vec((logical_gt(bp_n, 0) ? 1 : 0 )));
        std::vector<double> bp_sd((logical_gt(bp_n, 0) ? 1 : 0 ), double(0));
        size_t bp_sd_k_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < bp_sd_k_0_max__; ++k_0__) {
            bp_sd[k_0__] = vals_r__[pos__++];
        }
        size_t bp_sd_i_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
        for (size_t i_0__ = 0; i_0__ < bp_sd_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, bp_sd[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable bp_sd: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 609;
        if (!(context__.contains_r("bp_effects")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable bp_effects missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("bp_effects");
        pos__ = 0U;
        validate_non_negative_index("bp_effects", "bp_n", bp_n);
        context__.validate_dims("parameter initialization", "bp_effects", "double", context__.to_vec(bp_n));
        std::vector<double> bp_effects(bp_n, double(0));
        size_t bp_effects_k_0_max__ = bp_n;
        for (size_t k_0__ = 0; k_0__ < bp_effects_k_0_max__; ++k_0__) {
            bp_effects[k_0__] = vals_r__[pos__++];
        }
        size_t bp_effects_i_0_max__ = bp_n;
        for (size_t i_0__ = 0; i_0__ < bp_effects_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_unconstrain(bp_effects[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable bp_effects: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 611;
        if (!(context__.contains_r("delay_mean")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable delay_mean missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("delay_mean");
        pos__ = 0U;
        validate_non_negative_index("delay_mean", "delays", delays);
        context__.validate_dims("parameter initialization", "delay_mean", "double", context__.to_vec(delays));
        std::vector<double> delay_mean(delays, double(0));
        size_t delay_mean_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_mean_k_0_max__; ++k_0__) {
            delay_mean[k_0__] = vals_r__[pos__++];
        }
        size_t delay_mean_i_0_max__ = delays;
        for (size_t i_0__ = 0; i_0__ < delay_mean_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_unconstrain(delay_mean[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable delay_mean: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 612;
        if (!(context__.contains_r("delay_sd")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable delay_sd missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("delay_sd");
        pos__ = 0U;
        validate_non_negative_index("delay_sd", "delays", delays);
        context__.validate_dims("parameter initialization", "delay_sd", "double", context__.to_vec(delays));
        std::vector<double> delay_sd(delays, double(0));
        size_t delay_sd_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_sd_k_0_max__; ++k_0__) {
            delay_sd[k_0__] = vals_r__[pos__++];
        }
        size_t delay_sd_i_0_max__ = delays;
        for (size_t i_0__ = 0; i_0__ < delay_sd_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, delay_sd[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable delay_sd: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 613;
        if (!(context__.contains_r("day_of_week_simplex")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable day_of_week_simplex missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("day_of_week_simplex");
        pos__ = 0U;
        validate_non_negative_index("day_of_week_simplex", "week_effect", week_effect);
        context__.validate_dims("parameter initialization", "day_of_week_simplex", "vector_d", context__.to_vec(week_effect));
        Eigen::Matrix<double, Eigen::Dynamic, 1> day_of_week_simplex(week_effect);
        size_t day_of_week_simplex_j_1_max__ = week_effect;
        for (size_t j_1__ = 0; j_1__ < day_of_week_simplex_j_1_max__; ++j_1__) {
            day_of_week_simplex(j_1__) = vals_r__[pos__++];
        }
        try {
            writer__.simplex_unconstrain(day_of_week_simplex);
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable day_of_week_simplex: ") + e.what()), current_statement_begin__, prog_reader__());
        }

        current_statement_begin__ = 614;
        if (!(context__.contains_r("frac_obs")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable frac_obs missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("frac_obs");
        pos__ = 0U;
        validate_non_negative_index("frac_obs", "obs_scale", obs_scale);
        context__.validate_dims("parameter initialization", "frac_obs", "double", context__.to_vec(obs_scale));
        std::vector<double> frac_obs(obs_scale, double(0));
        size_t frac_obs_k_0_max__ = obs_scale;
        for (size_t k_0__ = 0; k_0__ < frac_obs_k_0_max__; ++k_0__) {
            frac_obs[k_0__] = vals_r__[pos__++];
        }
        size_t frac_obs_i_0_max__ = obs_scale;
        for (size_t i_0__ = 0; i_0__ < frac_obs_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, frac_obs[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable frac_obs: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 615;
        if (!(context__.contains_r("truncation_mean")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable truncation_mean missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("truncation_mean");
        pos__ = 0U;
        validate_non_negative_index("truncation_mean", "truncation", truncation);
        context__.validate_dims("parameter initialization", "truncation_mean", "double", context__.to_vec(truncation));
        std::vector<double> truncation_mean(truncation, double(0));
        size_t truncation_mean_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_mean_k_0_max__; ++k_0__) {
            truncation_mean[k_0__] = vals_r__[pos__++];
        }
        size_t truncation_mean_i_0_max__ = truncation;
        for (size_t i_0__ = 0; i_0__ < truncation_mean_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_unconstrain(truncation_mean[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable truncation_mean: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 616;
        if (!(context__.contains_r("truncation_sd")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable truncation_sd missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("truncation_sd");
        pos__ = 0U;
        validate_non_negative_index("truncation_sd", "truncation", truncation);
        context__.validate_dims("parameter initialization", "truncation_sd", "double", context__.to_vec(truncation));
        std::vector<double> truncation_sd(truncation, double(0));
        size_t truncation_sd_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_sd_k_0_max__; ++k_0__) {
            truncation_sd[k_0__] = vals_r__[pos__++];
        }
        size_t truncation_sd_i_0_max__ = truncation;
        for (size_t i_0__ = 0; i_0__ < truncation_sd_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, truncation_sd[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable truncation_sd: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 617;
        if (!(context__.contains_r("rep_phi")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable rep_phi missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("rep_phi");
        pos__ = 0U;
        validate_non_negative_index("rep_phi", "model_type", model_type);
        context__.validate_dims("parameter initialization", "rep_phi", "double", context__.to_vec(model_type));
        std::vector<double> rep_phi(model_type, double(0));
        size_t rep_phi_k_0_max__ = model_type;
        for (size_t k_0__ = 0; k_0__ < rep_phi_k_0_max__; ++k_0__) {
            rep_phi[k_0__] = vals_r__[pos__++];
        }
        size_t rep_phi_i_0_max__ = model_type;
        for (size_t i_0__ = 0; i_0__ < rep_phi_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_lb_unconstrain(0, rep_phi[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable rep_phi: ") + 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__ = 599;
            std::vector<local_scalar_t__> rho;
            size_t rho_d_0_max__ = (fixed ? 0 : 1 );
            rho.reserve(rho_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < rho_d_0_max__; ++d_0__) {
                if (jacobian__)
                    rho.push_back(in__.scalar_lub_constrain(ls_min, ls_max, lp__));
                else
                    rho.push_back(in__.scalar_lub_constrain(ls_min, ls_max));
            }

            current_statement_begin__ = 600;
            std::vector<local_scalar_t__> alpha;
            size_t alpha_d_0_max__ = (fixed ? 0 : 1 );
            alpha.reserve(alpha_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < alpha_d_0_max__; ++d_0__) {
                if (jacobian__)
                    alpha.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    alpha.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 601;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> eta;
            (void) eta;  // dummy to suppress unused var warning
            if (jacobian__)
                eta = in__.vector_constrain((fixed ? 0 : M ), lp__);
            else
                eta = in__.vector_constrain((fixed ? 0 : M ));

            current_statement_begin__ = 603;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> log_R;
            (void) log_R;  // dummy to suppress unused var warning
            if (jacobian__)
                log_R = in__.vector_constrain(estimate_r, lp__);
            else
                log_R = in__.vector_constrain(estimate_r);

            current_statement_begin__ = 604;
            std::vector<local_scalar_t__> initial_infections;
            size_t initial_infections_d_0_max__ = estimate_r;
            initial_infections.reserve(initial_infections_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < initial_infections_d_0_max__; ++d_0__) {
                if (jacobian__)
                    initial_infections.push_back(in__.scalar_constrain(lp__));
                else
                    initial_infections.push_back(in__.scalar_constrain());
            }

            current_statement_begin__ = 605;
            std::vector<local_scalar_t__> initial_growth;
            size_t initial_growth_d_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
            initial_growth.reserve(initial_growth_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < initial_growth_d_0_max__; ++d_0__) {
                if (jacobian__)
                    initial_growth.push_back(in__.scalar_constrain(lp__));
                else
                    initial_growth.push_back(in__.scalar_constrain());
            }

            current_statement_begin__ = 606;
            std::vector<local_scalar_t__> gt_mean;
            size_t gt_mean_d_0_max__ = estimate_r;
            gt_mean.reserve(gt_mean_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < gt_mean_d_0_max__; ++d_0__) {
                if (jacobian__)
                    gt_mean.push_back(in__.scalar_lub_constrain(0, max_gt, lp__));
                else
                    gt_mean.push_back(in__.scalar_lub_constrain(0, max_gt));
            }

            current_statement_begin__ = 607;
            std::vector<local_scalar_t__> gt_sd;
            size_t gt_sd_d_0_max__ = estimate_r;
            gt_sd.reserve(gt_sd_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < gt_sd_d_0_max__; ++d_0__) {
                if (jacobian__)
                    gt_sd.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    gt_sd.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 608;
            std::vector<local_scalar_t__> bp_sd;
            size_t bp_sd_d_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
            bp_sd.reserve(bp_sd_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < bp_sd_d_0_max__; ++d_0__) {
                if (jacobian__)
                    bp_sd.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    bp_sd.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 609;
            std::vector<local_scalar_t__> bp_effects;
            size_t bp_effects_d_0_max__ = bp_n;
            bp_effects.reserve(bp_effects_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < bp_effects_d_0_max__; ++d_0__) {
                if (jacobian__)
                    bp_effects.push_back(in__.scalar_constrain(lp__));
                else
                    bp_effects.push_back(in__.scalar_constrain());
            }

            current_statement_begin__ = 611;
            std::vector<local_scalar_t__> delay_mean;
            size_t delay_mean_d_0_max__ = delays;
            delay_mean.reserve(delay_mean_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < delay_mean_d_0_max__; ++d_0__) {
                if (jacobian__)
                    delay_mean.push_back(in__.scalar_constrain(lp__));
                else
                    delay_mean.push_back(in__.scalar_constrain());
            }

            current_statement_begin__ = 612;
            std::vector<local_scalar_t__> delay_sd;
            size_t delay_sd_d_0_max__ = delays;
            delay_sd.reserve(delay_sd_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < delay_sd_d_0_max__; ++d_0__) {
                if (jacobian__)
                    delay_sd.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    delay_sd.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 613;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> day_of_week_simplex;
            (void) day_of_week_simplex;  // dummy to suppress unused var warning
            if (jacobian__)
                day_of_week_simplex = in__.simplex_constrain(week_effect, lp__);
            else
                day_of_week_simplex = in__.simplex_constrain(week_effect);

            current_statement_begin__ = 614;
            std::vector<local_scalar_t__> frac_obs;
            size_t frac_obs_d_0_max__ = obs_scale;
            frac_obs.reserve(frac_obs_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < frac_obs_d_0_max__; ++d_0__) {
                if (jacobian__)
                    frac_obs.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    frac_obs.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 615;
            std::vector<local_scalar_t__> truncation_mean;
            size_t truncation_mean_d_0_max__ = truncation;
            truncation_mean.reserve(truncation_mean_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < truncation_mean_d_0_max__; ++d_0__) {
                if (jacobian__)
                    truncation_mean.push_back(in__.scalar_constrain(lp__));
                else
                    truncation_mean.push_back(in__.scalar_constrain());
            }

            current_statement_begin__ = 616;
            std::vector<local_scalar_t__> truncation_sd;
            size_t truncation_sd_d_0_max__ = truncation;
            truncation_sd.reserve(truncation_sd_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < truncation_sd_d_0_max__; ++d_0__) {
                if (jacobian__)
                    truncation_sd.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    truncation_sd.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 617;
            std::vector<local_scalar_t__> rep_phi;
            size_t rep_phi_d_0_max__ = model_type;
            rep_phi.reserve(rep_phi_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < rep_phi_d_0_max__; ++d_0__) {
                if (jacobian__)
                    rep_phi.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    rep_phi.push_back(in__.scalar_lb_constrain(0));
            }

            // transformed parameters
            current_statement_begin__ = 621;
            validate_non_negative_index("noise", "(fixed ? 0 : noise_terms )", (fixed ? 0 : noise_terms ));
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> noise((fixed ? 0 : noise_terms ));
            stan::math::initialize(noise, DUMMY_VAR__);
            stan::math::fill(noise, DUMMY_VAR__);

            current_statement_begin__ = 622;
            validate_non_negative_index("R", "(logical_gt(estimate_r, 0) ? ot_h : 0 )", (logical_gt(estimate_r, 0) ? ot_h : 0 ));
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> R((logical_gt(estimate_r, 0) ? ot_h : 0 ));
            stan::math::initialize(R, DUMMY_VAR__);
            stan::math::fill(R, DUMMY_VAR__);

            current_statement_begin__ = 623;
            validate_non_negative_index("infections", "t", t);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> infections(t);
            stan::math::initialize(infections, DUMMY_VAR__);
            stan::math::fill(infections, DUMMY_VAR__);

            current_statement_begin__ = 624;
            validate_non_negative_index("reports", "ot_h", ot_h);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> reports(ot_h);
            stan::math::initialize(reports, DUMMY_VAR__);
            stan::math::fill(reports, DUMMY_VAR__);

            current_statement_begin__ = 625;
            validate_non_negative_index("obs_reports", "ot", ot);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> obs_reports(ot);
            stan::math::initialize(obs_reports, DUMMY_VAR__);
            stan::math::fill(obs_reports, DUMMY_VAR__);

            // transformed parameters block statements
            current_statement_begin__ = 627;
            if (as_bool(logical_negation(fixed))) {

                current_statement_begin__ = 628;
                stan::math::assign(noise, update_gp(PHI, M, L, get_base1(alpha, 1, "alpha", 1), get_base1(rho, 1, "rho", 1), eta, gp_type, pstream__));
            }
            current_statement_begin__ = 631;
            if (as_bool(estimate_r)) {

                current_statement_begin__ = 633;
                stan::math::assign(R, update_Rt(R, get_base1(log_R, estimate_r, "log_R", 1), noise, breakpoints, bp_effects, stationary, pstream__));
                current_statement_begin__ = 634;
                stan::math::assign(infections, generate_infections(R, seeding_time, gt_mean, gt_sd, max_gt, initial_infections, initial_growth, pop, future_time, pstream__));
            } else {

                current_statement_begin__ = 639;
                stan::math::assign(infections, deconvolve_infections(shifted_cases, noise, fixed, backcalc_prior, pstream__));
            }
            current_statement_begin__ = 642;
            stan::math::assign(reports, convolve_to_report(infections, delay_mean, delay_sd, max_delay, seeding_time, pstream__));
            current_statement_begin__ = 644;
            if (as_bool(logical_gt(week_effect, 1))) {

                current_statement_begin__ = 645;
                stan::math::assign(reports, day_of_week_effect(reports, day_of_week, day_of_week_simplex, pstream__));
            }
            current_statement_begin__ = 648;
            if (as_bool(obs_scale)) {

                current_statement_begin__ = 649;
                stan::math::assign(reports, scale_obs(reports, get_base1(frac_obs, 1, "frac_obs", 1), pstream__));
            }
            current_statement_begin__ = 652;
            stan::math::assign(obs_reports, truncate(stan::model::rvalue(reports, stan::model::cons_list(stan::model::index_min_max(1, ot), stan::model::nil_index_list()), "reports"), truncation_mean, truncation_sd, max_truncation, 0, pstream__));

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

            current_statement_begin__ = 621;
            size_t noise_j_1_max__ = (fixed ? 0 : noise_terms );
            for (size_t j_1__ = 0; j_1__ < noise_j_1_max__; ++j_1__) {
                if (stan::math::is_uninitialized(noise(j_1__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: noise" << "(" << j_1__ << ")";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable noise: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            current_statement_begin__ = 622;
            size_t R_j_1_max__ = (logical_gt(estimate_r, 0) ? ot_h : 0 );
            for (size_t j_1__ = 0; j_1__ < R_j_1_max__; ++j_1__) {
                if (stan::math::is_uninitialized(R(j_1__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: R" << "(" << j_1__ << ")";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable R: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            check_less_or_equal(function__, "R", R, (10 * r_mean));

            current_statement_begin__ = 623;
            size_t infections_j_1_max__ = t;
            for (size_t j_1__ = 0; j_1__ < infections_j_1_max__; ++j_1__) {
                if (stan::math::is_uninitialized(infections(j_1__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: infections" << "(" << j_1__ << ")";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable infections: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            current_statement_begin__ = 624;
            size_t reports_j_1_max__ = ot_h;
            for (size_t j_1__ = 0; j_1__ < reports_j_1_max__; ++j_1__) {
                if (stan::math::is_uninitialized(reports(j_1__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: reports" << "(" << j_1__ << ")";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable reports: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            current_statement_begin__ = 625;
            size_t obs_reports_j_1_max__ = ot;
            for (size_t j_1__ = 0; j_1__ < obs_reports_j_1_max__; ++j_1__) {
                if (stan::math::is_uninitialized(obs_reports(j_1__))) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: obs_reports" << "(" << j_1__ << ")";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable obs_reports: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }

            // model body

            current_statement_begin__ = 658;
            if (as_bool(logical_negation(fixed))) {

                current_statement_begin__ = 659;
                gaussian_process_lp(get_base1(rho, 1, "rho", 1), get_base1(alpha, 1, "alpha", 1), eta, ls_meanlog, ls_sdlog, ls_min, ls_max, alpha_sd, lp__, lp_accum__, pstream__);
            }
            current_statement_begin__ = 663;
            delays_lp(delay_mean, delay_mean_mean, delay_mean_sd, delay_sd, delay_sd_mean, delay_sd_sd, t, lp__, lp_accum__, pstream__);
            current_statement_begin__ = 665;
            truncation_lp(truncation_mean, truncation_sd, trunc_mean_mean, trunc_mean_sd, trunc_sd_mean, trunc_sd_sd, lp__, lp_accum__, pstream__);
            current_statement_begin__ = 667;
            if (as_bool(estimate_r)) {

                current_statement_begin__ = 669;
                rt_lp(log_R, initial_infections, initial_growth, bp_effects, bp_sd, bp_n, seeding_time, r_logmean, r_logsd, prior_infections, prior_growth, lp__, lp_accum__, pstream__);
                current_statement_begin__ = 672;
                generation_time_lp(gt_mean, gt_mean_mean, gt_mean_sd, gt_sd, gt_sd_mean, gt_sd_sd, ot, lp__, lp_accum__, pstream__);
            }
            current_statement_begin__ = 675;
            if (as_bool(obs_scale)) {

                current_statement_begin__ = 676;
                lp_accum__.add(normal_log<propto__>(get_base1(frac_obs, 1, "frac_obs", 1), obs_scale_mean, obs_scale_sd));
                if (get_base1(frac_obs, 1, "frac_obs", 1) < 0) lp_accum__.add(-std::numeric_limits<double>::infinity());
                else lp_accum__.add(-normal_ccdf_log(0, obs_scale_mean, obs_scale_sd));
            }
            current_statement_begin__ = 679;
            if (as_bool(likelihood)) {

                current_statement_begin__ = 680;
                report_lp(cases, obs_reports, rep_phi, phi_mean, phi_sd, model_type, obs_weight, lp__, lp_accum__, pstream__);
            }

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

        lp_accum__.add(lp__);
        return lp_accum__.sum();

    } // log_prob()

    template <bool propto, bool jacobian, typename T_>
    T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
               std::ostream* pstream = 0) const {
      std::vector<T_> vec_params_r;
      vec_params_r.reserve(params_r.size());
      for (int i = 0; i < params_r.size(); ++i)
        vec_params_r.push_back(params_r(i));
      std::vector<int> vec_params_i;
      return log_prob<propto,jacobian,T_>(vec_params_r, vec_params_i, pstream);
    }


    void get_param_names(std::vector<std::string>& names__) const {
        names__.resize(0);
        names__.push_back("rho");
        names__.push_back("alpha");
        names__.push_back("eta");
        names__.push_back("log_R");
        names__.push_back("initial_infections");
        names__.push_back("initial_growth");
        names__.push_back("gt_mean");
        names__.push_back("gt_sd");
        names__.push_back("bp_sd");
        names__.push_back("bp_effects");
        names__.push_back("delay_mean");
        names__.push_back("delay_sd");
        names__.push_back("day_of_week_simplex");
        names__.push_back("frac_obs");
        names__.push_back("truncation_mean");
        names__.push_back("truncation_sd");
        names__.push_back("rep_phi");
        names__.push_back("noise");
        names__.push_back("R");
        names__.push_back("infections");
        names__.push_back("reports");
        names__.push_back("obs_reports");
        names__.push_back("imputed_reports");
        names__.push_back("gen_R");
        names__.push_back("r");
        names__.push_back("log_lik");
    }


    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((fixed ? 0 : 1 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((fixed ? 0 : 1 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((fixed ? 0 : M ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(estimate_r);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(estimate_r);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(estimate_r);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(estimate_r);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(bp_n, 0) ? 1 : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(bp_n);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(delays);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(delays);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(week_effect);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(obs_scale);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(truncation);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(truncation);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(model_type);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((fixed ? 0 : noise_terms ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(estimate_r, 0) ? ot_h : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(t);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(ot_h);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(ot);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(ot_h);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(estimate_r, 0) ? 0 : ot_h ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(ot_h);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(return_likelihood, 1) ? ot : 0 ));
        dimss__.push_back(dims__);
    }

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

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

        // read-transform, write parameters
        std::vector<double> rho;
        size_t rho_d_0_max__ = (fixed ? 0 : 1 );
        rho.reserve(rho_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < rho_d_0_max__; ++d_0__) {
            rho.push_back(in__.scalar_lub_constrain(ls_min, ls_max));
        }
        size_t rho_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < rho_k_0_max__; ++k_0__) {
            vars__.push_back(rho[k_0__]);
        }

        std::vector<double> alpha;
        size_t alpha_d_0_max__ = (fixed ? 0 : 1 );
        alpha.reserve(alpha_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < alpha_d_0_max__; ++d_0__) {
            alpha.push_back(in__.scalar_lb_constrain(0));
        }
        size_t alpha_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < alpha_k_0_max__; ++k_0__) {
            vars__.push_back(alpha[k_0__]);
        }

        Eigen::Matrix<double, Eigen::Dynamic, 1> eta = in__.vector_constrain((fixed ? 0 : M ));
        size_t eta_j_1_max__ = (fixed ? 0 : M );
        for (size_t j_1__ = 0; j_1__ < eta_j_1_max__; ++j_1__) {
            vars__.push_back(eta(j_1__));
        }

        Eigen::Matrix<double, Eigen::Dynamic, 1> log_R = in__.vector_constrain(estimate_r);
        size_t log_R_j_1_max__ = estimate_r;
        for (size_t j_1__ = 0; j_1__ < log_R_j_1_max__; ++j_1__) {
            vars__.push_back(log_R(j_1__));
        }

        std::vector<double> initial_infections;
        size_t initial_infections_d_0_max__ = estimate_r;
        initial_infections.reserve(initial_infections_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < initial_infections_d_0_max__; ++d_0__) {
            initial_infections.push_back(in__.scalar_constrain());
        }
        size_t initial_infections_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < initial_infections_k_0_max__; ++k_0__) {
            vars__.push_back(initial_infections[k_0__]);
        }

        std::vector<double> initial_growth;
        size_t initial_growth_d_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
        initial_growth.reserve(initial_growth_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < initial_growth_d_0_max__; ++d_0__) {
            initial_growth.push_back(in__.scalar_constrain());
        }
        size_t initial_growth_k_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < initial_growth_k_0_max__; ++k_0__) {
            vars__.push_back(initial_growth[k_0__]);
        }

        std::vector<double> gt_mean;
        size_t gt_mean_d_0_max__ = estimate_r;
        gt_mean.reserve(gt_mean_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < gt_mean_d_0_max__; ++d_0__) {
            gt_mean.push_back(in__.scalar_lub_constrain(0, max_gt));
        }
        size_t gt_mean_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_mean_k_0_max__; ++k_0__) {
            vars__.push_back(gt_mean[k_0__]);
        }

        std::vector<double> gt_sd;
        size_t gt_sd_d_0_max__ = estimate_r;
        gt_sd.reserve(gt_sd_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < gt_sd_d_0_max__; ++d_0__) {
            gt_sd.push_back(in__.scalar_lb_constrain(0));
        }
        size_t gt_sd_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_sd_k_0_max__; ++k_0__) {
            vars__.push_back(gt_sd[k_0__]);
        }

        std::vector<double> bp_sd;
        size_t bp_sd_d_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
        bp_sd.reserve(bp_sd_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < bp_sd_d_0_max__; ++d_0__) {
            bp_sd.push_back(in__.scalar_lb_constrain(0));
        }
        size_t bp_sd_k_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < bp_sd_k_0_max__; ++k_0__) {
            vars__.push_back(bp_sd[k_0__]);
        }

        std::vector<double> bp_effects;
        size_t bp_effects_d_0_max__ = bp_n;
        bp_effects.reserve(bp_effects_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < bp_effects_d_0_max__; ++d_0__) {
            bp_effects.push_back(in__.scalar_constrain());
        }
        size_t bp_effects_k_0_max__ = bp_n;
        for (size_t k_0__ = 0; k_0__ < bp_effects_k_0_max__; ++k_0__) {
            vars__.push_back(bp_effects[k_0__]);
        }

        std::vector<double> delay_mean;
        size_t delay_mean_d_0_max__ = delays;
        delay_mean.reserve(delay_mean_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < delay_mean_d_0_max__; ++d_0__) {
            delay_mean.push_back(in__.scalar_constrain());
        }
        size_t delay_mean_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_mean_k_0_max__; ++k_0__) {
            vars__.push_back(delay_mean[k_0__]);
        }

        std::vector<double> delay_sd;
        size_t delay_sd_d_0_max__ = delays;
        delay_sd.reserve(delay_sd_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < delay_sd_d_0_max__; ++d_0__) {
            delay_sd.push_back(in__.scalar_lb_constrain(0));
        }
        size_t delay_sd_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_sd_k_0_max__; ++k_0__) {
            vars__.push_back(delay_sd[k_0__]);
        }

        Eigen::Matrix<double, Eigen::Dynamic, 1> day_of_week_simplex = in__.simplex_constrain(week_effect);
        size_t day_of_week_simplex_j_1_max__ = week_effect;
        for (size_t j_1__ = 0; j_1__ < day_of_week_simplex_j_1_max__; ++j_1__) {
            vars__.push_back(day_of_week_simplex(j_1__));
        }

        std::vector<double> frac_obs;
        size_t frac_obs_d_0_max__ = obs_scale;
        frac_obs.reserve(frac_obs_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < frac_obs_d_0_max__; ++d_0__) {
            frac_obs.push_back(in__.scalar_lb_constrain(0));
        }
        size_t frac_obs_k_0_max__ = obs_scale;
        for (size_t k_0__ = 0; k_0__ < frac_obs_k_0_max__; ++k_0__) {
            vars__.push_back(frac_obs[k_0__]);
        }

        std::vector<double> truncation_mean;
        size_t truncation_mean_d_0_max__ = truncation;
        truncation_mean.reserve(truncation_mean_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < truncation_mean_d_0_max__; ++d_0__) {
            truncation_mean.push_back(in__.scalar_constrain());
        }
        size_t truncation_mean_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_mean_k_0_max__; ++k_0__) {
            vars__.push_back(truncation_mean[k_0__]);
        }

        std::vector<double> truncation_sd;
        size_t truncation_sd_d_0_max__ = truncation;
        truncation_sd.reserve(truncation_sd_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < truncation_sd_d_0_max__; ++d_0__) {
            truncation_sd.push_back(in__.scalar_lb_constrain(0));
        }
        size_t truncation_sd_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_sd_k_0_max__; ++k_0__) {
            vars__.push_back(truncation_sd[k_0__]);
        }

        std::vector<double> rep_phi;
        size_t rep_phi_d_0_max__ = model_type;
        rep_phi.reserve(rep_phi_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < rep_phi_d_0_max__; ++d_0__) {
            rep_phi.push_back(in__.scalar_lb_constrain(0));
        }
        size_t rep_phi_k_0_max__ = model_type;
        for (size_t k_0__ = 0; k_0__ < rep_phi_k_0_max__; ++k_0__) {
            vars__.push_back(rep_phi[k_0__]);
        }

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

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

        if (!include_tparams__ && !include_gqs__) return;

        try {
            // declare and define transformed parameters
            current_statement_begin__ = 621;
            validate_non_negative_index("noise", "(fixed ? 0 : noise_terms )", (fixed ? 0 : noise_terms ));
            Eigen::Matrix<double, Eigen::Dynamic, 1> noise((fixed ? 0 : noise_terms ));
            stan::math::initialize(noise, DUMMY_VAR__);
            stan::math::fill(noise, DUMMY_VAR__);

            current_statement_begin__ = 622;
            validate_non_negative_index("R", "(logical_gt(estimate_r, 0) ? ot_h : 0 )", (logical_gt(estimate_r, 0) ? ot_h : 0 ));
            Eigen::Matrix<double, Eigen::Dynamic, 1> R((logical_gt(estimate_r, 0) ? ot_h : 0 ));
            stan::math::initialize(R, DUMMY_VAR__);
            stan::math::fill(R, DUMMY_VAR__);

            current_statement_begin__ = 623;
            validate_non_negative_index("infections", "t", t);
            Eigen::Matrix<double, Eigen::Dynamic, 1> infections(t);
            stan::math::initialize(infections, DUMMY_VAR__);
            stan::math::fill(infections, DUMMY_VAR__);

            current_statement_begin__ = 624;
            validate_non_negative_index("reports", "ot_h", ot_h);
            Eigen::Matrix<double, Eigen::Dynamic, 1> reports(ot_h);
            stan::math::initialize(reports, DUMMY_VAR__);
            stan::math::fill(reports, DUMMY_VAR__);

            current_statement_begin__ = 625;
            validate_non_negative_index("obs_reports", "ot", ot);
            Eigen::Matrix<double, Eigen::Dynamic, 1> obs_reports(ot);
            stan::math::initialize(obs_reports, DUMMY_VAR__);
            stan::math::fill(obs_reports, DUMMY_VAR__);

            // do transformed parameters statements
            current_statement_begin__ = 627;
            if (as_bool(logical_negation(fixed))) {

                current_statement_begin__ = 628;
                stan::math::assign(noise, update_gp(PHI, M, L, get_base1(alpha, 1, "alpha", 1), get_base1(rho, 1, "rho", 1), eta, gp_type, pstream__));
            }
            current_statement_begin__ = 631;
            if (as_bool(estimate_r)) {

                current_statement_begin__ = 633;
                stan::math::assign(R, update_Rt(R, get_base1(log_R, estimate_r, "log_R", 1), noise, breakpoints, bp_effects, stationary, pstream__));
                current_statement_begin__ = 634;
                stan::math::assign(infections, generate_infections(R, seeding_time, gt_mean, gt_sd, max_gt, initial_infections, initial_growth, pop, future_time, pstream__));
            } else {

                current_statement_begin__ = 639;
                stan::math::assign(infections, deconvolve_infections(shifted_cases, noise, fixed, backcalc_prior, pstream__));
            }
            current_statement_begin__ = 642;
            stan::math::assign(reports, convolve_to_report(infections, delay_mean, delay_sd, max_delay, seeding_time, pstream__));
            current_statement_begin__ = 644;
            if (as_bool(logical_gt(week_effect, 1))) {

                current_statement_begin__ = 645;
                stan::math::assign(reports, day_of_week_effect(reports, day_of_week, day_of_week_simplex, pstream__));
            }
            current_statement_begin__ = 648;
            if (as_bool(obs_scale)) {

                current_statement_begin__ = 649;
                stan::math::assign(reports, scale_obs(reports, get_base1(frac_obs, 1, "frac_obs", 1), pstream__));
            }
            current_statement_begin__ = 652;
            stan::math::assign(obs_reports, truncate(stan::model::rvalue(reports, stan::model::cons_list(stan::model::index_min_max(1, ot), stan::model::nil_index_list()), "reports"), truncation_mean, truncation_sd, max_truncation, 0, pstream__));

            if (!include_gqs__ && !include_tparams__) return;
            // validate transformed parameters
            const char* function__ = "validate transformed params";
            (void) function__;  // dummy to suppress unused var warning

            current_statement_begin__ = 622;
            check_less_or_equal(function__, "R", R, (10 * r_mean));

            // write transformed parameters
            if (include_tparams__) {
                size_t noise_j_1_max__ = (fixed ? 0 : noise_terms );
                for (size_t j_1__ = 0; j_1__ < noise_j_1_max__; ++j_1__) {
                    vars__.push_back(noise(j_1__));
                }
                size_t R_j_1_max__ = (logical_gt(estimate_r, 0) ? ot_h : 0 );
                for (size_t j_1__ = 0; j_1__ < R_j_1_max__; ++j_1__) {
                    vars__.push_back(R(j_1__));
                }
                size_t infections_j_1_max__ = t;
                for (size_t j_1__ = 0; j_1__ < infections_j_1_max__; ++j_1__) {
                    vars__.push_back(infections(j_1__));
                }
                size_t reports_j_1_max__ = ot_h;
                for (size_t j_1__ = 0; j_1__ < reports_j_1_max__; ++j_1__) {
                    vars__.push_back(reports(j_1__));
                }
                size_t obs_reports_j_1_max__ = ot;
                for (size_t j_1__ = 0; j_1__ < obs_reports_j_1_max__; ++j_1__) {
                    vars__.push_back(obs_reports(j_1__));
                }
            }
            if (!include_gqs__) return;
            // declare and define generated quantities
            current_statement_begin__ = 686;
            validate_non_negative_index("imputed_reports", "ot_h", ot_h);
            std::vector<int> imputed_reports(ot_h, int(0));
            stan::math::fill(imputed_reports, std::numeric_limits<int>::min());

            current_statement_begin__ = 687;
            validate_non_negative_index("gen_R", "(logical_gt(estimate_r, 0) ? 0 : ot_h )", (logical_gt(estimate_r, 0) ? 0 : ot_h ));
            Eigen::Matrix<double, Eigen::Dynamic, 1> gen_R((logical_gt(estimate_r, 0) ? 0 : ot_h ));
            stan::math::initialize(gen_R, DUMMY_VAR__);
            stan::math::fill(gen_R, DUMMY_VAR__);

            current_statement_begin__ = 688;
            validate_non_negative_index("r", "ot_h", ot_h);
            std::vector<double> r(ot_h, double(0));
            stan::math::initialize(r, DUMMY_VAR__);
            stan::math::fill(r, DUMMY_VAR__);

            current_statement_begin__ = 689;
            validate_non_negative_index("log_lik", "(logical_gt(return_likelihood, 1) ? ot : 0 )", (logical_gt(return_likelihood, 1) ? ot : 0 ));
            Eigen::Matrix<double, Eigen::Dynamic, 1> log_lik((logical_gt(return_likelihood, 1) ? ot : 0 ));
            stan::math::initialize(log_lik, DUMMY_VAR__);
            stan::math::fill(log_lik, DUMMY_VAR__);

            // generated quantities statements
            current_statement_begin__ = 690;
            if (as_bool(estimate_r)) {

                current_statement_begin__ = 692;
                stan::math::assign(r, R_to_growth(R, get_base1(gt_mean, 1, "gt_mean", 1), get_base1(gt_sd, 1, "gt_sd", 1), pstream__));
            } else {
                {
                current_statement_begin__ = 695;
                local_scalar_t__ gt_mean_sample(DUMMY_VAR__);
                (void) gt_mean_sample;  // dummy to suppress unused var warning
                stan::math::initialize(gt_mean_sample, DUMMY_VAR__);
                stan::math::fill(gt_mean_sample, DUMMY_VAR__);
                stan::math::assign(gt_mean_sample,normal_rng(gt_mean_mean, gt_mean_sd, base_rng__));

                current_statement_begin__ = 696;
                local_scalar_t__ gt_sd_sample(DUMMY_VAR__);
                (void) gt_sd_sample;  // dummy to suppress unused var warning
                stan::math::initialize(gt_sd_sample, DUMMY_VAR__);
                stan::math::fill(gt_sd_sample, DUMMY_VAR__);
                stan::math::assign(gt_sd_sample,normal_rng(gt_sd_mean, gt_sd_sd, base_rng__));


                current_statement_begin__ = 698;
                stan::math::assign(gen_R, calculate_Rt(infections, seeding_time, gt_mean_sample, gt_sd_sample, max_gt, rt_half_window, pstream__));
                current_statement_begin__ = 701;
                stan::math::assign(r, R_to_growth(gen_R, gt_mean_sample, gt_sd_sample, pstream__));
                }
            }
            current_statement_begin__ = 704;
            stan::math::assign(imputed_reports, report_rng(reports, rep_phi, model_type, base_rng__, pstream__));
            current_statement_begin__ = 706;
            if (as_bool(return_likelihood)) {

                current_statement_begin__ = 707;
                stan::math::assign(log_lik, report_log_lik(cases, obs_reports, rep_phi, model_type, obs_weight, pstream__));
            }

            // validate, write generated quantities
            current_statement_begin__ = 686;
            size_t imputed_reports_k_0_max__ = ot_h;
            for (size_t k_0__ = 0; k_0__ < imputed_reports_k_0_max__; ++k_0__) {
                vars__.push_back(imputed_reports[k_0__]);
            }

            current_statement_begin__ = 687;
            size_t gen_R_j_1_max__ = (logical_gt(estimate_r, 0) ? 0 : ot_h );
            for (size_t j_1__ = 0; j_1__ < gen_R_j_1_max__; ++j_1__) {
                vars__.push_back(gen_R(j_1__));
            }

            current_statement_begin__ = 688;
            size_t r_k_0_max__ = ot_h;
            for (size_t k_0__ = 0; k_0__ < r_k_0_max__; ++k_0__) {
                vars__.push_back(r[k_0__]);
            }

            current_statement_begin__ = 689;
            size_t log_lik_j_1_max__ = (logical_gt(return_likelihood, 1) ? ot : 0 );
            for (size_t j_1__ = 0; j_1__ < log_lik_j_1_max__; ++j_1__) {
                vars__.push_back(log_lik(j_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 ***");
        }
    }

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

    std::string model_name() const {
        return "model_estimate_infections";
    }


    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 rho_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < rho_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "rho" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t alpha_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < alpha_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "alpha" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t eta_j_1_max__ = (fixed ? 0 : M );
        for (size_t j_1__ = 0; j_1__ < eta_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "eta" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t log_R_j_1_max__ = estimate_r;
        for (size_t j_1__ = 0; j_1__ < log_R_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "log_R" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t initial_infections_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < initial_infections_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "initial_infections" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t initial_growth_k_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < initial_growth_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "initial_growth" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t gt_mean_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_mean_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gt_mean" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t gt_sd_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gt_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bp_sd_k_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < bp_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bp_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bp_effects_k_0_max__ = bp_n;
        for (size_t k_0__ = 0; k_0__ < bp_effects_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bp_effects" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t delay_mean_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_mean_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "delay_mean" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t delay_sd_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "delay_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t day_of_week_simplex_j_1_max__ = week_effect;
        for (size_t j_1__ = 0; j_1__ < day_of_week_simplex_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "day_of_week_simplex" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t frac_obs_k_0_max__ = obs_scale;
        for (size_t k_0__ = 0; k_0__ < frac_obs_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "frac_obs" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t truncation_mean_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_mean_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "truncation_mean" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t truncation_sd_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "truncation_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t rep_phi_k_0_max__ = model_type;
        for (size_t k_0__ = 0; k_0__ < rep_phi_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "rep_phi" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
            size_t noise_j_1_max__ = (fixed ? 0 : noise_terms );
            for (size_t j_1__ = 0; j_1__ < noise_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "noise" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t R_j_1_max__ = (logical_gt(estimate_r, 0) ? ot_h : 0 );
            for (size_t j_1__ = 0; j_1__ < R_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "R" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t infections_j_1_max__ = t;
            for (size_t j_1__ = 0; j_1__ < infections_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "infections" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t reports_j_1_max__ = ot_h;
            for (size_t j_1__ = 0; j_1__ < reports_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "reports" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t obs_reports_j_1_max__ = ot;
            for (size_t j_1__ = 0; j_1__ < obs_reports_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "obs_reports" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__) return;
        size_t imputed_reports_k_0_max__ = ot_h;
        for (size_t k_0__ = 0; k_0__ < imputed_reports_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "imputed_reports" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t gen_R_j_1_max__ = (logical_gt(estimate_r, 0) ? 0 : ot_h );
        for (size_t j_1__ = 0; j_1__ < gen_R_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gen_R" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t r_k_0_max__ = ot_h;
        for (size_t k_0__ = 0; k_0__ < r_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "r" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t log_lik_j_1_max__ = (logical_gt(return_likelihood, 1) ? ot : 0 );
        for (size_t j_1__ = 0; j_1__ < log_lik_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "log_lik" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
    }


    void unconstrained_param_names(std::vector<std::string>& param_names__,
                                   bool include_tparams__ = true,
                                   bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        size_t rho_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < rho_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "rho" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t alpha_k_0_max__ = (fixed ? 0 : 1 );
        for (size_t k_0__ = 0; k_0__ < alpha_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "alpha" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t eta_j_1_max__ = (fixed ? 0 : M );
        for (size_t j_1__ = 0; j_1__ < eta_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "eta" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t log_R_j_1_max__ = estimate_r;
        for (size_t j_1__ = 0; j_1__ < log_R_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "log_R" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t initial_infections_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < initial_infections_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "initial_infections" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t initial_growth_k_0_max__ = ((primitive_value(estimate_r) && primitive_value(logical_gt(seeding_time, 1))) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < initial_growth_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "initial_growth" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t gt_mean_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_mean_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gt_mean" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t gt_sd_k_0_max__ = estimate_r;
        for (size_t k_0__ = 0; k_0__ < gt_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gt_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bp_sd_k_0_max__ = (logical_gt(bp_n, 0) ? 1 : 0 );
        for (size_t k_0__ = 0; k_0__ < bp_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bp_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bp_effects_k_0_max__ = bp_n;
        for (size_t k_0__ = 0; k_0__ < bp_effects_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bp_effects" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t delay_mean_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_mean_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "delay_mean" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t delay_sd_k_0_max__ = delays;
        for (size_t k_0__ = 0; k_0__ < delay_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "delay_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t day_of_week_simplex_j_1_max__ = (week_effect - 1);
        for (size_t j_1__ = 0; j_1__ < day_of_week_simplex_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "day_of_week_simplex" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t frac_obs_k_0_max__ = obs_scale;
        for (size_t k_0__ = 0; k_0__ < frac_obs_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "frac_obs" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t truncation_mean_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_mean_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "truncation_mean" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t truncation_sd_k_0_max__ = truncation;
        for (size_t k_0__ = 0; k_0__ < truncation_sd_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "truncation_sd" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t rep_phi_k_0_max__ = model_type;
        for (size_t k_0__ = 0; k_0__ < rep_phi_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "rep_phi" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
            size_t noise_j_1_max__ = (fixed ? 0 : noise_terms );
            for (size_t j_1__ = 0; j_1__ < noise_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "noise" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t R_j_1_max__ = (logical_gt(estimate_r, 0) ? ot_h : 0 );
            for (size_t j_1__ = 0; j_1__ < R_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "R" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t infections_j_1_max__ = t;
            for (size_t j_1__ = 0; j_1__ < infections_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "infections" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t reports_j_1_max__ = ot_h;
            for (size_t j_1__ = 0; j_1__ < reports_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "reports" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
            size_t obs_reports_j_1_max__ = ot;
            for (size_t j_1__ = 0; j_1__ < obs_reports_j_1_max__; ++j_1__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "obs_reports" << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }

        if (!include_gqs__) return;
        size_t imputed_reports_k_0_max__ = ot_h;
        for (size_t k_0__ = 0; k_0__ < imputed_reports_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "imputed_reports" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t gen_R_j_1_max__ = (logical_gt(estimate_r, 0) ? 0 : ot_h );
        for (size_t j_1__ = 0; j_1__ < gen_R_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "gen_R" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t r_k_0_max__ = ot_h;
        for (size_t k_0__ = 0; k_0__ < r_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "r" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t log_lik_j_1_max__ = (logical_gt(return_likelihood, 1) ? ot : 0 );
        for (size_t j_1__ = 0; j_1__ < log_lik_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "log_lik" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
    }

}; // model

}  // namespace

typedef model_estimate_infections_namespace::model_estimate_infections stan_model;

#ifndef USING_R

stan::model::model_base& new_model(
        stan::io::var_context& data_context,
        unsigned int seed,
        std::ostream* msg_stream) {
  stan_model* m = new stan_model(data_context, seed, msg_stream);
  return *m;
}

#endif