// generated with brms 2.14.4
functions {
  /* turn a vector into a matrix of defined dimension 
   * stores elements in row major order
   * Args: 
   *   X: a vector 
   *   N: first dimension of the desired matrix
   *   K: second dimension of the desired matrix 
   * Returns: 
   *   a matrix of dimension N x K 
   */ 
  matrix as_matrix(vector X, int N, int K) { 
    matrix[N, K] Y; 
    for (i in 1:N) {
      Y[i] = to_row_vector(X[((i - 1) * K + 1):(i * K)]); 
    }
    return Y; 
  } 
 /* compute correlated group-level effects
  * in the presence of a within-group covariance matrix
  * Args: 
  *   z: matrix of unscaled group-level effects
  *   SD: vector of standard deviation parameters
  *   L: cholesky factor correlation matrix
  *   Lcov: cholesky factor of within-group correlation matrix
  * Returns: 
  *   matrix of scaled group-level effects
  */ 
  matrix scale_r_cor_cov(matrix z, vector SD, matrix L, matrix Lcov) {
    vector[num_elements(z)] z_flat = to_vector(z);
    vector[num_elements(z)] r = rep_vector(0, num_elements(z));
    matrix[rows(L), cols(L)] LC = diag_pre_multiply(SD, L);
    int rows_z = rows(z);
    int rows_L = rows(L);
    // kronecker product of cholesky factors times a vector
    for (icov in 1:rows(Lcov)) {
      for (jcov in 1:icov) {
        if (Lcov[icov, jcov] > 1e-10) { 
          // avoid calculating products between unrelated individuals
          for (i in 1:rows_L) {
            for (j in 1:i) {
              // incremented element of the output vector
              int k = (rows_L * (icov - 1)) + i;
              // applied element of the input vector
              int l = (rows_L * (jcov - 1)) + j;
              r[k] = r[k] + Lcov[icov, jcov] * LC[i, j] * z_flat[l];
            }
          }
        }
      }
    }
    // r is returned in another dimension order than z
    return as_matrix(r, cols(z), rows(z));
  }
}
data {
  int<lower=1> N;  // total number of observations
  int<lower=1> N_bwt;  // number of observations
  vector[N_bwt] Y_bwt;  // response variable
  int<lower=1> N_tarsus;  // number of observations
  vector[N_tarsus] Y_tarsus;  // response variable
  int<lower=1> nresp;  // number of responses
  int nrescor;  // number of residual correlations
  // 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
  int<lower=1> J_1_bwt[N_bwt];  // grouping indicator per observation
  int<lower=1> J_1_tarsus[N_tarsus];  // grouping indicator per observation
  matrix[N_1, N_1] Lcov_1;  // cholesky factor of known covariance matrix
  // group-level predictor values
  vector[N_bwt] Z_1_bwt_1;
  vector[N_tarsus] Z_1_tarsus_2;
  int<lower=1> NC_1;  // number of group-level correlations
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
  vector[nresp] Y[N];  // response array
  for (n in 1:N) {
    Y[n] = transpose([Y_bwt[n], Y_tarsus[n]]);
  }
}
parameters {
  real Intercept_bwt;  // temporary intercept for centered predictors
  real<lower=0> sigma_bwt;  // residual SD
  real Intercept_tarsus;  // temporary intercept for centered predictors
  real<lower=0> sigma_tarsus;  // residual SD
  cholesky_factor_corr[nresp] Lrescor;  // parameters for multivariate linear models
  vector<lower=0>[M_1] sd_1;  // group-level standard deviations
  matrix[M_1, N_1] z_1;  // standardized group-level effects
  cholesky_factor_corr[M_1] L_1;  // cholesky factor of correlation matrix
}
transformed parameters {
  matrix[N_1, M_1] r_1;  // actual group-level effects
  // using vectors speeds up indexing in loops
  vector[N_1] r_1_bwt_1;
  vector[N_1] r_1_tarsus_2;
  // compute actual group-level effects
  r_1 = scale_r_cor_cov(z_1, sd_1, L_1, Lcov_1);
  r_1_bwt_1 = r_1[, 1];
  r_1_tarsus_2 = r_1[, 2];
}
model {
  // likelihood including all constants
  if (!prior_only) {
    // initialize linear predictor term
    vector[N_bwt] mu_bwt = Intercept_bwt + rep_vector(0.0, N_bwt);
    // initialize linear predictor term
    vector[N_tarsus] mu_tarsus = Intercept_tarsus + rep_vector(0.0, N_tarsus);
    // multivariate predictor array
    vector[nresp] Mu[N];
    vector[nresp] sigma = transpose([sigma_bwt, sigma_tarsus]);
    // cholesky factor of residual covariance matrix
    matrix[nresp, nresp] LSigma = diag_pre_multiply(sigma, Lrescor);
    for (n in 1:N_bwt) {
      // add more terms to the linear predictor
      mu_bwt[n] += r_1_bwt_1[J_1_bwt[n]] * Z_1_bwt_1[n];
    }
    for (n in 1:N_tarsus) {
      // add more terms to the linear predictor
      mu_tarsus[n] += r_1_tarsus_2[J_1_tarsus[n]] * Z_1_tarsus_2[n];
    }
    // combine univariate parameters
    for (n in 1:N) {
      Mu[n] = transpose([mu_bwt[n], mu_tarsus[n]]);
    }
    target += multi_normal_cholesky_lpdf(Y | Mu, LSigma);
  }
  // priors including all constants
  target += student_t_lpdf(Intercept_bwt | 3, 7.3, 2.8);
  target += student_t_lpdf(sigma_bwt | 3, 0, 2.8)
    - 1 * student_t_lccdf(0 | 3, 0, 2.8);
  target += student_t_lpdf(Intercept_tarsus | 3, 20.4, 5.2);
  target += student_t_lpdf(sigma_tarsus | 3, 0, 5.2)
    - 1 * student_t_lccdf(0 | 3, 0, 5.2);
  target += lkj_corr_cholesky_lpdf(Lrescor | 1);
  target += student_t_lpdf(sd_1[1] | 3, 0, 2.8)
    - 1 * student_t_lccdf(0 | 3, 0, 2.8);
  target += student_t_lpdf(sd_1[2] | 3, 0, 5.2)
    - 1 * student_t_lccdf(0 | 3, 0, 5.2);
  target += std_normal_lpdf(to_vector(z_1));
  target += lkj_corr_cholesky_lpdf(L_1 | 1);
}
generated quantities {
  // actual population-level intercept
  real b_bwt_Intercept = Intercept_bwt;
  // actual population-level intercept
  real b_tarsus_Intercept = Intercept_tarsus;
  // residual correlations
  corr_matrix[nresp] Rescor = multiply_lower_tri_self_transpose(Lrescor);
  vector<lower=-1,upper=1>[nrescor] rescor;
  // compute group-level correlations
  corr_matrix[M_1] Cor_1 = multiply_lower_tri_self_transpose(L_1);
  vector<lower=-1,upper=1>[NC_1] cor_1;
  // extract upper diagonal of correlation matrix
  for (k in 1:nresp) {
    for (j in 1:(k - 1)) {
      rescor[choose(k - 1, 2) + j] = Rescor[j, k];
    }
  }
  // extract upper diagonal of correlation matrix
  for (k in 1:M_1) {
    for (j in 1:(k - 1)) {
      cor_1[choose(k - 1, 2) + j] = Cor_1[j, k];
    }
  }
}