pairwise indexing arrays of varying effects

I'm not sure that this belongs here vs. in the math library---a fix will need to cross both.

Here's an example where it'd be nice to be able to vectorize:

data {
  int<lower=0> N, J, K;
  array[N] int<lower=1, upper=J> location;
  array[N] int<lower=1, upper=K> time;
  array[N] int<lower=0, upper=1> is_male;
  array[N] real y;
}
parameters {
  real alpha;
  sum_to_zero_vector[J] beta;
  sum_to_zero_vector[K] gamma;
  sum_to_zero_matrix[J, K] delta;
  real<lower=0> sigma;
}
model {
  alpha ~ normal(0, 5);
  beta ~ std_normal();
  gamma ~ std_normal();
  to_vector(delta) ~ std_normal();

  for (n in 1:N) {
    y[n] ~ normal(alpha + beta[location[n]] + gamma[time[n]]
                  + delta[location[n], time[n]], sigma);
  }
}

There's currently no way to write that loop as a one-liner. What I'd like to have this behavior:

zip_index[delta, location, time][n] == delta[location[n], time[n]]

It also naturally extends to three or more indexes. The loop then vectorizes as

  y ~ norma(alpha + beta[location] + gamma[time] + zip_index[delta, location, time], sigma);

I don't find this very transparent from a reading perspective.

The other solution I could imagine is zipping the location and time explicitly into an N x 2 array and defining the indexing that way.

zip(location, time)[n] == (location[n], time[n])

This requires a new allocation of an N x 2 array, though that can be done in transformed data or required in data. It also requires extending matrix indexing to support tuples, specifically defining

delta[(j, k)] == delta[j, k]

Then the loop simplifies to

y ~ normal(alpha + beta[location] + gamma[time] + delta[zip(location, time)], sigma);

This feels much more readable to me because the delta parallels the other variables.

Why it doesn't work already

If we have matrix[M, N] x and array[J] int a, and array[K] int b, then x[a, b] is a J x K matrix defined by x[a, b][m, n] = x[a[m], b[n]].

No files or tests are named. Start by locating the Stan indexing and vectorization implementation and the boundary with the math library, then compare the zip_index and tuple-indexing proposals. Done means one approach has agreed semantics for varying index arrays, including higher-dimensional cases, with compiler and math-library coverage.