// Saturated covariance
functions {
  vector cmp_probs(real alpha, real pa1, real pa2, vector thrRaw) {
    int nth = num_elements(thrRaw);
    vector[nth] thr = cumulative_sum(thrRaw);
    vector[1+nth*2] unsummed;
    real paDiff = pa2 - pa1;
    unsummed[1] = 0;
    for (tx in 1:nth) {
      real t1 = thr[nth+1-tx];
      unsummed[1+tx] = paDiff + t1;
      unsummed[1+2*nth+1-tx] = paDiff - t1;
    }
    return cumulative_sum(alpha * unsummed);
  }
}
data {
  // dimensions
  int<lower=1> NPA;             // physical activites
  int<lower=1> NFACETS;         // facets (i.e. characteristics)
  int<lower=1> NCMP;            // comparisons (sample size)
  int<lower=1> NTHRESH;         // number of thresholds
  int<lower=1> N;               // number of non-missing data
  // actual data follows
  int<lower=1, upper=NPA> pa1[NCMP];        // PA1 for observation N
  int<lower=1, upper=NPA> pa2[NCMP];        // PA2 for observation N
}
transformed data {
  int numPars = NPA*NFACETS + NFACETS*NTHRESH + NFACETS + NFACETS*NFACETS;
  matrix[NPA,NFACETS]     rawTheta_;
  matrix[NFACETS, NTHRESH] threshold_;
  vector<lower=0>[NFACETS] alpha_;
  cholesky_factor_corr[NFACETS] rawThetaCorChol_;
  matrix[NPA,NFACETS]     theta_;    // latent score of PA by facet
  int rcat[NCMP,NFACETS];

  rawThetaCorChol_ = lkj_corr_cholesky_rng(NFACETS, 2.0);
  for (pa in 1:NPA) {
    for (fx in 1:NFACETS) {
      rawTheta_[pa,fx] = normal_rng(0,1);
    }
  }
  for (fx in 1:NFACETS) {
    threshold_[fx,1] = normal_rng(0.4, 0.1);
    for (tx in 2:NTHRESH) {
      threshold_[fx,tx] = normal_rng(0.8, 0.1);
    }
    alpha_[fx] = lognormal_rng(log(1)-0.15^2/2.0, 0.15);
    rawTheta_[,fx] /= sd(rawTheta_[,fx]);
  }
  //print(threshold_);
  //print(sigma_);

  for (pa in 1:NPA) {
    theta_[pa,] = (rawThetaCorChol_ * rawTheta_[pa,]')';
  }
  //print(theta_);

  // generate data
  for (cmp in 1:NCMP) {
    for (ff in 1:NFACETS) {
      rcat[cmp,ff] = categorical_logit_rng(cmp_probs(alpha_[ff], theta_[pa1[cmp],ff],
                                                     theta_[pa2[cmp],ff], threshold_[ff,]'));
    }
  }
  //print(rcat);
}
parameters {
  matrix[NPA,NFACETS]     rawTheta;
  matrix[NFACETS,NTHRESH] threshold;
  vector<lower=0>[NFACETS] alpha;
  cholesky_factor_corr[NFACETS] rawThetaCorChol;
}
transformed parameters {
  // non-centered parameterization due to thin data
  matrix[NPA,NFACETS]     theta;    // latent score of PA by facet
  for (pa in 1:NPA) {
    theta[pa,] = (rawThetaCorChol * rawTheta[pa,]')';
  }
}
model {
  rawThetaCorChol ~ lkj_corr_cholesky(2.0);
  for (pa in 1:NPA) {
    rawTheta[pa,] ~ normal(0,1);
  }
  for (tx in 1:NTHRESH) {
    threshold[,tx] ~ normal(0,2);
  }
  alpha ~ lognormal(log(2)-.5^2/2.0, .5);
  for (cmp in 1:NCMP) {
    for (ff in 1:NFACETS) {
      if (rcat[cmp,ff] == 13) continue;  // special value to indicate missing
      rcat[cmp,ff] ~ categorical_logit(cmp_probs(alpha[ff], theta[pa1[cmp],ff],
                                                 theta[pa2[cmp],ff], threshold[ff,]'));
    }
  }
}
generated quantities {
  vector[numPars] pars_;
  int ranks_[numPars];
  int px=1;

  for (pa in 1:NPA) {
    for (fx in 1:NFACETS) {
      pars_[px] = rawTheta_[pa,fx];
      ranks_[px] = rawTheta[pa,fx] > rawTheta_[pa,fx];;
      px += 1;
    }
  }
  for (fx in 1:NFACETS) {
    for (tx in 1:NTHRESH) {
      pars_[px] = threshold_[fx,tx];
      ranks_[px] = threshold[fx,tx] > threshold_[fx,tx];
      px += 1;
    }
  }
  for (fx in 1:NFACETS) {
    pars_[px] = alpha_[fx];
    ranks_[px] = alpha[fx] > alpha_[fx];
    px += 1;
  }
  for (c in 1:(NFACETS)) {
    for (r in 1:NFACETS) {
      pars_[px] = rawThetaCorChol_[r,c];
      ranks_[px] = rawThetaCorChol[r,c] > rawThetaCorChol_[r,c];
      px += 1;
    }
  }
}