// Copyright (c) Facebook, Inc. and its affiliates.

// This source code is licensed under the MIT license found in the
// LICENSE file in the root directory of this source tree.

functions {
  real[ , ] get_changepoint_matrix(real[] t, real[] t_change, int T, int S) {
    // Assumes t and t_change are sorted.
    real A[T, S];
    real a_row[S];
    int cp_idx;

    // Start with an empty matrix.
    A = rep_array(0, T, S);
    a_row = rep_array(0, S);
    cp_idx = 1;

    // Fill in each row of A.
    for (i in 1:T) {
      while ((cp_idx <= S) && (t[i] >= t_change[cp_idx])) {
        a_row[cp_idx] = 1;
        cp_idx = cp_idx + 1;
      }
      A[i] = a_row;
    }
    return A;
  }

  // Logistic trend functions

  real[] logistic_gamma(real k, real m, real[] delta, real[] t_change, int S) {
    real gamma[S];  // adjusted offsets, for piecewise continuity
    real k_s[S + 1];  // actual rate in each segment
    real m_pr;

    // Compute the rate in each segment
    k_s[1] = k;
    for (i in 1:S) {
      k_s[i + 1] = k_s[i] + delta[i];
    }

    // Piecewise offsets
    m_pr = m; // The offset in the previous segment
      gamma[i] = (t_change[i] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
      m_pr = m_pr + gamma[i];  // update for the next segment
    }
    return gamma;
  }
  
  real[] logistic_trend(
    real k,
    real m,
    real[] delta,
    real[] t,
    real[] cap,
    real[ , ] A,
    real[] t_change,
    int S,
    int T
  ) {
    real gamma[S];
    real Y[T];

    gamma = logistic_gamma(k, m, delta, t_change, S);
      Y[i] = cap[i] / (1 + exp(-(k + dot_product(A[i], delta))
        * (t[i] - (m + dot_product(A[i], gamma)))));
    }
    return Y;
  }

  // Linear trend function

  real[] linear_trend(
    real k,
    real m,
    real[] t,
    real[ , ] A,
    int S,
    int T
  ) {
    real Y[T];

      gamma[i] = -t_change[i] * delta[i];
    }
      Y[i] = (k + dot_product(A[i], delta)) * t[i] + (
        m + dot_product(A[i], gamma));
    }
    return Y;
  }

   // Flat trend function

    real[] flat_trend(
    real m,
    int T
  ) {
    return rep_array(m, T);
  }

}

data {
  int T;                // Number of time periods
  int<lower=1> K;       // Number of regressors
  real t[T];            // Time
  real cap[T];          // Capacities for logistic trend
  real y[T];            // Time series
  int S;                // Number of changepoints
  real t_change[S];     // Times of trend changepoints
  real X[T,K];         // Regressors
  vector[K] sigmas;     // Scale on seasonality prior
  real<lower=0> tau;    // Scale on changepoints prior
  int trend_indicator;  // 0 for linear, 1 for logistic, 2 for flat
  real s_a[K];          // Indicator of additive features
  real s_m[K];          // Indicator of multiplicative features
}

transformed data {
  real A[T, S];
  A = get_changepoint_matrix(t, t_change, T, S);
}

parameters {
  real k;                   // Base trend growth rate
  real m;                   // Trend offset
  real delta[S];            // Trend rate adjustments
  real<lower=0> sigma_obs;  // Observation noise
  real beta[K];             // Regressor coefficients
}

transformed parameters {
  real trend[T];
  real Y[T];
  real beta_m[K];
  real beta_a[K];

  if (trend_indicator == 0) {
    trend = linear_trend(k, m, delta, t, A, t_change, S, T);
  } else if (trend_indicator == 1) {
    trend = logistic_trend(k, m, delta, t, cap, A, t_change, S, T);
  } else if (trend_indicator == 2){
    trend = flat_trend(m, T);
  }

  for (i in 1:K) {
    beta_m[i] = beta[i] * s_m[i];
    beta_a[i] = beta[i] * s_a[i];
  }

  for (i in 1:T) {
    Y[i] = (
      trend[i] * (1 + dot_product(X[i], beta_m)) + dot_product(X[i], beta_a)
    );
  }
}

model {
  //priors
  k ~ normal(0, 5);
  m ~ normal(0, 5);
  delta ~ double_exponential(0, tau);
  sigma_obs ~ normal(0, 0.5);
  beta ~ normal(0, sigmas);

  // Likelihood
  y ~ normal(Y, sigma_obs);
}

// Code generated by Stan version 2.19.1

#include <stan/model/model_header.hpp>

namespace anon_model_f5236004a3fd5b8429270d00efcc0cf9_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", "unknown file name");
    reader.add_event(177, 175, "end", "unknown file name");
    return reader;
}

template <typename T0__, typename T1__>
std::vector<std::vector<typename boost::math::tools::promote_args<T0__, T1__>::type> >
get_changepoint_matrix(const std::vector<T0__>& t,
                           const std::vector<T1__>& t_change,
                           const int& T,
                           const int& S, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__>::type local_scalar_t__;
    typedef local_scalar_t__ fun_return_scalar_t__;
    const static bool propto__ = true;
    (void) propto__;
        local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
        (void) DUMMY_VAR__;  // suppress unused var warning

    int current_statement_begin__ = -1;
    try {
        {
        current_statement_begin__ = 9;
        validate_non_negative_index("A", "T", T);
        validate_non_negative_index("A", "S", S);
        std::vector<std::vector<local_scalar_t__  >  > A(T, std::vector<local_scalar_t__>(S, local_scalar_t__(DUMMY_VAR__)));
        stan::math::initialize(A, DUMMY_VAR__);
        stan::math::fill(A, DUMMY_VAR__);

        current_statement_begin__ = 10;
        validate_non_negative_index("a_row", "S", S);
        std::vector<local_scalar_t__  > a_row(S, local_scalar_t__(DUMMY_VAR__));
        stan::math::initialize(a_row, DUMMY_VAR__);
        stan::math::fill(a_row, DUMMY_VAR__);

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

        current_statement_begin__ = 14;
        stan::math::assign(A, rep_array(0, T, S));
        current_statement_begin__ = 15;
        stan::math::assign(a_row, rep_array(0, S));
        current_statement_begin__ = 16;
        stan::math::assign(cp_idx, 1);
        current_statement_begin__ = 19;
        for (int i = 1; i <= T; ++i) {

            current_statement_begin__ = 20;
            while (as_bool((primitive_value(logical_lte(cp_idx, S)) && primitive_value(logical_gte(get_base1(t, i, "t", 1), get_base1(t_change, cp_idx, "t_change", 1)))))) {

                current_statement_begin__ = 21;
                stan::model::assign(a_row, 
                            stan::model::cons_list(stan::model::index_uni(cp_idx), stan::model::nil_index_list()), 
                            1, 
                            "assigning variable a_row");
                current_statement_begin__ = 22;
                stan::math::assign(cp_idx, (cp_idx + 1));
            }
            current_statement_begin__ = 24;
            stan::model::assign(A, 
                        stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                        a_row, 
                        "assigning variable A");
        }
        current_statement_begin__ = 26;
        return stan::math::promote_scalar<fun_return_scalar_t__>(A);
        }
    } 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 get_changepoint_matrix_functor__ {
    template <typename T0__, typename T1__>
        std::vector<std::vector<typename boost::math::tools::promote_args<T0__, T1__>::type> >
    operator()(const std::vector<T0__>& t,
                           const int& S, std::ostream* pstream__) const {
        return get_changepoint_matrix(t, t_change, T, S, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__>
std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__>::type>
logistic_gamma(const T0__& k,
                   const T1__& m,
                   const std::vector<T2__>& delta,
                   const std::vector<T3__>& t_change,
                   const int& S, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__>::type local_scalar_t__;

    try {
        {
        current_statement_begin__ = 32;
        validate_non_negative_index("gamma", "S", S);
        std::vector<local_scalar_t__  > gamma(S, local_scalar_t__(DUMMY_VAR__));
        stan::math::initialize(gamma, DUMMY_VAR__);
        stan::math::fill(gamma, DUMMY_VAR__);

        current_statement_begin__ = 33;
        validate_non_negative_index("k_s", "(S + 1)", (S + 1));
        std::vector<local_scalar_t__  > k_s((S + 1), local_scalar_t__(DUMMY_VAR__));
        stan::math::initialize(k_s, DUMMY_VAR__);
        stan::math::fill(k_s, DUMMY_VAR__);

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

        current_statement_begin__ = 37;
        stan::model::assign(k_s, 
                    stan::model::cons_list(stan::model::index_uni(1), stan::model::nil_index_list()), 
                    k, 
                    "assigning variable k_s");
        current_statement_begin__ = 38;
        for (int i = 1; i <= S; ++i) {

            current_statement_begin__ = 39;
            stan::model::assign(k_s, 
                        stan::model::cons_list(stan::model::index_uni((i + 1)), stan::model::nil_index_list()), 
                        (get_base1(k_s, i, "k_s", 1) + get_base1(delta, i, "delta", 1)), 
                        "assigning variable k_s");
        }
        current_statement_begin__ = 43;
        stan::math::assign(m_pr, m);
        current_statement_begin__ = 44;

            current_statement_begin__ = 45;
            stan::model::assign(gamma, 
                        ((get_base1(t_change, i, "t_change", 1) - m_pr) * (1 - (get_base1(k_s, i, "k_s", 1) / get_base1(k_s, (i + 1), "k_s", 1)))), 
                        "assigning variable gamma");
            current_statement_begin__ = 46;
            stan::math::assign(m_pr, (m_pr + get_base1(gamma, i, "gamma", 1)));
        }
        current_statement_begin__ = 48;
        return stan::math::promote_scalar<fun_return_scalar_t__>(gamma);
        }
    }
}

struct logistic_gamma_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__>
        std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__>::type>
    operator()(const T0__& k,
                   const int& S, std::ostream* pstream__) const {
        return logistic_gamma(k, m, delta, t_change, S, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T6__>
std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T6__>::type>::type>
logistic_trend(const T0__& k,
                   const std::vector<T3__>& t,
                   const std::vector<T4__>& cap,
                   const std::vector<std::vector<T5__> >& A,
                   const std::vector<T6__>& t_change,
                   const int& S,
                   const int& T, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T6__>::type>::type local_scalar_t__;

    try {
        {
        current_statement_begin__ = 62;

        current_statement_begin__ = 63;
        validate_non_negative_index("Y", "T", T);
        std::vector<local_scalar_t__  > Y(T, local_scalar_t__(DUMMY_VAR__));
        stan::math::initialize(Y, DUMMY_VAR__);
        stan::math::fill(Y, DUMMY_VAR__);

        current_statement_begin__ = 65;
        stan::math::assign(gamma, logistic_gamma(k, m, delta, t_change, S, pstream__));
        current_statement_begin__ = 66;

            current_statement_begin__ = 67;
            stan::model::assign(Y, 
                        (get_base1(cap, i, "cap", 1) / (1 + stan::math::exp((-((k + dot_product(get_base1(A, i, "A", 1), delta))) * (get_base1(t, i, "t", 1) - (m + dot_product(get_base1(A, i, "A", 1), gamma))))))), 
                        "assigning variable Y");
        }
        current_statement_begin__ = 70;
        return stan::math::promote_scalar<fun_return_scalar_t__>(Y);
        }
    }
}

struct logistic_trend_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__, typename T6__>
        std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__, T6__>::type>::type>
                   const int& T, std::ostream* pstream__) const {
        return logistic_trend(k, m, delta, t, cap, A, t_change, S, T, pstream__);
    }
};

template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__>
std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type>
linear_trend(const T0__& k,
                 const T1__& m,
                 const std::vector<T2__>& delta,
                 const std::vector<T3__>& t,
                 const std::vector<std::vector<T4__> >& A,
                 const std::vector<T5__>& t_change,
                 const int& S,
                 const int& T, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type local_scalar_t__;

    try {
        {
        current_statement_begin__ = 85;

        current_statement_begin__ = 86;

        current_statement_begin__ = 88;

            current_statement_begin__ = 89;
                        (-(get_base1(t_change, i, "t_change", 1)) * get_base1(delta, i, "delta", 1)), 
        }
        current_statement_begin__ = 91;

            current_statement_begin__ = 92;
                        (((k + dot_product(get_base1(A, i, "A", 1), delta)) * get_base1(t, i, "t", 1)) + (m + dot_product(get_base1(A, i, "A", 1), gamma))), 
        }
        current_statement_begin__ = 95;
        }
    }
}

struct linear_trend_functor__ {
    template <typename T0__, typename T1__, typename T2__, typename T3__, typename T4__, typename T5__>
        std::vector<typename boost::math::tools::promote_args<T0__, T1__, T2__, T3__, typename boost::math::tools::promote_args<T4__, T5__>::type>::type>
                 const int& T, std::ostream* pstream__) const {
        return linear_trend(k, m, delta, t, A, t_change, S, T, pstream__);
    }
};

template <typename T0__>
std::vector<typename boost::math::tools::promote_args<T0__>::type>
flat_trend(const T0__& m,
               const int& T, std::ostream* pstream__) {
    typedef typename boost::math::tools::promote_args<T0__>::type local_scalar_t__;

    try {

        current_statement_begin__ = 104;
        return stan::math::promote_scalar<fun_return_scalar_t__>(rep_array(m, T));
    }
}

struct flat_trend_functor__ {
    template <typename T0__>
        std::vector<typename boost::math::tools::promote_args<T0__>::type>
    operator()(const T0__& m,
               const int& T, std::ostream* pstream__) const {
        return flat_trend(m, T, pstream__);
    }
};

class anon_model_f5236004a3fd5b8429270d00efcc0cf9 : public prob_grad {
private:
        int T;
        int K;
        std::vector<double> t;
        std::vector<double> cap;
        std::vector<double> y;
        int S;
        std::vector<double> t_change;
        std::vector<std::vector<double> > X;
        vector_d sigmas;
        double tau;
        int trend_indicator;
        std::vector<double> s_a;
        std::vector<double> s_m;
        std::vector<std::vector<double> > A;
public:
    anon_model_f5236004a3fd5b8429270d00efcc0cf9(stan::io::var_context& context__,
        std::ostream* pstream__ = 0)
        : prob_grad(0) {
        ctor_body(context__, 0, pstream__);
    }

        unsigned int random_seed__,
        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__ = "anon_model_f5236004a3fd5b8429270d00efcc0cf9_namespace::anon_model_f5236004a3fd5b8429270d00efcc0cf9";
        (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__;

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

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

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

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

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

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

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

            current_statement_begin__ = 118;
            validate_non_negative_index("X", "T", T);
            validate_non_negative_index("X", "K", K);
            context__.validate_dims("data initialization", "X", "double", context__.to_vec(T,K));
            X = std::vector<std::vector<double> >(T, std::vector<double>(K, double(0)));
            vals_r__ = context__.vals_r("X");
            pos__ = 0;
            size_t X_k_0_max__ = T;
            size_t X_k_1_max__ = K;
            for (size_t k_1__ = 0; k_1__ < X_k_1_max__; ++k_1__) {
                for (size_t k_0__ = 0; k_0__ < X_k_0_max__; ++k_0__) {
                    X[k_0__][k_1__] = vals_r__[pos__++];
                }
            }

            current_statement_begin__ = 119;
            validate_non_negative_index("sigmas", "K", K);
            context__.validate_dims("data initialization", "sigmas", "vector_d", context__.to_vec(K));
            sigmas = Eigen::Matrix<double, Eigen::Dynamic, 1>(K);
            vals_r__ = context__.vals_r("sigmas");
            pos__ = 0;
            size_t sigmas_j_1_max__ = K;
            for (size_t j_1__ = 0; j_1__ < sigmas_j_1_max__; ++j_1__) {
                sigmas(j_1__) = vals_r__[pos__++];
            }

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

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

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

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

            // initialize transformed data variables
            current_statement_begin__ = 127;
            validate_non_negative_index("A", "T", T);
            validate_non_negative_index("A", "S", S);
            A = std::vector<std::vector<double> >(T, std::vector<double>(S, double(0)));
            stan::math::fill(A, DUMMY_VAR__);

            // execute transformed data statements
            current_statement_begin__ = 128;
            stan::math::assign(A, get_changepoint_matrix(t, t_change, T, S, pstream__));

            // validate transformed data

            // validate, set parameter ranges
            num_params_r__ = 0U;
            param_ranges_i__.clear();
            current_statement_begin__ = 132;
            num_params_r__ += 1;
            current_statement_begin__ = 133;
            current_statement_begin__ = 134;
            validate_non_negative_index("delta", "S", S);
            num_params_r__ += (1 * S);
            current_statement_begin__ = 135;
            current_statement_begin__ = 136;
            validate_non_negative_index("beta", "K", K);
            num_params_r__ += (1 * K);
        } 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 ***");
        }
    }

    ~anon_model_f5236004a3fd5b8429270d00efcc0cf9() { }

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

        current_statement_begin__ = 132;
        if (!(context__.contains_r("k")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable k missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("k");
        pos__ = 0U;
        context__.validate_dims("parameter initialization", "k", "double", context__.to_vec());
        double k(0);
        k = vals_r__[pos__++];
        try {
            writer__.scalar_unconstrain(k);
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable k: ") + e.what()), current_statement_begin__, prog_reader__());
        }

        current_statement_begin__ = 133;
        if (!(context__.contains_r("m")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable m missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("m");
        pos__ = 0U;
        context__.validate_dims("parameter initialization", "m", "double", context__.to_vec());
        double m(0);
        m = vals_r__[pos__++];
        try {
            writer__.scalar_unconstrain(m);
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable m: ") + e.what()), current_statement_begin__, prog_reader__());
        }

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

        current_statement_begin__ = 135;
        if (!(context__.contains_r("sigma_obs")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable sigma_obs missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("sigma_obs");
        pos__ = 0U;
        context__.validate_dims("parameter initialization", "sigma_obs", "double", context__.to_vec());
        double sigma_obs(0);
        sigma_obs = vals_r__[pos__++];
        try {
            writer__.scalar_lb_unconstrain(0, sigma_obs);
            stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable sigma_obs: ") + e.what()), current_statement_begin__, prog_reader__());
        }

        current_statement_begin__ = 136;
        if (!(context__.contains_r("beta")))
            stan::lang::rethrow_located(std::runtime_error(std::string("Variable beta missing")), current_statement_begin__, prog_reader__());
        vals_r__ = context__.vals_r("beta");
        pos__ = 0U;
        validate_non_negative_index("beta", "K", K);
        context__.validate_dims("parameter initialization", "beta", "double", context__.to_vec(K));
        std::vector<double> beta(K, double(0));
        size_t beta_k_0_max__ = K;
        for (size_t k_0__ = 0; k_0__ < beta_k_0_max__; ++k_0__) {
            beta[k_0__] = vals_r__[pos__++];
        }
        size_t beta_i_0_max__ = K;
        for (size_t i_0__ = 0; i_0__ < beta_i_0_max__; ++i_0__) {
            try {
                writer__.scalar_unconstrain(beta[i_0__]);
                stan::lang::rethrow_located(std::runtime_error(std::string("Error transforming variable beta: ") + 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::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__;

        (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
            local_scalar_t__ k;
            (void) k;  // dummy to suppress unused var warning
            if (jacobian__)
                k = in__.scalar_constrain(lp__);
            else
                k = in__.scalar_constrain();

            local_scalar_t__ m;
            (void) m;  // dummy to suppress unused var warning
                m = in__.scalar_constrain(lp__);
            else
                m = in__.scalar_constrain();

            std::vector<local_scalar_t__> delta;
            size_t delta_d_0_max__ = S;
            delta.reserve(delta_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < delta_d_0_max__; ++d_0__) {
                if (jacobian__)
                    delta.push_back(in__.scalar_constrain(lp__));
                else
                    delta.push_back(in__.scalar_constrain());
            }

            local_scalar_t__ sigma_obs;
            (void) sigma_obs;  // dummy to suppress unused var warning
                sigma_obs = in__.scalar_lb_constrain(0, lp__);
            else
                sigma_obs = in__.scalar_lb_constrain(0);

            std::vector<local_scalar_t__> beta;
            size_t beta_d_0_max__ = K;
            beta.reserve(beta_d_0_max__);
            for (size_t d_0__ = 0; d_0__ < beta_d_0_max__; ++d_0__) {
                    beta.push_back(in__.scalar_constrain(lp__));
                else
                    beta.push_back(in__.scalar_constrain());
            }

            // transformed parameters
            current_statement_begin__ = 140;
            validate_non_negative_index("trend", "T", T);
            std::vector<local_scalar_t__> trend(T, local_scalar_t__(0));
            stan::math::initialize(trend, DUMMY_VAR__);
            stan::math::fill(trend, DUMMY_VAR__);

            current_statement_begin__ = 141;
            validate_non_negative_index("Y", "T", T);
            std::vector<local_scalar_t__> Y(T, local_scalar_t__(0));
            stan::math::initialize(Y, DUMMY_VAR__);
            stan::math::fill(Y, DUMMY_VAR__);

            current_statement_begin__ = 142;
            validate_non_negative_index("beta_m", "K", K);
            std::vector<local_scalar_t__> beta_m(K, local_scalar_t__(0));
            stan::math::initialize(beta_m, DUMMY_VAR__);
            stan::math::fill(beta_m, DUMMY_VAR__);

            current_statement_begin__ = 143;
            validate_non_negative_index("beta_a", "K", K);
            std::vector<local_scalar_t__> beta_a(K, local_scalar_t__(0));
            stan::math::initialize(beta_a, DUMMY_VAR__);
            stan::math::fill(beta_a, DUMMY_VAR__);

            // transformed parameters block statements
            current_statement_begin__ = 145;
            if (as_bool(logical_eq(trend_indicator, 0))) {

                current_statement_begin__ = 146;
                stan::math::assign(trend, linear_trend(k, m, delta, t, A, t_change, S, T, pstream__));
            } else if (as_bool(logical_eq(trend_indicator, 1))) {

                current_statement_begin__ = 148;
                stan::math::assign(trend, logistic_trend(k, m, delta, t, cap, A, t_change, S, T, pstream__));
            } else if (as_bool(logical_eq(trend_indicator, 2))) {

                current_statement_begin__ = 150;
                stan::math::assign(trend, flat_trend(m, T, pstream__));
            }
            current_statement_begin__ = 153;
            for (int i = 1; i <= K; ++i) {

                current_statement_begin__ = 154;
                stan::model::assign(beta_m, 
                            stan::model::cons_list(stan::model::index_uni(i), stan::model::nil_index_list()), 
                            (get_base1(beta, i, "beta", 1) * get_base1(s_m, i, "s_m", 1)), 
                            "assigning variable beta_m");
                current_statement_begin__ = 155;
                stan::model::assign(beta_a, 
                            (get_base1(beta, i, "beta", 1) * get_base1(s_a, i, "s_a", 1)), 
                            "assigning variable beta_a");
            }
            current_statement_begin__ = 158;
            for (int i = 1; i <= T; ++i) {

                current_statement_begin__ = 159;
                stan::model::assign(Y, 
                            ((get_base1(trend, i, "trend", 1) * (1 + dot_product(get_base1(X, i, "X", 1), beta_m))) + dot_product(get_base1(X, i, "X", 1), beta_a)), 
                            "assigning variable Y");
            }

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

            size_t trend_k_0_max__ = T;
            for (size_t k_0__ = 0; k_0__ < trend_k_0_max__; ++k_0__) {
                if (stan::math::is_uninitialized(trend[k_0__])) {
                    std::stringstream msg__;
                    msg__ << "Undefined transformed parameter: trend" << "[" << k_0__ << "]";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable trend: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            size_t Y_k_0_max__ = T;
            for (size_t k_0__ = 0; k_0__ < Y_k_0_max__; ++k_0__) {
                if (stan::math::is_uninitialized(Y[k_0__])) {
                    msg__ << "Undefined transformed parameter: Y" << "[" << k_0__ << "]";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable Y: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            size_t beta_m_k_0_max__ = K;
            for (size_t k_0__ = 0; k_0__ < beta_m_k_0_max__; ++k_0__) {
                if (stan::math::is_uninitialized(beta_m[k_0__])) {
                    msg__ << "Undefined transformed parameter: beta_m" << "[" << k_0__ << "]";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable beta_m: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }
            size_t beta_a_k_0_max__ = K;
            for (size_t k_0__ = 0; k_0__ < beta_a_k_0_max__; ++k_0__) {
                if (stan::math::is_uninitialized(beta_a[k_0__])) {
                    msg__ << "Undefined transformed parameter: beta_a" << "[" << k_0__ << "]";
                    stan::lang::rethrow_located(std::runtime_error(std::string("Error initializing variable beta_a: ") + msg__.str()), current_statement_begin__, prog_reader__());
                }
            }

            // model body

            current_statement_begin__ = 167;
            lp_accum__.add(normal_log<propto__>(k, 0, 5));
            current_statement_begin__ = 168;
            lp_accum__.add(normal_log<propto__>(m, 0, 5));
            current_statement_begin__ = 169;
            lp_accum__.add(double_exponential_log<propto__>(delta, 0, tau));
            current_statement_begin__ = 170;
            lp_accum__.add(normal_log<propto__>(sigma_obs, 0, 0.5));
            current_statement_begin__ = 171;
            lp_accum__.add(normal_log<propto__>(beta, 0, sigmas));
            current_statement_begin__ = 174;
            lp_accum__.add(normal_log<propto__>(y, Y, sigma_obs));

        }

        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());
        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("k");
        names__.push_back("m");
        names__.push_back("delta");
        names__.push_back("sigma_obs");
        names__.push_back("beta");
        names__.push_back("trend");
        names__.push_back("Y");
        names__.push_back("beta_m");
        names__.push_back("beta_a");
    }

    void get_dims(std::vector<std::vector<size_t> >& dimss__) const {
        dimss__.resize(0);
        std::vector<size_t> dims__;
        dims__.resize(0);
        dimss__.push_back(dims__);
        dims__.push_back(S);
        dims__.push_back(K);
        dims__.push_back(T);
    }

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

        vars__.resize(0);
        stan::io::reader<local_scalar_t__> in__(params_r__, params_i__);
        static const char* function__ = "anon_model_f5236004a3fd5b8429270d00efcc0cf9_namespace::write_array";

        // read-transform, write parameters
        double k = in__.scalar_constrain();
        vars__.push_back(k);

        double m = in__.scalar_constrain();
        vars__.push_back(m);

        std::vector<double> delta;
        size_t delta_d_0_max__ = S;
        delta.reserve(delta_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < delta_d_0_max__; ++d_0__) {
            delta.push_back(in__.scalar_constrain());
        }
            vars__.push_back(delta[k_0__]);
        }

        double sigma_obs = in__.scalar_lb_constrain(0);
        vars__.push_back(sigma_obs);

        std::vector<double> beta;
        size_t beta_d_0_max__ = K;
        beta.reserve(beta_d_0_max__);
        for (size_t d_0__ = 0; d_0__ < beta_d_0_max__; ++d_0__) {
            beta.push_back(in__.scalar_constrain());
        }
            vars__.push_back(beta[k_0__]);
        }

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

        if (!include_tparams__ && !include_gqs__) return;

        try {
            // declare and define transformed parameters
            std::vector<double> trend(T, double(0));

            std::vector<double> Y(T, double(0));

            std::vector<double> beta_m(K, double(0));

            std::vector<double> beta_a(K, double(0));

            // do transformed parameters statements

            }

            }

            }

            if (!include_gqs__ && !include_tparams__) return;

            // write transformed parameters
            if (include_tparams__) {
                size_t trend_k_0_max__ = T;
                for (size_t k_0__ = 0; k_0__ < trend_k_0_max__; ++k_0__) {
                    vars__.push_back(trend[k_0__]);
                }
                size_t Y_k_0_max__ = T;
                for (size_t k_0__ = 0; k_0__ < Y_k_0_max__; ++k_0__) {
                    vars__.push_back(Y[k_0__]);
                }
                size_t beta_m_k_0_max__ = K;
                for (size_t k_0__ = 0; k_0__ < beta_m_k_0_max__; ++k_0__) {
                    vars__.push_back(beta_m[k_0__]);
                }
                size_t beta_a_k_0_max__ = K;
                for (size_t k_0__ = 0; k_0__ < beta_a_k_0_max__; ++k_0__) {
                    vars__.push_back(beta_a[k_0__]);
                }
            }
            if (!include_gqs__) return;
        }
    }

    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());
        params_r_vec[i] = params_r(i);
      std::vector<double> vars_vec;
      write_array(base_rng, params_r_vec, params_i_vec, vars_vec, include_tparams, include_gqs, pstream);
      vars.resize(vars_vec.size());
      for (int i = 0; i < vars.size(); ++i)
        vars(i) = vars_vec[i];
    }

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

    void constrained_param_names(std::vector<std::string>& param_names__,
                                 bool include_tparams__ = true,
                                 bool include_gqs__ = true) const {
        std::stringstream param_name_stream__;
        param_name_stream__.str(std::string());
        param_name_stream__ << "k";
        param_names__.push_back(param_name_stream__.str());
        param_name_stream__ << "m";
            param_name_stream__.str(std::string());
            param_name_stream__ << "delta" << '.' << k_0__ + 1;
            param_names__.push_back(param_name_stream__.str());
        }
        param_name_stream__ << "sigma_obs";
            param_name_stream__ << "beta" << '.' << k_0__ + 1;
        }

        if (!include_gqs__ && !include_tparams__) return;

        if (include_tparams__) {
                param_name_stream__.str(std::string());
                param_name_stream__ << "trend" << '.' << k_0__ + 1;
                param_names__.push_back(param_name_stream__.str());
            }
                param_name_stream__ << "Y" << '.' << k_0__ + 1;
            }
                param_name_stream__ << "beta_m" << '.' << k_0__ + 1;
            }
                param_name_stream__ << "beta_a" << '.' << k_0__ + 1;
            }
        }

        if (!include_gqs__) return;
    }

    void unconstrained_param_names(std::vector<std::string>& param_names__,
                                   bool include_tparams__ = true,
                                   bool include_gqs__ = true) const {
        }
        }

            }
            }
            }
            }
        }

    }

}; // model

}  // namespace

typedef anon_model_f5236004a3fd5b8429270d00efcc0cf9_namespace::anon_model_f5236004a3fd5b8429270d00efcc0cf9 stan_model;