Ordered matrix along both rows and columns

Modeling
3

That issue gives Stan code only for scalar bounds but an ordered matrix needs the bounds to be vectors. As discussed in the issue, vector bounds are much more difficult to implement.
Fortunately, the really difficult case happens when both upper and lower bounds are vectors–with only a lower bound a smooth constraining transform is possible.

functions {
  matrix ordered_matrix_lp(vector spine, matrix submatrix) {
    int N = rows(submatrix) + 1;
    int M = cols(submatrix) + 1;
    matrix[N,M] output;
    output[1,:] = spine[1:M]';
    for (r in 2:N) {
      output[r,M] = spine[M-1 + r];
      for (c in 1:M-1) {
        int j = M-c;
        real diff = output[r,j+1] - output[r-1,j];
        output[r,j] = output[r-1,j] + diff*submatrix[r-1,j];
        target += log(diff);
      }
    }
    return output;
  }
}
data {
  int<lower=1> N;
  int<lower=1> M;
}
parameters {
  ordered[N+M-1] spine;
  matrix<lower=0,upper=1>[N-1,M-1] submatrix;
}
transformed parameters {
  matrix[N,M] mat = ordered_matrix_lp(spine, submatrix);
}