[python] Use 31-bit-friendly default shape for ingest

[johnkerl]

Issue and/or context: As discussed this week and last week. The room-for-growth is good. Two billion room-for-growth in each dimension is a quite roomy default, and works better for R interop.

Tracking issue: #1445

Changes: The default domain in Python is room for plenty of growth after creation. Here we retain that -- still plenty of room for growth -- just 32-bit friendly for better ability for R code to read Python-created SOMA data.

Notes for Reviewer: This is a non-breaking change at the API level.

[codecov-commenter]

Codecov Report

Patch coverage has no change and project coverage change: +26.87 🎉

Comparison is base (5b781f1) 64.44% compared to head (9e3ab57) 91.32%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files

@@             Coverage Diff             @@
##             main    #1440       +/-   ##
===========================================
+ Coverage   64.44%   91.32%   +26.87%     
===========================================
  Files         102       30       -72     
  Lines        8336     2547     -5789     
===========================================
- Hits         5372     2326     -3046     
+ Misses       2964      221     -2743     

Flag	Coverage Δ
python	91.32% <ø> (+0.54%)	⬆️
r	?

Flags with carried forward coverage won't be shown. Click here to find out more.

see 86 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[bkmartinjr]

One bug to be resolved. The other question: in some cases, this is a breaking change, so do we want to ensure it is only released when we have the remainder of the issues resolved?

[johnkerl]

Other related issues to be proposed by @mojaveazure , as we discussed today.

I'm happy to wait on merging this PR until then, if folks prefer.

[pablo-gar]

• Requested change: The fixtures in ./apis/python/src/tiledbsoma ./apis/python/notebooks/data need to be re-generated.
• Nice-to-have: I think letting the user to specify a shape would be ideal, with the default set to 32-bit. See for reference this https://github.com/single-cell-data/TileDB-SOMA/pull/1358/files#diff-dc1f1c04136350fc3134f40bd69db3715831b77579de8a425a1eacac03effcc7

[johnkerl]

@pablo-gar two points about passing user-specified shape all the way from top level (see also #1445):

(1) I want to do that on a separate PR -- not this one

(2) You point to a diff for passing the desired shape to a single array. But tiledbsoma.io does not deal with single arrays -- it deals with all the arrays in an experiment, and we would need to offer specification for different shapes for all of them. Here for example is pbmc3k_processed:

pbmc3k_processed/obs
pbmc3k_processed/ms/RNA/var
pbmc3k_processed/ms/RNA/X/data
pbmc3k_processed/ms/RNA/obsm/X_umap
pbmc3k_processed/ms/RNA/obsm/X_tsne
pbmc3k_processed/ms/RNA/obsm/X_pca
pbmc3k_processed/ms/RNA/obsm/X_draw_graph_fr
pbmc3k_processed/ms/RNA/obsp/distances
pbmc3k_processed/ms/RNA/obsp/connectivities
pbmc3k_processed/ms/RNA/varm/PCs
pbmc3k_processed/ms/raw/X/data
pbmc3k_processed/ms/raw/var

There are 12 different arrays here, with varying shapes.

What parameterization do you propose?

tiledbsoma.io.from_anndata(
  'pbmc3k_processed',
  'pbmc3k_processed.h5ad',
  'RNA',
  shapes={
    'obs': (2638,)
    'ms/RNA/var': (1838,),
    'ms/RNA/X/data': (2638, 1838),
    'ms/RNA/obsm/X_umap': (2638, 2),
    'ms/RNA/obsm/X_tsne': (2638, 2),
    'ms/RNA/obsm/X_pca': (2638, 50),
    'ms/RNA/obsm/X_draw_graph_fr': (2638, 2),
    'ms/RNA/obsp/distances': (2638, 2638),
    'ms/RNA/obsp/connectivities': (2638, 2638),
    'ms/RNA/varm/PCs': (1838, 1838),
    'ms/raw/X/data': (2638, 13714),
    'ms/raw/var': (13714,),
  }
)

Something like the above? Or something else?

[pablo-gar]

@johnkerl Thanks! Important details on the many different arrays. As I think about it, there two main use cases we want to fulfill with this:

• Future addition of data to the tiledbsoma object -- here it makes sense to have the shapes as the maximum 32-bit capacity.
• Read-only tiledbsoma object. For this there are two options:
  • tiledbsoma object from one h5ad -- here it would make sense to have defaults as the original matrix shapes
  • tiledbsoma object from multiple h5ads (like we do in the Census) -- here it would make sense to be able to specify individual shapes for each matrix, like you have in your example.

Not sure I have the right answer since non-defaults add ux complexity.

[johnkerl]

Interop CI is failing with Error: Process completed with exit code 143 which, Google suggests to me, may be indicative of our running out of memory/CPU/time resources within the CI run ...

[johnkerl]

Interop CI is failing with Error: Process completed with exit code 143 which, Google suggests to me, may be indicative of our running out of memory/CPU/time resources within the CI run ...

It passes on my laptop with this branch checked out, so, it needs to be a CI-resource-config problem to be solved ... 🤔

[johnkerl]

Not sure I have the right answer since non-defaults add ux complexity.

@pablo-gar my position is that the above

tiledbsoma.io.from_anndata(
  'pbmc3k_processed',
  'pbmc3k_processed.h5ad',
  'RNA',
  shapes={
    'obs': (2638,)
    'ms/RNA/var': (1838,),
    'ms/RNA/X/data': (2638, 1838),
    'ms/RNA/obsm/X_umap': (2638, 2),
    'ms/RNA/obsm/X_tsne': (2638, 2),
    'ms/RNA/obsm/X_pca': (2638, 50),
    'ms/RNA/obsm/X_draw_graph_fr': (2638, 2),
    'ms/RNA/obsp/distances': (2638, 2638),
    'ms/RNA/obsp/connectivities': (2638, 2638),
    'ms/RNA/varm/PCs': (1838, 1838),
    'ms/raw/X/data': (2638, 13714),
    'ms/raw/var': (13714,),
  }
)

is prohibitively complex to require people to type. Just not OK.

Another option is we require them to do

tiledbsoma.io.from_anndata(
  'pbmc3k_processed',
  'pbmc3k_processed.h5ad',
  'RNA',
  nobs=2638, nvar, 1838, nrawvar = 13714,
)

or some such -- let the num-PCAs axes etc be tight.

The best option, I believe, is:

• If not appendable: use tight shapes
• If appendable (see also [python] Design and implement append-mode logic #1278) then use 2^31-1 defaults
• Let them specify nobs nvar nrawvar as above if they want to

[johnkerl]

Summary

Re the CI fail:

Interop CI is failing with Error: Process completed with exit code 143 which, Google suggests to me, may be indicative of our running out of memory/CPU/time resources within the CI run ...

This is due to increased runtime and a combination of several factors:

• On this PR from Python's from_h5ad/from_anndata we replace 63-bit shape with 31-bit shape
• Both cause problems with R reading that, but at different points
• Muffle-warning logic in tiledbsoma-r hides this
• With 31-bit shape, it runs longer (in particular, too long for CI) when the Python-written shapes are 31-bit

Details

• Here is a branch: https://github.com/single-cell-data/TileDB-SOMA/compare/kerl/r-alloc-trace
  • That's a branch off today's main
  • What I vary in my experiment is whether a SOMA experiment was ingested with 63-bit shape (../python/pre) or 31-bit shape (../python/post, pre-prepared on my laptop)
  • Then in my experiment I used the above branch -- today's main but with some print statements added

Output for the pre case:

ℹ Loading tiledbsoma
Warning: argument must be coercible to non-negative integer
Warning: argument must be coercible to non-negative integer
Warning: argument must be coercible to non-negative integer
Warning: argument must be coercible to non-negative integer
SEAQ SPARSE_MATRIX PRE
Warning: Array dimensions must not exceed '.Machine$integer.max'
SEAQ SPARSE_MATRIX PRE
Warning: Array dimensions must not exceed '.Machine$integer.max'
An object of class Seurat
1838 features across 2638 samples within 1 assay
Active assay: RNA (1838 features, 0 variable features)
Warning messages:
1: In as.integer.integer64(embed$shape()[2L]) :
  NAs produced by integer overflow
2: In as.integer.integer64(embed$shape()[2L]) :
  NAs produced by integer overflow
3: In as.integer.integer64(embed$shape()[2L]) :
  NAs produced by integer overflow
4: In as.integer.integer64(embed$shape()[2L]) :
  NAs produced by integer overflow

Output for the post case:

ℹ Loading tiledbsoma
Warning: vector memory exhausted (limit reached?)
Warning: vector memory exhausted (limit reached?)
Warning: vector memory exhausted (limit reached?)
Warning: vector memory exhausted (limit reached?)
SEAQ SPARSE_MATRIX PRE
SEAQ SPARSE_MATRIX POST
SEAQ MAT IDX IDX PRE
Warning: Cholmod error 'problem too large' at file ../Core/cholmod_memory.c, line 135
SEAQ SPARSE_MATRIX PRE
SEAQ SPARSE_MATRIX POST
SEAQ MAT IDX IDX PRE
Warning: Cholmod error 'problem too large' at file ../Core/cholmod_memory.c, line 135
An object of class Seurat
1838 features across 2638 samples within 1 assay
Active assay: RNA (1838 features, 0 variable features)

From the above program output we see that in the pre case:

• This line is reached https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L719
• This line is not https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L721
• The muffleWarning here means that even though line 721 isn't reached, the test fails to fail: https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L368

From the above program output we see that in the post case:

• This line is reached https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L728
• This line is not https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L739
• The muffleWarning here means that even though line 721 isn't reached, the test fails to fail: https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L368
• Yet the code runs a bit longer than in the pre case, hence the timeout
• This is why it fails to fail on my laptop -- no CI artificial resource limits -- although execution times goes up significantly, and still, we have the issue noted above of failure to reach line 739, and hence, incomplete data processing both before and after this PR

Action items

• Due to R idiosyncracies (tiledbsoma-py has zero problems reading or writing 63-bit shapes, or 31-bit shapes) we are encountering either failures in R or extreme slowness in R trying to write even 31-bit shapes written by Python
  • Therefore we need to delay this PR until @mojaveazure 's bounding-box logic on [r/python] Parent/tracking issue for shape #1445
• I don't like the muffleWarning at all -- this makes previously failing cases appear to succeed
  • I'd like to get rid of our use of muffleWarning
  • Failing that, I think we should use options(warn=2) in our unit-test cases -- as we've just seen, we had incomplete I/O even before this PR, and the R/Python interop tests weren't catching them

cc @mojaveazure @ebezzi @pablo-gar @aaronwolen

[eddelbuettel]

Can you isolate where the cholmod error comes from? That is from inside package Matrix and indicative that we masy still do something wrong.

[johnkerl]

Can you isolate where the cholmod error comes from? That is from inside package Matrix and indicative that we masy still do something wrong.

@eddelbuettel it's at this line:
https://github.com/single-cell-data/TileDB-SOMA/blob/kerl/r-alloc-trace/apis/r/R/SOMAExperimentAxisQuery.R#L729
where mat is of type dgCMatrix

[eddelbuettel]

Good. If you put a line with browser() there we can evaluate idx. Is it outside the matrix dims? For "normal" requests it behaves:

> set.seed(123); N <- 100; M <- as(rsparsematrix(N, N, 0.2), "dgCMatrix"); M[N,N]
[1] 0
> set.seed(123); N <- 100; M <- as(rsparsematrix(N, N, 0.2), "dgCMatrix"); M[N+1,N+1]
Error in intI(i, n = d[1L], dn[[1L]], give.dn = FALSE) : 
  index larger than maximal 100
> 
> for (i in 4:7) { set.seed(123); N <- 10^i; M <- as(rsparsematrix(N, N, 1000/(N*N)), "dgCMatrix"); print(M[N,N]) }
[1] 0
[1] 0
[1] 0
[1] 0
> 

[johnkerl]

@eddelbuettel https://gist.github.com/johnkerl/567bdfc1a6c750be1b7fec8792605eb3

I believe the 31-bit shape, not the idx itself, is to blame ...

Browse[1]> mat
2147483646 x 2147483646 sparse Matrix of class "dgCMatrix"
Error in subCsp_cols(x, j, drop = drop) :
  Cholmod error 'problem too large' at file ../Core/cholmod_memory.c, line 135

Browse[1]> idx
   [1]    1    2    3    4    5    6    7    8    9   10   11   12   13   14
  [15]   15   16   17   18   19   20   21   22   23   24   25   26   27   28
...
[2619] 2619 2620 2621 2622 2623 2624 2625 2626 2627 2628 2629 2630 2631 2632
[2633] 2633 2634 2635 2636 2637 2638
Browse[1]>

[johnkerl]

I just rebased on now-merged #1504

[johnkerl]

I just rebased on now-merged #1504

@aaronwolen @mojaveazure

• Before that rebase, as narrated in thorough detail above
  • The test pytest apis/system/tests/test_write_anndata_read_seurat.py was passing (with warnings) on my laptop with very long runtime
  • It was 143'ing in CI (out of resources) -- which I (mistakenly) attributed to the high runtime
• After that rebase:
  • The test pytest apis/system/tests/test_write_anndata_read_seurat.py is passing (with warnings) on my laptop with very good runtime (just a few seconds)
  • But still 143'ing in CI -- which I am now hypothesizing is due to RAM needs within the CI container

Not sure where to go next from here.

I'm not willing to merge a PR which turns interop CI (or any CI for that matter) from green to red.

Files were regenerated. Extra non-default-shape arguments are a great idea, and outside the scope of this PR.

[johnkerl]

@aaronwolen @pablo-gar after rebase on #1504 #1508 #1521 this is now a simple, green-CI PR

[aaronwolen]

🚀