`tapply` cannot produce 2D table from greta arrays

Hi, this is another case showing that `tapply` is not yet fully compatible with greta arrays. Currently `tapply` works well when the `INDEX` argument is a vector of IDs, but not when we supply a list of vectors to produce >1D output. Would it be possible to get this working in an upcoming release? Not sure if this is a low hanging fruit or not. I'm personally using this to generate a model matrix / covariates from weighted sums. Thanks!

``` r
library(greta)
#>
#> Attaching package: 'greta'
#> The following objects are masked from 'package:stats':
#>
#> binomial, cov2cor, poisson
#> The following objects are masked from 'package:base':
#>
#> %*%, apply, backsolve, beta, chol2inv, colMeans, colSums, diag,
#> eigen, forwardsolve, gamma, identity, rowMeans, rowSums, sweep,
#> tapply

set.seed(123)

n <- 100

row <- sample(1:5, n, replace = TRUE)
col <- sample(1:2, n, replace = TRUE)

x <- rnorm(n, 0, 1)

beta <- normal(0, 1)
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!

beta_x <- beta * x

tapply(x, list(row, col), sum) # works
#> 1 2
#> 1 -4.3243935 -3.01086731
#> 2 -2.5329235 0.08524451
#> 3 -0.8380712 1.61988790
#> 4 3.1535623 -2.85590318
#> 5 1.2308146 -0.20006750

tapply(beta_x, row, "sum") # works
#> greta array
#> [,1]
#> [1,] ?
#> [2,] ?
#> [3,] ?
#> [4,] ?
#> [5,] ?
#>

tapply(beta_x, list(row, col), "sum") # doesn't work
#> Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...): 'x' must be atomic
```

^{Created on 2025-03-26 with [reprex v2.1.1](https://reprex.tidyverse.org)}