Cbind or subsetting greta_array with fixed zeros

hrlai · March 23, 2022, 11:04pm

Hi, not sure if this deserves an issue on GitHub, so posting here. I encountered some strange behaviour while cbind to or subsetting from greta arrays with fixed zeroes, as generated from zeros(). A minimal example to illustrate:

library(greta)

# let's say we want to create a square matrix,
# where some columns are fixed to zeroes
n_row <- 10
n_col <- 10

# Approach 1: generate non-zero and zero columns separately, then cbind
x1 <- normal(0, 1, dim = c(n_row, 2))
x2 <- zeros(n_row, n_col - 2)
X  <- cbind(x1, x2)
X   # but columns 3-10 aren't fixed at zero...

# Approach 2: generate the target matrix with zeroes, then fill in the unknowns
# this works (as used in most greta example codes)
X <- zeros(n_row, n_col)
X[, 1:2] <- x1
X

# but when subsetting...
X[, 1:2] 
X[, 1:3]  # 3rd column aren't zeros
X[, 1:4]  # likewise...

Just wondering if this is intended? The X matrix from method 2 still works for mcmc etc., just that I feel slightly unease by the subsetting output… (e.g., not sure if it’s fine if someone uses X[, 1:3] as a generated quantities for post-hoc analyses).

njtierney · March 25, 2022, 1:20am

Hi there!

Thanks for posting!

Just wanted to confirm that I get your results

Excuse some of the startup messages below in the example, I’m currently trying to get greta working on M1 and these are some of the tricks I’m trying out.

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

# let's say we want to create a square matrix,
# where some columns are fixed to zeroes
n_row <- 10
n_col <- 10

# Approach 1: generate non-zero and zero columns separately, then cbind
x1 <- normal(0, 1, dim = c(n_row, 2))
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✓ Initialising python and checking dependencies ... done!
#> 
#> eager flag value is FALSE
#> ℹ Disabling TF eager execution with: `tf$compat$v1$disable_eager_execution()`

x1
#> greta array (variable following a normal distribution)
#> 
#>       [,1] [,2]
#>  [1,]  ?    ?  
#>  [2,]  ?    ?  
#>  [3,]  ?    ?  
#>  [4,]  ?    ?  
#>  [5,]  ?    ?  
#>  [6,]  ?    ?  
#>  [7,]  ?    ?  
#>  [8,]  ?    ?  
#>  [9,]  ?    ?  
#> [10,]  ?    ?

x2 <- zeros(n_row, n_col - 2)

x2
#> greta array (data)
#> 
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#>  [1,]    0    0    0    0    0    0    0    0
#>  [2,]    0    0    0    0    0    0    0    0
#>  [3,]    0    0    0    0    0    0    0    0
#>  [4,]    0    0    0    0    0    0    0    0
#>  [5,]    0    0    0    0    0    0    0    0
#>  [6,]    0    0    0    0    0    0    0    0
#>  [7,]    0    0    0    0    0    0    0    0
#>  [8,]    0    0    0    0    0    0    0    0
#>  [9,]    0    0    0    0    0    0    0    0
#> [10,]    0    0    0    0    0    0    0    0

X  <- cbind(x1, x2)
X   # but columns 3-10 aren't fixed at zero...
#> greta array (operation)
#> 
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [2,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [3,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [4,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [5,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [6,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [7,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [8,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#>  [9,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?   
#> [10,]  ?    ?    ?    ?    ?    ?    ?    ?    ?    ?

# Approach 2: generate the target matrix with zeroes, then fill in the unknowns
# this works (as used in most greta example codes)
X <- zeros(n_row, n_col)
X[, 1:2] <- x1
X
#> greta array (operation)
#> 
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [2,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [3,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [4,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [5,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [6,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [7,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [8,]  ?    ?   0    0    0    0    0    0    0    0    
#>  [9,]  ?    ?   0    0    0    0    0    0    0    0    
#> [10,]  ?    ?   0    0    0    0    0    0    0    0

# but when subsetting...
X[, 1:2] 
#> greta array (operation)
#> 
#>       [,1] [,2]
#>  [1,]  ?    ?  
#>  [2,]  ?    ?  
#>  [3,]  ?    ?  
#>  [4,]  ?    ?  
#>  [5,]  ?    ?  
#>  [6,]  ?    ?  
#>  [7,]  ?    ?  
#>  [8,]  ?    ?  
#>  [9,]  ?    ?  
#> [10,]  ?    ?
X[, 1:3]  # 3rd column aren't zeros
#> greta array (operation)
#> 
#>       [,1] [,2] [,3]
#>  [1,]  ?    ?    ?  
#>  [2,]  ?    ?    ?  
#>  [3,]  ?    ?    ?  
#>  [4,]  ?    ?    ?  
#>  [5,]  ?    ?    ?  
#>  [6,]  ?    ?    ?  
#>  [7,]  ?    ?    ?  
#>  [8,]  ?    ?    ?  
#>  [9,]  ?    ?    ?  
#> [10,]  ?    ?    ?
X[, 1:4]  # likewise...
#> greta array (operation)
#> 
#>       [,1] [,2] [,3] [,4]
#>  [1,]  ?    ?    ?    ?  
#>  [2,]  ?    ?    ?    ?  
#>  [3,]  ?    ?    ?    ?  
#>  [4,]  ?    ?    ?    ?  
#>  [5,]  ?    ?    ?    ?  
#>  [6,]  ?    ?    ?    ?  
#>  [7,]  ?    ?    ?    ?  
#>  [8,]  ?    ?    ?    ?  
#>  [9,]  ?    ?    ?    ?  
#> [10,]  ?    ?    ?    ?

^{Created on 2022-03-25 by the reprex package (v2.0.1)}

Session info

sessioninfo::session_info()
#> ─ Session info  🇫🇯  ⚽  🇲🇾   ─────────────────────────────────────────────────
#>  hash: flag: Fiji, soccer ball, flag: Malaysia
#> 
#>  setting  value
#>  version  R version 4.1.3 (2022-03-10)
#>  os       macOS Big Sur 11.2.2
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_AU.UTF-8
#>  ctype    en_AU.UTF-8
#>  tz       Australia/Perth
#>  date     2022-03-25
#>  pandoc   2.16.1 @ /Applications/RStudio.app/Contents/MacOS/quarto/bin/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version    date (UTC) lib source
#>  backports     1.4.1      2021-12-13 [1] CRAN (R 4.1.1)
#>  base64enc     0.1-3      2015-07-28 [1] CRAN (R 4.1.0)
#>  callr         3.7.0      2021-04-20 [1] CRAN (R 4.1.0)
#>  cli           3.2.0      2022-02-14 [1] CRAN (R 4.1.1)
#>  coda          0.19-4     2020-09-30 [1] CRAN (R 4.1.0)
#>  codetools     0.2-18     2020-11-04 [1] CRAN (R 4.1.3)
#>  crayon        1.5.0      2022-02-14 [1] CRAN (R 4.1.1)
#>  digest        0.6.29     2021-12-01 [1] CRAN (R 4.1.1)
#>  ellipsis      0.3.2      2021-04-29 [1] CRAN (R 4.1.0)
#>  evaluate      0.15       2022-02-18 [1] CRAN (R 4.1.1)
#>  fansi         1.0.2      2022-01-14 [1] CRAN (R 4.1.1)
#>  fastmap       1.1.0      2021-01-25 [1] CRAN (R 4.1.0)
#>  fs            1.5.2      2021-12-08 [1] CRAN (R 4.1.1)
#>  future        1.24.0     2022-02-19 [1] CRAN (R 4.1.1)
#>  globals       0.14.0     2020-11-22 [1] CRAN (R 4.1.0)
#>  glue          1.6.2      2022-02-24 [1] CRAN (R 4.1.1)
#>  greta       * 0.4.2.9000 2022-03-24 [1] local
#>  here          1.0.1      2020-12-13 [1] CRAN (R 4.1.0)
#>  highr         0.9        2021-04-16 [1] CRAN (R 4.1.0)
#>  hms           1.1.1      2021-09-26 [1] CRAN (R 4.1.1)
#>  htmltools     0.5.2      2021-08-25 [1] CRAN (R 4.1.1)
#>  jsonlite      1.8.0      2022-02-22 [1] CRAN (R 4.1.1)
#>  knitr         1.37       2021-12-16 [1] CRAN (R 4.1.1)
#>  lattice       0.20-45    2021-09-22 [1] CRAN (R 4.1.3)
#>  lifecycle     1.0.1      2021-09-24 [1] CRAN (R 4.1.1)
#>  listenv       0.8.0      2019-12-05 [1] CRAN (R 4.1.0)
#>  magrittr      2.0.2      2022-01-26 [1] CRAN (R 4.1.1)
#>  Matrix        1.4-0      2021-12-08 [1] CRAN (R 4.1.3)
#>  parallelly    1.30.0     2021-12-17 [1] CRAN (R 4.1.1)
#>  pillar        1.7.0      2022-02-01 [1] CRAN (R 4.1.1)
#>  pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.1.0)
#>  png           0.1-7      2013-12-03 [1] CRAN (R 4.1.0)
#>  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.1.0)
#>  processx      3.5.2      2021-04-30 [1] CRAN (R 4.1.0)
#>  progress      1.2.2      2019-05-16 [1] CRAN (R 4.1.0)
#>  ps            1.6.0      2021-02-28 [1] CRAN (R 4.1.0)
#>  purrr         0.3.4      2020-04-17 [1] CRAN (R 4.1.0)
#>  R.cache       0.15.0     2021-04-30 [1] CRAN (R 4.1.0)
#>  R.methodsS3   1.8.1      2020-08-26 [1] CRAN (R 4.1.0)
#>  R.oo          1.24.0     2020-08-26 [1] CRAN (R 4.1.0)
#>  R.utils       2.11.0     2021-09-26 [1] CRAN (R 4.1.1)
#>  R6            2.5.1      2021-08-19 [1] CRAN (R 4.1.1)
#>  Rcpp          1.0.8      2022-01-13 [1] CRAN (R 4.1.1)
#>  reprex        2.0.1      2021-08-05 [1] CRAN (R 4.1.1)
#>  reticulate    1.24-9000  2022-01-28 [1] Github (rstudio/reticulate@8347417)
#>  rlang         1.0.2      2022-03-04 [1] CRAN (R 4.1.1)
#>  rmarkdown     2.11       2021-09-14 [1] CRAN (R 4.1.1)
#>  rprojroot     2.0.2      2020-11-15 [1] CRAN (R 4.1.0)
#>  rstudioapi    0.13       2020-11-12 [1] CRAN (R 4.1.0)
#>  sessioninfo   1.2.1      2021-11-02 [1] CRAN (R 4.1.1)
#>  stringi       1.7.6      2021-11-29 [1] CRAN (R 4.1.1)
#>  stringr       1.4.0      2019-02-10 [1] CRAN (R 4.1.1)
#>  styler        1.6.2      2021-09-23 [1] CRAN (R 4.1.1)
#>  tensorflow    2.7.0.9000 2022-01-28 [1] Github (rstudio/tensorflow@ff9eaeb)
#>  tfruns        1.5.0      2021-02-26 [1] CRAN (R 4.1.0)
#>  tibble        3.1.6      2021-11-07 [1] CRAN (R 4.1.1)
#>  utf8          1.2.2      2021-07-24 [1] CRAN (R 4.1.0)
#>  vctrs         0.3.8      2021-04-29 [1] CRAN (R 4.1.0)
#>  whisker       0.4        2019-08-28 [1] CRAN (R 4.1.0)
#>  withr         2.5.0      2022-03-03 [1] CRAN (R 4.1.1)
#>  xfun          0.30       2022-03-02 [1] CRAN (R 4.1.1)
#>  yaml          2.3.5      2022-02-21 [1] CRAN (R 4.1.1)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library
#> 
#> ─ Python configuration ───────────────────────────────────────────────────────
#>  python:         /Users/nick/Library/r-miniconda-arm64/envs/greta-env/bin/python
#>  libpython:      /Users/nick/Library/r-miniconda-arm64/envs/greta-env/lib/libpython3.8.dylib
#>  pythonhome:     /Users/nick/Library/r-miniconda-arm64/envs/greta-env:/Users/nick/Library/r-miniconda-arm64/envs/greta-env
#>  version:        3.8.12 | packaged by conda-forge | (default, Jan 30 2022, 23:13:24)  [Clang 11.1.0 ]
#>  numpy:          /Users/nick/Library/r-miniconda-arm64/envs/greta-env/lib/python3.8/site-packages/numpy
#>  numpy_version:  1.19.5
#>  tensorflow:     /Users/nick/Library/r-miniconda-arm64/envs/greta-env/lib/python3.8/site-packages/tensorflow
#>  
#>  NOTE: Python version was forced by use_python function
#> 
#> ──────────────────────────────────────────────────────────────────────────────

My understanding (which could be incorrect) is that the matrix gets filled with unknowns, but these will get computed as zero, or normal, in subsequent operations. Are subsequent calculations behaving as you’d expect? Or is everything is coerced to normal? Or just more so that the print method is not what you’d expect?

I agree it is interesting that you can get the zeros to appear when doing replacement in place (X[, 1:2] <- x1).

Thanks again for posting!

hrlai · March 25, 2022, 2:05am

Hmmm I was too afraid to sample the posterior using Approach 1, so I didn’t do it. But hearing you said that it could simply be a print issue, let’s sample the prior:

library(greta)

x <- normal(0, 1, dim = 2)
y <- zeros(1)
z <- c(x, z)  # 3rd row element prints "?"

# sample from prior to check z[3, 1]
sim <- calculate(z, nsim = 1000)

# print(sim) and you'll see that the z[3, 1] samples are all zeroes
# but let's plot it
par(mfrow = c(3, 1), mar = c(4, 4, 0, 0))
plot(sim$z[, 1, 1], type = "l")
plot(sim$z[, 2, 1], type = "l")
plot(sim$z[, 3, 1], type = "l")

It seems that you were right! It’s just the print method not displaying what I expected. No big deal, but it’d be great if this could be one of the future enhancement

njtierney · March 28, 2022, 3:01am

OK, awesome! Glad to hear it is just a printing issue and there isn’t a horrible hidden bug there.

To help discuss this new printing feature/enhancement in greta, if you’re able to post it as an issue on the greta github then I can allocate it to a milestone to make sure it gets done!

njtierney · March 29, 2022, 12:51am

I’ve just added this issue here: https://github.com/greta-dev/greta/issues/512