iterate_matrix transposes transition matrix but docs don't reflect this

[jdyen]

The examples use a matrix with columns moving to rows (matrix(A) vec(x)) but tf_iterate_matrix calculates rows to columns (rowvec(x) matrix(A)): https://github.com/greta-dev/greta.dynamics/blob/aa09ee2250ddb08981dbffe98943904f972611d9/R/iterate_matrix.R#L239

Not sure which format you prefer??

[jdyen]

Closed: see #7

[jdyen]

Re-opened because the multisite version is now wrong. With a k x k matrix, tf$matmul needs transpose_a = FALSE, but the same matrix expanded in a n x k x k array needs transpose_a = TRUE.

I prefer this matrix structure (k x k, columns move to rows) but don't like the idea of having to change matrix formats for the multisite version. I will try making arrays as k x k x n and see if that works.

[goldingn]

should be fixable. Do you have some code that we could use as a test for this?

[jdyen]

test script that demonstrates issue is here: https://gist.github.com/jdyen/6d7983c2d050b0547840618516b41f35

[jdyen]

the tf$matmul function is correct; the issue is in the greta arrays themselves. Can see this if you run the following code after the gist linked above:

# pull out values of single matrix
tmat_tmp <- calculate(tmat, list(survival = survival_fixed,
                                    stasis = stasis_fixed,
                                  recruit = recruit_fixed))
# pull out values of array
tmats_tmp <- calculate(tmats, list(survival = survival_fixed,
                                    stasis = stasis_fixed,
                                   recruit = recruit_fixed)) 
all.equal(tmat_tmp, tmats_tmp[1, , ])

[goldingn]

oh right, it's probably just the dim argument to op() that needs changing then. That determines the dimensions of the greta array that's created

[jdyen]

is it the dim argument or is it how greta recreates a 2D matrix from a slice of a 3D array?

[goldingn]

Is this still an issue? IF so, is it just the transpose_a argument here that needs to be turned off?

[jdyen]

I'll dive back in and have a look.

[jdyen]

I think this is still an issue:

# current
mat <- matrix(c(0.1, 0.5, 0.0, 0.0, 0.6, 0.4, 1.5, 0.0, 0.8), nrow = 3)
init <- rep(1, 3)
iterates <- iterate_matrix(mat,
                           initial_state = init,
                           niter = 1,
                           tol = -10)
calculate(iterates$all_states)

# transposed
iterates <- iterate_matrix(t(mat),
                           initial_state = init,
                           niter = 1,
                           tol = -10)
calculate(iterates$all_states)

# what it should be
mat %*% init

I would just turn off the transpose where you pointed but could alternatively update docs to match this. Depends which format you'd prefer.

[jdyen]

The r tests here treat this the same way because states %*% matrix is the same as t(matrix) %*% states, just with a transposed output.

[goldingn]

Yeah, that makes sense. So the new way is more standard, right?

[jdyen]

This is what I learnt (from Caswell's book):

People with a math/eng/economics background may well do something different. The example for ?iterate_matrix is structured in the Caswell way (transitions on lower diag, recruits in top row).

[goldingn]

Right, yes this makes sense!

So will just need to change that transpose argument, change the R version in the tests, let a couple of people know the behaviour has changed, and maybe be more explicit about this in the docs.

[jdyen]

I can make these changes and open a PR if you'd like.

[goldingn]

Yes please!