// Code generated by Stan version 2.21.0

#include <stan/model/model_header.hpp>

namespace model_historylm_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_historylm");
    reader.add_event(283, 281, "end", "model_historylm");
    return reader;
}

template <typename T1__, typename T2__, typename T3__, typename T4__>
Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__, T3__, T4__>::type, Eigen::Dynamic, 1>
expand_history_param_to_individuals(const int& option,
                                        const T1__& fixed_value,
                                        const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& mu,
                                        const Eigen::Matrix<T3__, Eigen::Dynamic, 1>& sigma,
                                        const Eigen::Matrix<T4__, Eigen::Dynamic, 1>& rnd,
                                        const int& randomN,
                                        const int& link_function, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T1__, T2__, T3__, 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__ = 19;
        validate_non_negative_index("ind", "randomN", randomN);
        Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> ind(randomN);
        stan::math::initialize(ind, DUMMY_VAR__);
        stan::math::fill(ind, DUMMY_VAR__);


        current_statement_begin__ = 20;
        if (as_bool(logical_eq(option, 1))) {

            current_statement_begin__ = 21;
            stan::math::assign(ind, rep_vector(fixed_value, randomN));
            current_statement_begin__ = 22;
            return stan::math::promote_scalar<fun_return_scalar_t__>(ind);
        } else if (as_bool(logical_eq(option, 2))) {

            current_statement_begin__ = 25;
            stan::math::assign(ind, rep_vector(get_base1(mu, 1, "mu", 1), randomN));
        } else if (as_bool(logical_eq(option, 3))) {

            current_statement_begin__ = 28;
            stan::math::assign(ind, rnd);
        } else {

            current_statement_begin__ = 32;
            stan::math::assign(ind, add(get_base1(mu, 1, "mu", 1), multiply(get_base1(sigma, 1, "sigma", 1), rnd)));
        }
        current_statement_begin__ = 36;
        if (as_bool(logical_eq(link_function, 2))) {

            current_statement_begin__ = 37;
            return stan::math::promote_scalar<fun_return_scalar_t__>(stan::math::exp(ind));
        } else if (as_bool(logical_eq(link_function, 3))) {

            current_statement_begin__ = 40;
            return stan::math::promote_scalar<fun_return_scalar_t__>(inv_logit(ind));
        }
        current_statement_begin__ = 43;
        return stan::math::promote_scalar<fun_return_scalar_t__>(ind);
        }
    } 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 expand_history_param_to_individuals_functor__ {
    template <typename T1__, typename T2__, typename T3__, typename T4__>
        Eigen::Matrix<typename boost::math::tools::promote_args<T1__, T2__, T3__, T4__>::type, Eigen::Dynamic, 1>
    operator()(const int& option,
                                        const T1__& fixed_value,
                                        const Eigen::Matrix<T2__, Eigen::Dynamic, 1>& mu,
                                        const Eigen::Matrix<T3__, Eigen::Dynamic, 1>& sigma,
                                        const Eigen::Matrix<T4__, Eigen::Dynamic, 1>& rnd,
                                        const int& randomN,
                                        const int& link_function, std::ostream* pstream__) const {
        return expand_history_param_to_individuals(option, fixed_value, mu, sigma, rnd, randomN, link_function, pstream__);
    }
};

class model_historylm
  : public stan::model::model_base_crtp<model_historylm> {
private:
        int family;
        int rowsN;
        std::vector<double> duration;
        std::vector<int> istate;
        std::vector<int> is_used;
        std::vector<int> run_start;
        std::vector<double> session_tmean;
        int clearN;
        std::vector<double> clear_duration;
        int randomN;
        std::vector<int> irandom;
        int fixedN;
        matrix_d fixed;
        std::vector<double> history_starting_values;
        int tau_option;
        double fixed_tau;
        int tau_mu_size;
        int tau_sigma_size;
        int tau_rnd_size;
        std::vector<double> tau_prior;
        int mixed_state_option;
        double fixed_mixed_state;
        int mixed_state_mu_size;
        int mixed_state_sigma_size;
        int mixed_state_rnd_size;
        std::vector<double> mixed_state_prior;
        int lmN;
        int varianceN;
        std::vector<std::vector<double> > a_prior;
        std::vector<std::vector<double> > bH_prior;
        matrix_d fixed_priors;
        int lGamma;
        int lLogNormal;
        int lNormal;
        int oConstant;
        int oSingle;
        int oIndependent;
        int oPooled;
        int fSingle;
        int fPooled;
        int lfIdentity;
        int lfLog;
        int lfLogit;
public:
    model_historylm(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : model_base_crtp(0) {
        ctor_body(context__, 0, pstream__);
    }

    model_historylm(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_historylm_namespace::model_historylm";
        (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__ = 51;
            context__.validate_dims("data initialization", "family", "int", context__.to_vec());
            family = int(0);
            vals_i__ = context__.vals_i("family");
            pos__ = 0;
            family = vals_i__[pos__++];
            check_greater_or_equal(function__, "family", family, 1);
            check_less_or_equal(function__, "family", family, 3);

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

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

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

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

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

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

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

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

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

            current_statement_begin__ = 67;
            validate_non_negative_index("irandom", "rowsN", rowsN);
            context__.validate_dims("data initialization", "irandom", "int", context__.to_vec(rowsN));
            irandom = std::vector<int>(rowsN, int(0));
            vals_i__ = context__.vals_i("irandom");
            pos__ = 0;
            size_t irandom_k_0_max__ = rowsN;
            for (size_t k_0__ = 0; k_0__ < irandom_k_0_max__; ++k_0__) {
                irandom[k_0__] = vals_i__[pos__++];
            }
            size_t irandom_i_0_max__ = rowsN;
            for (size_t i_0__ = 0; i_0__ < irandom_i_0_max__; ++i_0__) {
                check_greater_or_equal(function__, "irandom[i_0__]", irandom[i_0__], 1);
                check_less_or_equal(function__, "irandom[i_0__]", irandom[i_0__], randomN);
            }

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

            current_statement_begin__ = 71;
            validate_non_negative_index("fixed", "rowsN", rowsN);
            validate_non_negative_index("fixed", "(logical_gt(fixedN, 0) ? fixedN : 0 )", (logical_gt(fixedN, 0) ? fixedN : 0 ));
            context__.validate_dims("data initialization", "fixed", "matrix_d", context__.to_vec(rowsN,(logical_gt(fixedN, 0) ? fixedN : 0 )));
            fixed = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(rowsN, (logical_gt(fixedN, 0) ? fixedN : 0 ));
            vals_r__ = context__.vals_r("fixed");
            pos__ = 0;
            size_t fixed_j_2_max__ = (logical_gt(fixedN, 0) ? fixedN : 0 );
            size_t fixed_j_1_max__ = rowsN;
            for (size_t j_2__ = 0; j_2__ < fixed_j_2_max__; ++j_2__) {
                for (size_t j_1__ = 0; j_1__ < fixed_j_1_max__; ++j_1__) {
                    fixed(j_1__, j_2__) = vals_r__[pos__++];
                }
            }

            current_statement_begin__ = 74;
            validate_non_negative_index("history_starting_values", "2", 2);
            context__.validate_dims("data initialization", "history_starting_values", "double", context__.to_vec(2));
            history_starting_values = std::vector<double>(2, double(0));
            vals_r__ = context__.vals_r("history_starting_values");
            pos__ = 0;
            size_t history_starting_values_k_0_max__ = 2;
            for (size_t k_0__ = 0; k_0__ < history_starting_values_k_0_max__; ++k_0__) {
                history_starting_values[k_0__] = vals_r__[pos__++];
            }
            size_t history_starting_values_i_0_max__ = 2;
            for (size_t i_0__ = 0; i_0__ < history_starting_values_i_0_max__; ++i_0__) {
                check_greater_or_equal(function__, "history_starting_values[i_0__]", history_starting_values[i_0__], 0);
                check_less_or_equal(function__, "history_starting_values[i_0__]", history_starting_values[i_0__], 1);
            }

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

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

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

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

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

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

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

            current_statement_begin__ = 86;
            context__.validate_dims("data initialization", "fixed_mixed_state", "double", context__.to_vec());
            fixed_mixed_state = double(0);
            vals_r__ = context__.vals_r("fixed_mixed_state");
            pos__ = 0;
            fixed_mixed_state = vals_r__[pos__++];
            check_greater_or_equal(function__, "fixed_mixed_state", fixed_mixed_state, 0);
            check_less_or_equal(function__, "fixed_mixed_state", fixed_mixed_state, 1);

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

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

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

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

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

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

            current_statement_begin__ = 97;
            validate_non_negative_index("a_prior", "lmN", lmN);
            validate_non_negative_index("a_prior", "2", 2);
            context__.validate_dims("data initialization", "a_prior", "double", context__.to_vec(lmN,2));
            a_prior = std::vector<std::vector<double> >(lmN, std::vector<double>(2, double(0)));
            vals_r__ = context__.vals_r("a_prior");
            pos__ = 0;
            size_t a_prior_k_0_max__ = lmN;
            size_t a_prior_k_1_max__ = 2;
            for (size_t k_1__ = 0; k_1__ < a_prior_k_1_max__; ++k_1__) {
                for (size_t k_0__ = 0; k_0__ < a_prior_k_0_max__; ++k_0__) {
                    a_prior[k_0__][k_1__] = vals_r__[pos__++];
                }
            }

            current_statement_begin__ = 100;
            validate_non_negative_index("bH_prior", "lmN", lmN);
            validate_non_negative_index("bH_prior", "2", 2);
            context__.validate_dims("data initialization", "bH_prior", "double", context__.to_vec(lmN,2));
            bH_prior = std::vector<std::vector<double> >(lmN, std::vector<double>(2, double(0)));
            vals_r__ = context__.vals_r("bH_prior");
            pos__ = 0;
            size_t bH_prior_k_0_max__ = lmN;
            size_t bH_prior_k_1_max__ = 2;
            for (size_t k_1__ = 0; k_1__ < bH_prior_k_1_max__; ++k_1__) {
                for (size_t k_0__ = 0; k_0__ < bH_prior_k_0_max__; ++k_0__) {
                    bH_prior[k_0__][k_1__] = vals_r__[pos__++];
                }
            }

            current_statement_begin__ = 103;
            validate_non_negative_index("fixed_priors", "(logical_gt(fixedN, 0) ? fixedN : 0 )", (logical_gt(fixedN, 0) ? fixedN : 0 ));
            validate_non_negative_index("fixed_priors", "2", 2);
            context__.validate_dims("data initialization", "fixed_priors", "matrix_d", context__.to_vec((logical_gt(fixedN, 0) ? fixedN : 0 ),2));
            fixed_priors = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>((logical_gt(fixedN, 0) ? fixedN : 0 ), 2);
            vals_r__ = context__.vals_r("fixed_priors");
            pos__ = 0;
            size_t fixed_priors_j_2_max__ = 2;
            size_t fixed_priors_j_1_max__ = (logical_gt(fixedN, 0) ? fixedN : 0 );
            for (size_t j_2__ = 0; j_2__ < fixed_priors_j_2_max__; ++j_2__) {
                for (size_t j_1__ = 0; j_1__ < fixed_priors_j_1_max__; ++j_1__) {
                    fixed_priors(j_1__, j_2__) = vals_r__[pos__++];
                }
            }


            // initialize transformed data variables
            current_statement_begin__ = 107;
            lGamma = int(0);
            stan::math::fill(lGamma, std::numeric_limits<int>::min());
            stan::math::assign(lGamma,1);

            current_statement_begin__ = 108;
            lLogNormal = int(0);
            stan::math::fill(lLogNormal, std::numeric_limits<int>::min());
            stan::math::assign(lLogNormal,2);

            current_statement_begin__ = 109;
            lNormal = int(0);
            stan::math::fill(lNormal, std::numeric_limits<int>::min());
            stan::math::assign(lNormal,3);

            current_statement_begin__ = 112;
            oConstant = int(0);
            stan::math::fill(oConstant, std::numeric_limits<int>::min());
            stan::math::assign(oConstant,1);

            current_statement_begin__ = 113;
            oSingle = int(0);
            stan::math::fill(oSingle, std::numeric_limits<int>::min());
            stan::math::assign(oSingle,2);

            current_statement_begin__ = 114;
            oIndependent = int(0);
            stan::math::fill(oIndependent, std::numeric_limits<int>::min());
            stan::math::assign(oIndependent,3);

            current_statement_begin__ = 115;
            oPooled = int(0);
            stan::math::fill(oPooled, std::numeric_limits<int>::min());
            stan::math::assign(oPooled,4);

            current_statement_begin__ = 118;
            fSingle = int(0);
            stan::math::fill(fSingle, std::numeric_limits<int>::min());
            stan::math::assign(fSingle,1);

            current_statement_begin__ = 119;
            fPooled = int(0);
            stan::math::fill(fPooled, std::numeric_limits<int>::min());
            stan::math::assign(fPooled,2);

            current_statement_begin__ = 122;
            lfIdentity = int(0);
            stan::math::fill(lfIdentity, std::numeric_limits<int>::min());
            stan::math::assign(lfIdentity,1);

            current_statement_begin__ = 123;
            lfLog = int(0);
            stan::math::fill(lfLog, std::numeric_limits<int>::min());
            stan::math::assign(lfLog,2);

            current_statement_begin__ = 124;
            lfLogit = int(0);
            stan::math::fill(lfLogit, std::numeric_limits<int>::min());
            stan::math::assign(lfLogit,3);

            // execute transformed data statements

            // validate transformed data

            // validate, set parameter ranges
            num_params_r__ = 0U;
            param_ranges_i__.clear();
            current_statement_begin__ = 129;
            validate_non_negative_index("tau_mu", "tau_mu_size", tau_mu_size);
            num_params_r__ += tau_mu_size;
            current_statement_begin__ = 130;
            validate_non_negative_index("tau_sigma", "tau_sigma_size", tau_sigma_size);
            num_params_r__ += tau_sigma_size;
            current_statement_begin__ = 131;
            validate_non_negative_index("tau_rnd", "tau_rnd_size", tau_rnd_size);
            num_params_r__ += tau_rnd_size;
            current_statement_begin__ = 134;
            validate_non_negative_index("mixed_state_mu", "mixed_state_mu_size", mixed_state_mu_size);
            num_params_r__ += mixed_state_mu_size;
            current_statement_begin__ = 135;
            validate_non_negative_index("mixed_state_sigma", "mixed_state_sigma_size", mixed_state_sigma_size);
            num_params_r__ += mixed_state_sigma_size;
            current_statement_begin__ = 136;
            validate_non_negative_index("mixed_state_rnd", "mixed_state_rnd_size", mixed_state_rnd_size);
            num_params_r__ += mixed_state_rnd_size;
            current_statement_begin__ = 140;
            validate_non_negative_index("a", "randomN", randomN);
            validate_non_negative_index("a", "lmN", lmN);
            num_params_r__ += (randomN * lmN);
            current_statement_begin__ = 143;
            validate_non_negative_index("bH_mu", "lmN", lmN);
            num_params_r__ += (1 * lmN);
            current_statement_begin__ = 144;
            validate_non_negative_index("bH_sigma", "(logical_gt(randomN, 1) ? lmN : 0 )", (logical_gt(randomN, 1) ? lmN : 0 ));
            num_params_r__ += (1 * (logical_gt(randomN, 1) ? lmN : 0 ));
            current_statement_begin__ = 145;
            validate_non_negative_index("bH_rnd", "(logical_gt(randomN, 1) ? randomN : 0 )", (logical_gt(randomN, 1) ? randomN : 0 ));
            validate_non_negative_index("bH_rnd", "(logical_gt(randomN, 1) ? lmN : 0 )", (logical_gt(randomN, 1) ? lmN : 0 ));
            num_params_r__ += ((logical_gt(randomN, 1) ? randomN : 0 ) * (logical_gt(randomN, 1) ? lmN : 0 ));
            current_statement_begin__ = 148;
            validate_non_negative_index("bF", "fixedN", fixedN);
            validate_non_negative_index("bF", "(logical_gt(fixedN, 0) ? lmN : 0 )", (logical_gt(fixedN, 0) ? lmN : 0 ));
            num_params_r__ += (fixedN * (logical_gt(fixedN, 0) ? lmN : 0 ));
            current_statement_begin__ = 151;
            validate_non_negative_index("sigma", "varianceN", varianceN);
            num_params_r__ += varianceN;
        } 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_historylm() { }


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

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

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

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

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

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

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

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

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

        current_statement_begin__ = 145;
        if (!(context__.contains_r("bH_rnd")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable bH_rnd missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("bH_rnd");
        pos__ = 0U;
        validate_non_negative_index("bH_rnd", "(logical_gt(randomN, 1) ? randomN : 0 )", (logical_gt(randomN, 1) ? randomN : 0 ));
        validate_non_negative_index("bH_rnd", "(logical_gt(randomN, 1) ? lmN : 0 )", (logical_gt(randomN, 1) ? lmN : 0 ));
        context__.validate_dims("parameter initialization", "bH_rnd", "vector_d", context__.to_vec((logical_gt(randomN, 1) ? lmN : 0 ),(logical_gt(randomN, 1) ? randomN : 0 )));
        std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1> > bH_rnd((logical_gt(randomN, 1) ? lmN : 0 ), Eigen::Matrix<double, Eigen::Dynamic, 1>((logical_gt(randomN, 1) ? randomN : 0 )));
        size_t bH_rnd_j_1_max__ = (logical_gt(randomN, 1) ? randomN : 0 );
        size_t bH_rnd_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bH_rnd_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bH_rnd_k_0_max__; ++k_0__) {
                bH_rnd[k_0__](j_1__) = vals_r__[pos__++];
            }
        }
        size_t bH_rnd_i_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t i_0__ = 0; i_0__ < bH_rnd_i_0_max__; ++i_0__) {
            try {
                writer__.vector_unconstrain(bH_rnd[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable bH_rnd: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 148;
        if (!(context__.contains_r("bF")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable bF missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("bF");
        pos__ = 0U;
        validate_non_negative_index("bF", "fixedN", fixedN);
        validate_non_negative_index("bF", "(logical_gt(fixedN, 0) ? lmN : 0 )", (logical_gt(fixedN, 0) ? lmN : 0 ));
        context__.validate_dims("parameter initialization", "bF", "row_vector_d", context__.to_vec((logical_gt(fixedN, 0) ? lmN : 0 ),fixedN));
        std::vector<Eigen::Matrix<double, 1, Eigen::Dynamic> > bF((logical_gt(fixedN, 0) ? lmN : 0 ), Eigen::Matrix<double, 1, Eigen::Dynamic>(fixedN));
        size_t bF_j_1_max__ = fixedN;
        size_t bF_k_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bF_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bF_k_0_max__; ++k_0__) {
                bF[k_0__](j_1__) = vals_r__[pos__++];
            }
        }
        size_t bF_i_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
        for (size_t i_0__ = 0; i_0__ < bF_i_0_max__; ++i_0__) {
            try {
                writer__.row_vector_unconstrain(bF[i_0__]);
            } catch (const std::exception& e) {
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable bF: ") + e.what()), current_statement_begin__, prog_reader__());
            }
        }

        current_statement_begin__ = 151;
        if (!(context__.contains_r("sigma")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable sigma missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("sigma");
        pos__ = 0U;
        validate_non_negative_index("sigma", "varianceN", varianceN);
        context__.validate_dims("parameter initialization", "sigma", "vector_d", context__.to_vec(varianceN));
        Eigen::Matrix<double, Eigen::Dynamic, 1> sigma(varianceN);
        size_t sigma_j_1_max__ = varianceN;
        for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
            sigma(j_1__) = vals_r__[pos__++];
        }
        try {
            writer__.vector_lb_unconstrain(0, sigma);
        } catch (const std::exception& e) {
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable sigma: ") + 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__ = 129;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> tau_mu;
            (void) tau_mu;  // dummy to suppress unused var warning
            if (jacobian__)
                tau_mu = in__.vector_constrain(tau_mu_size, lp__);
            else
                tau_mu = in__.vector_constrain(tau_mu_size);

            current_statement_begin__ = 130;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> tau_sigma;
            (void) tau_sigma;  // dummy to suppress unused var warning
            if (jacobian__)
                tau_sigma = in__.vector_lb_constrain(0, tau_sigma_size, lp__);
            else
                tau_sigma = in__.vector_lb_constrain(0, tau_sigma_size);

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

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

            current_statement_begin__ = 135;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> mixed_state_sigma;
            (void) mixed_state_sigma;  // dummy to suppress unused var warning
            if (jacobian__)
                mixed_state_sigma = in__.vector_lb_constrain(0, mixed_state_sigma_size, lp__);
            else
                mixed_state_sigma = in__.vector_lb_constrain(0, mixed_state_sigma_size);

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

            current_statement_begin__ = 140;
            std::vector<Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> > a;
            size_t a_d_0_max__ = lmN;
            a.reserve(a_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < a_d_0_max__; ++d_0__) {
                if (jacobian__)
                    a.push_back(in__.vector_constrain(randomN, lp__));
                else
                    a.push_back(in__.vector_constrain(randomN));
            }

            current_statement_begin__ = 143;
            std::vector<local_scalar_t__> bH_mu;
            size_t bH_mu_d_0_max__ = lmN;
            bH_mu.reserve(bH_mu_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < bH_mu_d_0_max__; ++d_0__) {
                if (jacobian__)
                    bH_mu.push_back(in__.scalar_constrain(lp__));
                else
                    bH_mu.push_back(in__.scalar_constrain());
            }

            current_statement_begin__ = 144;
            std::vector<local_scalar_t__> bH_sigma;
            size_t bH_sigma_d_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
            bH_sigma.reserve(bH_sigma_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < bH_sigma_d_0_max__; ++d_0__) {
                if (jacobian__)
                    bH_sigma.push_back(in__.scalar_lb_constrain(0, lp__));
                else
                    bH_sigma.push_back(in__.scalar_lb_constrain(0));
            }

            current_statement_begin__ = 145;
            std::vector<Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> > bH_rnd;
            size_t bH_rnd_d_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
            bH_rnd.reserve(bH_rnd_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < bH_rnd_d_0_max__; ++d_0__) {
                if (jacobian__)
                    bH_rnd.push_back(in__.vector_constrain((logical_gt(randomN, 1) ? randomN : 0 ), lp__));
                else
                    bH_rnd.push_back(in__.vector_constrain((logical_gt(randomN, 1) ? randomN : 0 )));
            }

            current_statement_begin__ = 148;
            std::vector<Eigen::Matrix<local_scalar_t__, 1, Eigen::Dynamic> > bF;
            size_t bF_d_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
            bF.reserve(bF_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < bF_d_0_max__; ++d_0__) {
                if (jacobian__)
                    bF.push_back(in__.row_vector_constrain(fixedN, lp__));
                else
                    bF.push_back(in__.row_vector_constrain(fixedN));
            }

            current_statement_begin__ = 151;
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> sigma;
            (void) sigma;  // dummy to suppress unused var warning
            if (jacobian__)
                sigma = in__.vector_lb_constrain(0, varianceN, lp__);
            else
                sigma = in__.vector_lb_constrain(0, varianceN);

            // transformed parameters
            current_statement_begin__ = 154;
            validate_non_negative_index("lm_param", "clearN", clearN);
            validate_non_negative_index("lm_param", "lmN", lmN);
            std::vector<Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> > lm_param(lmN, Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1>(clearN));
            stan::math::initialize(lm_param, DUMMY_VAR__);
            stan::math::fill(lm_param, DUMMY_VAR__);

            // transformed parameters block statements
            {
            current_statement_begin__ = 157;
            validate_non_negative_index("level", "2", 2);
            validate_non_negative_index("level", "3", 3);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> level(2, 3);
            stan::math::initialize(level, DUMMY_VAR__);
            stan::math::fill(level, DUMMY_VAR__);

            current_statement_begin__ = 158;
            validate_non_negative_index("current_history", "2", 2);
            std::vector<local_scalar_t__  > current_history(2, local_scalar_t__(DUMMY_VAR__));
            stan::math::initialize(current_history, DUMMY_VAR__);
            stan::math::fill(current_history, DUMMY_VAR__);

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

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

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

            current_statement_begin__ = 166;
            validate_non_negative_index("tau_ind", "randomN", randomN);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> tau_ind(randomN);
            stan::math::initialize(tau_ind, DUMMY_VAR__);
            stan::math::fill(tau_ind, DUMMY_VAR__);
            stan::math::assign(tau_ind,expand_history_param_to_individuals(tau_option, fixed_tau, tau_mu, tau_sigma, tau_rnd, randomN, lfLog, pstream__));

            current_statement_begin__ = 167;
            validate_non_negative_index("mixed_state_ind", "randomN", randomN);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> mixed_state_ind(randomN);
            stan::math::initialize(mixed_state_ind, DUMMY_VAR__);
            stan::math::fill(mixed_state_ind, DUMMY_VAR__);
            stan::math::assign(mixed_state_ind,expand_history_param_to_individuals(mixed_state_option, fixed_mixed_state, mixed_state_mu, mixed_state_sigma, mixed_state_rnd, randomN, lfLogit, pstream__));

            current_statement_begin__ = 168;
            validate_non_negative_index("bH_ind", "randomN", randomN);
            validate_non_negative_index("bH_ind", "lmN", lmN);
            std::vector<Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1>  > bH_ind(lmN, Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1>(randomN));
            stan::math::initialize(bH_ind, DUMMY_VAR__);
            stan::math::fill(bH_ind, DUMMY_VAR__);


            current_statement_begin__ = 169;
            for (int iLM = 1; iLM <= lmN; ++iLM) {

                current_statement_begin__ = 170;
                if (as_bool(logical_eq(randomN, 1))) {

                    current_statement_begin__ = 171;
                    stan::model::assign(bH_ind, 
                                stan::model::cons_list(stan::model::index_uni(iLM), stan::model::nil_index_list()), 
                                rep_vector(get_base1(bH_mu, iLM, "bH_mu", 1), randomN), 
                                "assigning variable bH_ind");
                } else {

                    current_statement_begin__ = 174;
                    stan::model::assign(bH_ind, 
                                stan::model::cons_list(stan::model::index_uni(iLM), stan::model::nil_index_list()), 
                                add(get_base1(bH_mu, iLM, "bH_mu", 1), multiply(get_base1(bH_sigma, iLM, "bH_sigma", 1), get_base1(bH_rnd, iLM, "bH_rnd", 1))), 
                                "assigning variable bH_ind");
                }
            }
            current_statement_begin__ = 178;
            for (int iT = 1; iT <= rowsN; ++iT) {

                current_statement_begin__ = 180;
                if (as_bool(get_base1(run_start, iT, "run_start", 1))) {

                    current_statement_begin__ = 181;
                    stan::math::assign(current_history, history_starting_values);
                    current_statement_begin__ = 182;
                    stan::math::assign(tau, (get_base1(session_tmean, iT, "session_tmean", 1) * get_base1(tau_ind, get_base1(irandom, iT, "irandom", 1), "tau_ind", 1)));
                    current_statement_begin__ = 185;
                    stan::math::assign(level, stan::math::to_matrix(stan::math::array_builder<Eigen::Matrix<local_scalar_t__, 1, Eigen::Dynamic> >().add(stan::math::to_row_vector(stan::math::array_builder<local_scalar_t__ >().add(1).add(0).add(get_base1(mixed_state_ind, get_base1(irandom, iT, "irandom", 1), "mixed_state_ind", 1)).array())).add(stan::math::to_row_vector(stan::math::array_builder<local_scalar_t__ >().add(0).add(1).add(get_base1(mixed_state_ind, get_base1(irandom, iT, "irandom", 1), "mixed_state_ind", 1)).array())).array()));
                }
                current_statement_begin__ = 190;
                if (as_bool(logical_eq(get_base1(is_used, iT, "is_used", 1), 1))) {

                    current_statement_begin__ = 192;
                    stan::math::assign(dH, (get_base1(current_history, (3 - get_base1(istate, iT, "istate", 1)), "current_history", 1) - get_base1(current_history, get_base1(istate, iT, "istate", 1), "current_history", 1)));
                    current_statement_begin__ = 195;
                    for (int iLM = 1; iLM <= lmN; ++iLM) {

                        current_statement_begin__ = 196;
                        stan::model::assign(lm_param, 
                                    stan::model::cons_list(stan::model::index_uni(iLM), stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list())), 
                                    (get_base1(get_base1(a, iLM, "a", 1), get_base1(irandom, iT, "irandom", 1), "a", 2) + (get_base1(get_base1(bH_ind, iLM, "bH_ind", 1), get_base1(irandom, iT, "irandom", 1), "bH_ind", 2) * dH)), 
                                    "assigning variable lm_param");
                        current_statement_begin__ = 197;
                        if (as_bool(logical_gt(fixedN, 0))) {

                            current_statement_begin__ = 198;
                            stan::model::assign(lm_param, 
                                        stan::model::cons_list(stan::model::index_uni(iLM), stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list())), 
                                        (stan::model::rvalue(lm_param, stan::model::cons_list(stan::model::index_uni(iLM), stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list())), "lm_param") + sum(elt_multiply(get_base1(bF, iLM, "bF", 1), get_base1(fixed, iT, "fixed", 1)))), 
                                        "assigning variable lm_param");
                        }
                    }
                    current_statement_begin__ = 202;
                    stan::math::assign(iC, (iC + 1));
                }
                current_statement_begin__ = 206;
                for (int iState = 1; iState <= 2; ++iState) {

                    current_statement_begin__ = 207;
                    stan::model::assign(current_history, 
                                stan::model::cons_list(stan::model::index_uni(iState), stan::model::nil_index_list()), 
                                (get_base1(level, iState, get_base1(istate, iT, "istate", 1), "level", 1) + ((get_base1(current_history, iState, "current_history", 1) - get_base1(level, iState, get_base1(istate, iT, "istate", 1), "level", 1)) * stan::math::exp((-(get_base1(duration, iT, "duration", 1)) / tau)))), 
                                "assigning variable current_history");
                }
            }
            }

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

            current_statement_begin__ = 154;
            size_t lm_param_k_0_max__ = lmN;
            size_t lm_param_j_1_max__ = clearN;
            for (size_t k_0__ = 0; k_0__ < lm_param_k_0_max__; ++k_0__) {
                for (size_t j_1__ = 0; j_1__ < lm_param_j_1_max__; ++j_1__) {
                    if (stan::math::is_uninitialized(lm_param[k_0__](j_1__))) {
                        std::stringstream msg__;
                        msg__ << "Undefined transformed parameter: lm_param" << "[" << k_0__ << "]" << "(" << j_1__ << ")";
                        stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable lm_param: ") + msg__.str()), current_statement_begin__, prog_reader__());
                    }
                }
            }

            // model body


            current_statement_begin__ = 214;
            if (as_bool(logical_eq(tau_option, oSingle))) {

                current_statement_begin__ = 215;
                lp_accum__.add(normal_log<propto__>(tau_mu, get_base1(tau_prior, 1, "tau_prior", 1), get_base1(tau_prior, 2, "tau_prior", 1)));
            } else if (as_bool(logical_eq(tau_option, oIndependent))) {

                current_statement_begin__ = 218;
                lp_accum__.add(normal_log<propto__>(tau_rnd, get_base1(tau_prior, 1, "tau_prior", 1), get_base1(tau_prior, 2, "tau_prior", 1)));
            } else if (as_bool(logical_eq(tau_option, oPooled))) {

                current_statement_begin__ = 221;
                lp_accum__.add(normal_log<propto__>(tau_mu, get_base1(tau_prior, 1, "tau_prior", 1), get_base1(tau_prior, 2, "tau_prior", 1)));
                current_statement_begin__ = 222;
                lp_accum__.add(exponential_log<propto__>(tau_sigma, 10));
                current_statement_begin__ = 223;
                lp_accum__.add(normal_log<propto__>(tau_rnd, 0, 1));
            }

            current_statement_begin__ = 227;
            if (as_bool(logical_eq(mixed_state_option, oSingle))) {

                current_statement_begin__ = 228;
                lp_accum__.add(normal_log<propto__>(mixed_state_mu, get_base1(mixed_state_prior, 1, "mixed_state_prior", 1), get_base1(mixed_state_prior, 2, "mixed_state_prior", 1)));
            } else if (as_bool(logical_eq(mixed_state_option, oIndependent))) {

                current_statement_begin__ = 231;
                lp_accum__.add(normal_log<propto__>(mixed_state_rnd, get_base1(mixed_state_prior, 1, "mixed_state_prior", 1), get_base1(mixed_state_prior, 2, "mixed_state_prior", 1)));
            } else if (as_bool(logical_eq(mixed_state_option, oPooled))) {

                current_statement_begin__ = 234;
                lp_accum__.add(normal_log<propto__>(mixed_state_mu, get_base1(mixed_state_prior, 1, "mixed_state_prior", 1), get_base1(mixed_state_prior, 2, "mixed_state_prior", 1)));
                current_statement_begin__ = 235;
                lp_accum__.add(exponential_log<propto__>(mixed_state_sigma, 10));
                current_statement_begin__ = 236;
                lp_accum__.add(normal_log<propto__>(mixed_state_rnd, 0, 1));
            }
            current_statement_begin__ = 241;
            for (int iLM = 1; iLM <= lmN; ++iLM) {

                current_statement_begin__ = 243;
                lp_accum__.add(normal_log<propto__>(get_base1(a, iLM, "a", 1), get_base1(get_base1(a_prior, iLM, "a_prior", 1), 1, "a_prior", 2), get_base1(get_base1(a_prior, iLM, "a_prior", 1), 2, "a_prior", 2)));
                current_statement_begin__ = 246;
                lp_accum__.add(normal_log<propto__>(get_base1(bH_mu, iLM, "bH_mu", 1), get_base1(get_base1(bH_prior, iLM, "bH_prior", 1), 1, "bH_prior", 2), get_base1(get_base1(bH_prior, iLM, "bH_prior", 1), 2, "bH_prior", 2)));
                current_statement_begin__ = 247;
                if (as_bool(logical_gt(randomN, 1))) {

                    current_statement_begin__ = 248;
                    lp_accum__.add(exponential_log<propto__>(get_base1(bH_sigma, iLM, "bH_sigma", 1), 1));
                    current_statement_begin__ = 249;
                    lp_accum__.add(normal_log<propto__>(get_base1(bH_rnd, iLM, "bH_rnd", 1), 0, 1));
                }
                current_statement_begin__ = 252;
                if (as_bool(logical_gt(fixedN, 0))) {

                    current_statement_begin__ = 253;
                    for (int iF = 1; iF <= fixedN; ++iF) {
                        current_statement_begin__ = 253;
                        lp_accum__.add(normal_log<propto__>(get_base1(get_base1(bF, iLM, "bF", 1), iF, "bF", 2), get_base1(fixed_priors, iF, 1, "fixed_priors", 1), get_base1(fixed_priors, iF, 2, "fixed_priors", 1)));
                    }
                }
            }
            current_statement_begin__ = 258;
            if (as_bool(logical_gt(varianceN, 0))) {

                current_statement_begin__ = 259;
                lp_accum__.add(exponential_log<propto__>(sigma, 1));
            }
            current_statement_begin__ = 263;
            if (as_bool(logical_eq(family, lGamma))) {

                current_statement_begin__ = 264;
                lp_accum__.add(gamma_log<propto__>(clear_duration, stan::math::exp(get_base1(lm_param, 1, "lm_param", 1)), elt_divide(1.0, stan::math::exp(get_base1(lm_param, 2, "lm_param", 1)))));
            } else if (as_bool(logical_eq(family, lLogNormal))) {

                current_statement_begin__ = 266;
                lp_accum__.add(lognormal_log<propto__>(clear_duration, stan::math::exp(get_base1(lm_param, 1, "lm_param", 1)), get_base1(sigma, 1, "sigma", 1)));
            } else if (as_bool(logical_eq(family, lNormal))) {

                current_statement_begin__ = 268;
                lp_accum__.add(normal_log<propto__>(clear_duration, get_base1(lm_param, 1, "lm_param", 1), get_base1(sigma, 1, "sigma", 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 ***");
        }

        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("tau_mu");
        names__.push_back("tau_sigma");
        names__.push_back("tau_rnd");
        names__.push_back("mixed_state_mu");
        names__.push_back("mixed_state_sigma");
        names__.push_back("mixed_state_rnd");
        names__.push_back("a");
        names__.push_back("bH_mu");
        names__.push_back("bH_sigma");
        names__.push_back("bH_rnd");
        names__.push_back("bF");
        names__.push_back("sigma");
        names__.push_back("lm_param");
        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(tau_mu_size);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(tau_sigma_size);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(tau_rnd_size);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(mixed_state_mu_size);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(mixed_state_sigma_size);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(mixed_state_rnd_size);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(lmN);
        dims__.push_back(randomN);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(lmN);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(randomN, 1) ? lmN : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(randomN, 1) ? lmN : 0 ));
        dims__.push_back((logical_gt(randomN, 1) ? randomN : 0 ));
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back((logical_gt(fixedN, 0) ? lmN : 0 ));
        dims__.push_back(fixedN);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(varianceN);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(lmN);
        dims__.push_back(clearN);
        dimss__.push_back(dims__);
        dims__.resize(0);
        dims__.push_back(clearN);
        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_historylm_namespace::write_array";
        (void) function__;  // dummy to suppress unused var warning

        // read-transform, write parameters
        Eigen::Matrix<double, Eigen::Dynamic, 1> tau_mu = in__.vector_constrain(tau_mu_size);
        size_t tau_mu_j_1_max__ = tau_mu_size;
        for (size_t j_1__ = 0; j_1__ < tau_mu_j_1_max__; ++j_1__) {
            vars__.push_back(tau_mu(j_1__));
        }

        Eigen::Matrix<double, Eigen::Dynamic, 1> tau_sigma = in__.vector_lb_constrain(0, tau_sigma_size);
        size_t tau_sigma_j_1_max__ = tau_sigma_size;
        for (size_t j_1__ = 0; j_1__ < tau_sigma_j_1_max__; ++j_1__) {
            vars__.push_back(tau_sigma(j_1__));
        }

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

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

        Eigen::Matrix<double, Eigen::Dynamic, 1> mixed_state_sigma = in__.vector_lb_constrain(0, mixed_state_sigma_size);
        size_t mixed_state_sigma_j_1_max__ = mixed_state_sigma_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_sigma_j_1_max__; ++j_1__) {
            vars__.push_back(mixed_state_sigma(j_1__));
        }

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

        std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1> > a;
        size_t a_d_0_max__ = lmN;
        a.reserve(a_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < a_d_0_max__; ++d_0__) {
            a.push_back(in__.vector_constrain(randomN));
        }
        size_t a_j_1_max__ = randomN;
        size_t a_k_0_max__ = lmN;
        for (size_t j_1__ = 0; j_1__ < a_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < a_k_0_max__; ++k_0__) {
                vars__.push_back(a[k_0__](j_1__));
            }
        }

        std::vector<double> bH_mu;
        size_t bH_mu_d_0_max__ = lmN;
        bH_mu.reserve(bH_mu_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < bH_mu_d_0_max__; ++d_0__) {
            bH_mu.push_back(in__.scalar_constrain());
        }
        size_t bH_mu_k_0_max__ = lmN;
        for (size_t k_0__ = 0; k_0__ < bH_mu_k_0_max__; ++k_0__) {
            vars__.push_back(bH_mu[k_0__]);
        }

        std::vector<double> bH_sigma;
        size_t bH_sigma_d_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        bH_sigma.reserve(bH_sigma_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < bH_sigma_d_0_max__; ++d_0__) {
            bH_sigma.push_back(in__.scalar_lb_constrain(0));
        }
        size_t bH_sigma_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t k_0__ = 0; k_0__ < bH_sigma_k_0_max__; ++k_0__) {
            vars__.push_back(bH_sigma[k_0__]);
        }

        std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1> > bH_rnd;
        size_t bH_rnd_d_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        bH_rnd.reserve(bH_rnd_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < bH_rnd_d_0_max__; ++d_0__) {
            bH_rnd.push_back(in__.vector_constrain((logical_gt(randomN, 1) ? randomN : 0 )));
        }
        size_t bH_rnd_j_1_max__ = (logical_gt(randomN, 1) ? randomN : 0 );
        size_t bH_rnd_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bH_rnd_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bH_rnd_k_0_max__; ++k_0__) {
                vars__.push_back(bH_rnd[k_0__](j_1__));
            }
        }

        std::vector<Eigen::Matrix<double, 1, Eigen::Dynamic> > bF;
        size_t bF_d_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
        bF.reserve(bF_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < bF_d_0_max__; ++d_0__) {
            bF.push_back(in__.row_vector_constrain(fixedN));
        }
        size_t bF_j_1_max__ = fixedN;
        size_t bF_k_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bF_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bF_k_0_max__; ++k_0__) {
                vars__.push_back(bF[k_0__](j_1__));
            }
        }

        Eigen::Matrix<double, Eigen::Dynamic, 1> sigma = in__.vector_lb_constrain(0, varianceN);
        size_t sigma_j_1_max__ = varianceN;
        for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
            vars__.push_back(sigma(j_1__));
        }

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

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

        if (!include_tparams__ && !include_gqs__) return;

        try {
            // declare and define transformed parameters
            current_statement_begin__ = 154;
            validate_non_negative_index("lm_param", "clearN", clearN);
            validate_non_negative_index("lm_param", "lmN", lmN);
            std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1> > lm_param(lmN, Eigen::Matrix<double, Eigen::Dynamic, 1>(clearN));
            stan::math::initialize(lm_param, DUMMY_VAR__);
            stan::math::fill(lm_param, DUMMY_VAR__);

            // do transformed parameters statements
            {
            current_statement_begin__ = 157;
            validate_non_negative_index("level", "2", 2);
            validate_non_negative_index("level", "3", 3);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, Eigen::Dynamic> level(2, 3);
            stan::math::initialize(level, DUMMY_VAR__);
            stan::math::fill(level, DUMMY_VAR__);

            current_statement_begin__ = 158;
            validate_non_negative_index("current_history", "2", 2);
            std::vector<local_scalar_t__  > current_history(2, local_scalar_t__(DUMMY_VAR__));
            stan::math::initialize(current_history, DUMMY_VAR__);
            stan::math::fill(current_history, DUMMY_VAR__);

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

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

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

            current_statement_begin__ = 166;
            validate_non_negative_index("tau_ind", "randomN", randomN);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> tau_ind(randomN);
            stan::math::initialize(tau_ind, DUMMY_VAR__);
            stan::math::fill(tau_ind, DUMMY_VAR__);
            stan::math::assign(tau_ind,expand_history_param_to_individuals(tau_option, fixed_tau, tau_mu, tau_sigma, tau_rnd, randomN, lfLog, pstream__));

            current_statement_begin__ = 167;
            validate_non_negative_index("mixed_state_ind", "randomN", randomN);
            Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1> mixed_state_ind(randomN);
            stan::math::initialize(mixed_state_ind, DUMMY_VAR__);
            stan::math::fill(mixed_state_ind, DUMMY_VAR__);
            stan::math::assign(mixed_state_ind,expand_history_param_to_individuals(mixed_state_option, fixed_mixed_state, mixed_state_mu, mixed_state_sigma, mixed_state_rnd, randomN, lfLogit, pstream__));

            current_statement_begin__ = 168;
            validate_non_negative_index("bH_ind", "randomN", randomN);
            validate_non_negative_index("bH_ind", "lmN", lmN);
            std::vector<Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1>  > bH_ind(lmN, Eigen::Matrix<local_scalar_t__, Eigen::Dynamic, 1>(randomN));
            stan::math::initialize(bH_ind, DUMMY_VAR__);
            stan::math::fill(bH_ind, DUMMY_VAR__);


            current_statement_begin__ = 169;
            for (int iLM = 1; iLM <= lmN; ++iLM) {

                current_statement_begin__ = 170;
                if (as_bool(logical_eq(randomN, 1))) {

                    current_statement_begin__ = 171;
                    stan::model::assign(bH_ind, 
                                stan::model::cons_list(stan::model::index_uni(iLM), stan::model::nil_index_list()), 
                                rep_vector(get_base1(bH_mu, iLM, "bH_mu", 1), randomN), 
                                "assigning variable bH_ind");
                } else {

                    current_statement_begin__ = 174;
                    stan::model::assign(bH_ind, 
                                stan::model::cons_list(stan::model::index_uni(iLM), stan::model::nil_index_list()), 
                                add(get_base1(bH_mu, iLM, "bH_mu", 1), multiply(get_base1(bH_sigma, iLM, "bH_sigma", 1), get_base1(bH_rnd, iLM, "bH_rnd", 1))), 
                                "assigning variable bH_ind");
                }
            }
            current_statement_begin__ = 178;
            for (int iT = 1; iT <= rowsN; ++iT) {

                current_statement_begin__ = 180;
                if (as_bool(get_base1(run_start, iT, "run_start", 1))) {

                    current_statement_begin__ = 181;
                    stan::math::assign(current_history, history_starting_values);
                    current_statement_begin__ = 182;
                    stan::math::assign(tau, (get_base1(session_tmean, iT, "session_tmean", 1) * get_base1(tau_ind, get_base1(irandom, iT, "irandom", 1), "tau_ind", 1)));
                    current_statement_begin__ = 185;
                    stan::math::assign(level, stan::math::to_matrix(stan::math::array_builder<Eigen::Matrix<local_scalar_t__, 1, Eigen::Dynamic> >().add(stan::math::to_row_vector(stan::math::array_builder<local_scalar_t__ >().add(1).add(0).add(get_base1(mixed_state_ind, get_base1(irandom, iT, "irandom", 1), "mixed_state_ind", 1)).array())).add(stan::math::to_row_vector(stan::math::array_builder<local_scalar_t__ >().add(0).add(1).add(get_base1(mixed_state_ind, get_base1(irandom, iT, "irandom", 1), "mixed_state_ind", 1)).array())).array()));
                }
                current_statement_begin__ = 190;
                if (as_bool(logical_eq(get_base1(is_used, iT, "is_used", 1), 1))) {

                    current_statement_begin__ = 192;
                    stan::math::assign(dH, (get_base1(current_history, (3 - get_base1(istate, iT, "istate", 1)), "current_history", 1) - get_base1(current_history, get_base1(istate, iT, "istate", 1), "current_history", 1)));
                    current_statement_begin__ = 195;
                    for (int iLM = 1; iLM <= lmN; ++iLM) {

                        current_statement_begin__ = 196;
                        stan::model::assign(lm_param, 
                                    stan::model::cons_list(stan::model::index_uni(iLM), stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list())), 
                                    (get_base1(get_base1(a, iLM, "a", 1), get_base1(irandom, iT, "irandom", 1), "a", 2) + (get_base1(get_base1(bH_ind, iLM, "bH_ind", 1), get_base1(irandom, iT, "irandom", 1), "bH_ind", 2) * dH)), 
                                    "assigning variable lm_param");
                        current_statement_begin__ = 197;
                        if (as_bool(logical_gt(fixedN, 0))) {

                            current_statement_begin__ = 198;
                            stan::model::assign(lm_param, 
                                        stan::model::cons_list(stan::model::index_uni(iLM), stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list())), 
                                        (stan::model::rvalue(lm_param, stan::model::cons_list(stan::model::index_uni(iLM), stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list())), "lm_param") + sum(elt_multiply(get_base1(bF, iLM, "bF", 1), get_base1(fixed, iT, "fixed", 1)))), 
                                        "assigning variable lm_param");
                        }
                    }
                    current_statement_begin__ = 202;
                    stan::math::assign(iC, (iC + 1));
                }
                current_statement_begin__ = 206;
                for (int iState = 1; iState <= 2; ++iState) {

                    current_statement_begin__ = 207;
                    stan::model::assign(current_history, 
                                stan::model::cons_list(stan::model::index_uni(iState), stan::model::nil_index_list()), 
                                (get_base1(level, iState, get_base1(istate, iT, "istate", 1), "level", 1) + ((get_base1(current_history, iState, "current_history", 1) - get_base1(level, iState, get_base1(istate, iT, "istate", 1), "level", 1)) * stan::math::exp((-(get_base1(duration, iT, "duration", 1)) / tau)))), 
                                "assigning variable current_history");
                }
            }
            }

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

            // write transformed parameters
            if (include_tparams__) {
                size_t lm_param_j_1_max__ = clearN;
                size_t lm_param_k_0_max__ = lmN;
                for (size_t j_1__ = 0; j_1__ < lm_param_j_1_max__; ++j_1__) {
                    for (size_t k_0__ = 0; k_0__ < lm_param_k_0_max__; ++k_0__) {
                        vars__.push_back(lm_param[k_0__](j_1__));
                    }
                }
            }
            if (!include_gqs__) return;
            // declare and define generated quantities
            current_statement_begin__ = 272;
            validate_non_negative_index("log_lik", "clearN", clearN);
            Eigen::Matrix<double, Eigen::Dynamic, 1> log_lik(clearN);
            stan::math::initialize(log_lik, DUMMY_VAR__);
            stan::math::fill(log_lik, DUMMY_VAR__);

            // generated quantities statements
            current_statement_begin__ = 274;
            if (as_bool(logical_eq(family, lGamma))) {

                current_statement_begin__ = 275;
                for (int iC = 1; iC <= clearN; ++iC) {
                    current_statement_begin__ = 275;
                    stan::model::assign(log_lik, 
                                stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list()), 
                                gamma_log(get_base1(clear_duration, iC, "clear_duration", 1), stan::math::exp(get_base1(get_base1(lm_param, 1, "lm_param", 1), iC, "lm_param", 2)), (1 / stan::math::exp(get_base1(get_base1(lm_param, 2, "lm_param", 1), iC, "lm_param", 2)))), 
                                "assigning variable log_lik");
                }
            } else if (as_bool(logical_eq(family, lLogNormal))) {

                current_statement_begin__ = 277;
                for (int iC = 1; iC <= clearN; ++iC) {
                    current_statement_begin__ = 277;
                    stan::model::assign(log_lik, 
                                stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list()), 
                                lognormal_log(get_base1(clear_duration, iC, "clear_duration", 1), stan::math::exp(get_base1(get_base1(lm_param, 1, "lm_param", 1), iC, "lm_param", 2)), get_base1(sigma, 1, "sigma", 1)), 
                                "assigning variable log_lik");
                }
            } else if (as_bool(logical_eq(family, lNormal))) {

                current_statement_begin__ = 279;
                for (int iC = 1; iC <= clearN; ++iC) {
                    current_statement_begin__ = 279;
                    stan::model::assign(log_lik, 
                                stan::model::cons_list(stan::model::index_uni(iC), stan::model::nil_index_list()), 
                                normal_log(get_base1(clear_duration, iC, "clear_duration", 1), get_base1(get_base1(lm_param, 1, "lm_param", 1), iC, "lm_param", 2), get_base1(sigma, 1, "sigma", 1)), 
                                "assigning variable log_lik");
                }
            }

            // validate, write generated quantities
            current_statement_begin__ = 272;
            size_t log_lik_j_1_max__ = clearN;
            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_historylm";
    }


    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 tau_mu_j_1_max__ = tau_mu_size;
        for (size_t j_1__ = 0; j_1__ < tau_mu_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau_mu" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t tau_sigma_j_1_max__ = tau_sigma_size;
        for (size_t j_1__ = 0; j_1__ < tau_sigma_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau_sigma" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t tau_rnd_j_1_max__ = tau_rnd_size;
        for (size_t j_1__ = 0; j_1__ < tau_rnd_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau_rnd" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t mixed_state_mu_j_1_max__ = mixed_state_mu_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_mu_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "mixed_state_mu" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t mixed_state_sigma_j_1_max__ = mixed_state_sigma_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_sigma_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "mixed_state_sigma" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t mixed_state_rnd_j_1_max__ = mixed_state_rnd_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_rnd_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "mixed_state_rnd" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t a_j_1_max__ = randomN;
        size_t a_k_0_max__ = lmN;
        for (size_t j_1__ = 0; j_1__ < a_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < a_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "a" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t bH_mu_k_0_max__ = lmN;
        for (size_t k_0__ = 0; k_0__ < bH_mu_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bH_mu" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bH_sigma_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t k_0__ = 0; k_0__ < bH_sigma_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bH_sigma" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bH_rnd_j_1_max__ = (logical_gt(randomN, 1) ? randomN : 0 );
        size_t bH_rnd_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bH_rnd_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bH_rnd_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "bH_rnd" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t bF_j_1_max__ = fixedN;
        size_t bF_k_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bF_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bF_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "bF" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t sigma_j_1_max__ = varianceN;
        for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "sigma" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
            size_t lm_param_j_1_max__ = clearN;
            size_t lm_param_k_0_max__ = lmN;
            for (size_t j_1__ = 0; j_1__ < lm_param_j_1_max__; ++j_1__) {
                for (size_t k_0__ = 0; k_0__ < lm_param_k_0_max__; ++k_0__) {
                    param_name_stream__.str(std::string());
                    param_name_stream__ << "lm_param" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                    param_names__.push_back(param_name_stream__.str());
                }
            }
        }

        if (!include_gqs__) return;
        size_t log_lik_j_1_max__ = clearN;
        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 tau_mu_j_1_max__ = tau_mu_size;
        for (size_t j_1__ = 0; j_1__ < tau_mu_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau_mu" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t tau_sigma_j_1_max__ = tau_sigma_size;
        for (size_t j_1__ = 0; j_1__ < tau_sigma_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau_sigma" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t tau_rnd_j_1_max__ = tau_rnd_size;
        for (size_t j_1__ = 0; j_1__ < tau_rnd_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "tau_rnd" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t mixed_state_mu_j_1_max__ = mixed_state_mu_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_mu_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "mixed_state_mu" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t mixed_state_sigma_j_1_max__ = mixed_state_sigma_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_sigma_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "mixed_state_sigma" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t mixed_state_rnd_j_1_max__ = mixed_state_rnd_size;
        for (size_t j_1__ = 0; j_1__ < mixed_state_rnd_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "mixed_state_rnd" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t a_j_1_max__ = randomN;
        size_t a_k_0_max__ = lmN;
        for (size_t j_1__ = 0; j_1__ < a_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < a_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "a" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t bH_mu_k_0_max__ = lmN;
        for (size_t k_0__ = 0; k_0__ < bH_mu_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bH_mu" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bH_sigma_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t k_0__ = 0; k_0__ < bH_sigma_k_0_max__; ++k_0__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "bH_sigma" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        size_t bH_rnd_j_1_max__ = (logical_gt(randomN, 1) ? randomN : 0 );
        size_t bH_rnd_k_0_max__ = (logical_gt(randomN, 1) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bH_rnd_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bH_rnd_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "bH_rnd" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t bF_j_1_max__ = fixedN;
        size_t bF_k_0_max__ = (logical_gt(fixedN, 0) ? lmN : 0 );
        for (size_t j_1__ = 0; j_1__ < bF_j_1_max__; ++j_1__) {
            for (size_t k_0__ = 0; k_0__ < bF_k_0_max__; ++k_0__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "bF" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
        }
        size_t sigma_j_1_max__ = varianceN;
        for (size_t j_1__ = 0; j_1__ < sigma_j_1_max__; ++j_1__) {
            param_name_stream__.str(std::string());
            param_name_stream__ << "sigma" << '.' << j_1__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
            size_t lm_param_j_1_max__ = clearN;
            size_t lm_param_k_0_max__ = lmN;
            for (size_t j_1__ = 0; j_1__ < lm_param_j_1_max__; ++j_1__) {
                for (size_t k_0__ = 0; k_0__ < lm_param_k_0_max__; ++k_0__) {
                    param_name_stream__.str(std::string());
                    param_name_stream__ << "lm_param" << '.' << k_0__ + 1 << '.' << j_1__ + 1;
                    param_names__.push_back(param_name_stream__.str());
                }
            }
        }

        if (!include_gqs__) return;
        size_t log_lik_j_1_max__ = clearN;
        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_historylm_namespace::model_historylm 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