// generated with brms 2.19.0
functions {
  /* compute a latent Gaussian process
   * Args:
   *   x: array of continuous predictor values
   *   sdgp: marginal SD parameter
   *   lscale: length-scale parameter
   *   zgp: vector of independent standard normal variables
   * Returns:
   *   a vector to be added to the linear predictor
   */
  vector gp(data array[] vector x, real sdgp, vector lscale, vector zgp) {
    int Dls = rows(lscale);
    int N = size(x);
    matrix[N, N] cov;
    if (Dls == 1) {
      // one dimensional or isotropic GP
      cov = gp_exp_quad_cov(x, sdgp, lscale[1]);
    } else {
      // multi-dimensional non-isotropic GP
      cov = gp_exp_quad_cov(x[ : , 1], sdgp, lscale[1]);
      for (d in 2 : Dls) {
        cov = cov .* gp_exp_quad_cov(x[ : , d], 1, lscale[d]);
      }
    }
    for (n in 1 : N) {
      // deal with numerical non-positive-definiteness
      cov[n, n] += 1e-12;
    }
    return cholesky_decompose(cov) * zgp;
  }
  /* Spectral density function of a Gaussian process
   * with squared exponential covariance kernel
   * Args:
   *   x: array of numeric values of dimension NB x D
   *   sdgp: marginal SD parameter
   *   lscale: vector of length-scale parameters
   * Returns:
   *   numeric values of the function evaluated at 'x'
   */
  vector spd_cov_exp_quad(data array[] vector x, real sdgp, vector lscale) {
    int NB = dims(x)[1];
    int D = dims(x)[2];
    int Dls = rows(lscale);
    vector[NB] out;
    if (Dls == 1) {
      // one dimensional or isotropic GP
      real constant = square(sdgp) * (sqrt(2 * pi()) * lscale[1]) ^ D;
      real neg_half_lscale2 = -0.5 * square(lscale[1]);
      for (m in 1 : NB) {
        out[m] = constant * exp(neg_half_lscale2 * dot_self(x[m]));
      }
    } else {
      // multi-dimensional non-isotropic GP
      real constant = square(sdgp) * sqrt(2 * pi()) ^ D * prod(lscale);
      vector[Dls] neg_half_lscale2 = -0.5 * square(lscale);
      for (m in 1 : NB) {
        out[m] = constant * exp(dot_product(neg_half_lscale2, square(x[m])));
      }
    }
    return out;
  }
  /* how many elements are in a range of integers?
   * Args:
   *   x: an integer array
   *   start: start of the range (inclusive)
   *   end: end of the range (inclusive)
   * Returns:
   *   a scalar integer
   */
  int size_range(array[] int x, int start, int end) {
    int out = 0;
    for (i in 1 : size(x)) {
      out += x[i] >= start && x[i] <= end;
    }
    return out;
  }
  /* which elements are in a range of integers?
   * Args:
   *   x: an integer array
   *   start: start of the range (inclusive)
   *   end: end of the range (inclusive)
   * Returns:
   *   an integer array
   */
  array[] int which_range(array[] int x, int start, int end) {
    array[size_range(x, start, end)] int out;
    int j = 1;
    for (i in 1 : size(x)) {
      if (x[i] >= start && x[i] <= end) {
        out[j] = i;
        j += 1;
      }
    }
    return out;
  }
  /* adjust array values to x - start + 1
   * Args:
   *   x: an integer array
   *   start: start of the range of values in x (inclusive)
   * Returns:
   *   an integer array
   */
  array[] int start_at_one(array[] int x, int start) {
    array[size(x)] int out;
    for (i in 1 : size(x)) {
      out[i] = x[i] - start + 1;
    }
    return out;
  }
}
data {
  int<lower=1> N; // total number of observations
  vector[N] Y; // response variable
  // data related to GPs
  int<lower=1> Kgp_1; // number of sub-GPs (equal to 1 unless 'by' was used)
  int<lower=1> Dgp_1; // GP dimension
  // number of latent GP groups
  int<lower=1> Nsubgp_1;
  // indices of latent GP groups per observation
  array[N] int<lower=1> Jgp_1;
  array[Nsubgp_1] vector[Dgp_1] Xgp_1; // covariates of the GP
  // data for group-level effects of ID 1
  int<lower=1> N_1; // number of grouping levels
  int<lower=1> M_1; // number of coefficients per level
  array[N] int<lower=1> J_1; // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_1_1;
  // data for group-level effects of ID 2
  int<lower=1> N_2; // number of grouping levels
  int<lower=1> M_2; // number of coefficients per level
  array[N] int<lower=1> J_2; // grouping indicator per observation
  // group-level predictor values
  vector[N] Z_2_1;
  int prior_only; // should the likelihood be ignored?
}
transformed data {
  
}
parameters {
  real Intercept; // temporary intercept for centered predictors
  vector<lower=0>[Kgp_1] sdgp_1; // GP standard deviation parameters
  array[Kgp_1] vector<lower=0>[1] lscale_1; // GP length-scale parameters
  vector[Nsubgp_1] zgp_1; // latent variables of the GP
  real<lower=0> sigma; // dispersion parameter
  vector<lower=0>[M_1] sd_1; // group-level standard deviations
  array[M_1] vector[N_1] z_1; // standardized group-level effects
  vector<lower=0>[M_2] sd_2; // group-level standard deviations
  array[M_2] vector[N_2] z_2; // standardized group-level effects
}
transformed parameters {
  vector[N_1] r_1_1; // actual group-level effects
  vector[N_2] r_2_1; // actual group-level effects
  real lprior = 0; // prior contributions to the log posterior
  r_1_1 = sd_1[1] * z_1[1];
  r_2_1 = sd_2[1] * z_2[1];
  lprior += student_t_lpdf(Intercept | 3, 4.6, 2.5);
  lprior += student_t_lpdf(sdgp_1 | 3, 0, 2.5)
            - 1 * student_t_lccdf(0 | 3, 0, 2.5);
  lprior += inv_gamma_lpdf(lscale_1[1][1] | 2.86035, 0.388975);
  lprior += student_t_lpdf(sigma | 3, 0, 2.5)
            - 1 * student_t_lccdf(0 | 3, 0, 2.5);
  lprior += student_t_lpdf(sd_1 | 3, 0, 2.5)
            - 1 * student_t_lccdf(0 | 3, 0, 2.5);
  lprior += student_t_lpdf(sd_2 | 3, 0, 2.5)
            - 1 * student_t_lccdf(0 | 3, 0, 2.5);
}
model {
  // likelihood including constants
  if (!prior_only) {
    vector[Nsubgp_1] gp_pred_1 = gp(Xgp_1, sdgp_1[1], lscale_1[1], zgp_1);
    // initialize linear predictor term
    vector[N] mu = rep_vector(0.0, N);
    mu += Intercept + gp_pred_1[Jgp_1];
    for (n in 1 : N) {
      // add more terms to the linear predictor
      mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
    }
    target += lognormal_lpdf(Y | mu, sigma);
  }
  // priors including constants
  target += lprior;
  target += std_normal_lpdf(zgp_1);
  target += std_normal_lpdf(z_1[1]);
  target += std_normal_lpdf(z_2[1]);
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
}