|
| 1 | +//! Mock “training finished” path: write a tiny `.h5ad`, patch `var` with `mean_lasso_r2` |
| 2 | +//! (same HDF5 path as real training), then read back with **anndata-rs** and optionally **h5py**. |
| 3 | +//! |
| 4 | +//! Run from repo root: |
| 5 | +//! ```text |
| 6 | +//! cargo run --example mock_h5ad_roundtrip |
| 7 | +//! ``` |
| 8 | +//! |
| 9 | +//! If `python3` + `h5py` are available, the script also probes `var/_index` (usually readable) |
| 10 | +//! and `var/mean_lasso_r2` (may fail on some arm64 `h5py` wheels when data use Blosc). |
| 11 | +
|
| 12 | +use anndata::data::ArrayData; |
| 13 | +use anndata::{AnnData, AnnDataOp, Backend}; |
| 14 | +use anndata_hdf5::H5; |
| 15 | +use ndarray::Array2; |
| 16 | +use polars::prelude::{DataFrame, NamedFrom, Series}; |
| 17 | +use spacetravlr::spatial_estimator::{MeanLassoR2Accum, dense_to_csr_f64, patch_adata_var_mean_lasso_r2}; |
| 18 | +use std::collections::HashMap; |
| 19 | +use std::path::Path; |
| 20 | +use std::process::Command; |
| 21 | +use std::sync::atomic::AtomicU64; |
| 22 | +use std::sync::Arc; |
| 23 | + |
| 24 | +fn main() -> anyhow::Result<()> { |
| 25 | + let dir = std::env::temp_dir().join(format!( |
| 26 | + "spacetravlr_mock_roundtrip_{}", |
| 27 | + std::process::id() |
| 28 | + )); |
| 29 | + let _ = std::fs::remove_dir_all(&dir); |
| 30 | + std::fs::create_dir_all(&dir)?; |
| 31 | + let p = dir.join("mock.h5ad"); |
| 32 | + |
| 33 | + println!("→ write {}", p.display()); |
| 34 | + let a = AnnData::<H5>::new(&p)?; |
| 35 | + a.set_obs_names(vec!["c0".into(), "c1".into()].into())?; |
| 36 | + a.set_var_names(vec!["G0".into(), "G1".into()].into())?; |
| 37 | + let obs = DataFrame::new(vec![Series::new( |
| 38 | + "cell_type".into(), |
| 39 | + vec!["t".to_string(), "t".to_string()], |
| 40 | + ) |
| 41 | + .into()])?; |
| 42 | + a.set_obs(obs)?; |
| 43 | + let var = DataFrame::new(vec![Series::new("gene_ids".into(), vec!["g0", "g1"]).into()])?; |
| 44 | + a.set_var(var)?; |
| 45 | + let dense = Array2::from_elem((2, 2), 0.5f64); |
| 46 | + let csr = dense_to_csr_f64(&dense)?; |
| 47 | + a.set_x(ArrayData::from(csr))?; |
| 48 | + a.close()?; |
| 49 | + |
| 50 | + println!("→ patch_adata_var_mean_lasso_r2 (adds mean_lasso_r2 column)"); |
| 51 | + let mut m: HashMap<String, usize> = HashMap::new(); |
| 52 | + m.insert("G0".into(), 0); |
| 53 | + m.insert("G1".into(), 1); |
| 54 | + let scores = Arc::new(vec![ |
| 55 | + AtomicU64::new(0.25f64.to_bits()), |
| 56 | + AtomicU64::new(f64::NAN.to_bits()), |
| 57 | + ]); |
| 58 | + let accum = MeanLassoR2Accum { |
| 59 | + gene_to_idx: Arc::new(m), |
| 60 | + scores, |
| 61 | + }; |
| 62 | + patch_adata_var_mean_lasso_r2(&p, &accum)?; |
| 63 | + |
| 64 | + println!("→ read back with anndata-rs"); |
| 65 | + let a2 = AnnData::<H5>::open(H5::open(&p)?)?; |
| 66 | + let v = a2.read_var()?; |
| 67 | + let r2 = v.column("mean_lasso_r2")?.f64()?; |
| 68 | + anyhow::ensure!((r2.get(0).unwrap() - 0.25).abs() < 1e-9); |
| 69 | + anyhow::ensure!(r2.get(1).unwrap().is_nan()); |
| 70 | + a2.close()?; |
| 71 | + println!(" OK: anndata-rs read_var + mean_lasso_r2"); |
| 72 | + |
| 73 | + try_h5py_probe(&p); |
| 74 | + |
| 75 | + let _ = std::fs::remove_dir_all(&dir); |
| 76 | + println!("→ done (temp dir removed)"); |
| 77 | + Ok(()) |
| 78 | +} |
| 79 | + |
| 80 | +fn try_h5py_probe(p: &Path) { |
| 81 | + let path = p.to_string_lossy().replace('\\', "\\\\"); |
| 82 | + let script = format!( |
| 83 | + r#" |
| 84 | +import sys |
| 85 | +path = r"{path}" |
| 86 | +try: |
| 87 | + import h5py |
| 88 | +except ImportError: |
| 89 | + print(" (skip h5py: not installed)") |
| 90 | + sys.exit(0) |
| 91 | +f = h5py.File(path, "r") |
| 92 | +vg = f["var"] |
| 93 | +names = list(vg.keys()) |
| 94 | +print(" h5py var keys:", names) |
| 95 | +for k in names: |
| 96 | + obj = vg[k] |
| 97 | + if not hasattr(obj, "shape"): |
| 98 | + continue |
| 99 | + try: |
| 100 | + sl = obj[:2] if getattr(obj, "shape", (0,))[0] != 0 else obj[()] |
| 101 | + print(f" OK: h5py read var/{{k}} (sample {{type(sl).__name__}})") |
| 102 | + except Exception as e: |
| 103 | + print(f" WARN: h5py cannot read var/{{k}}:", type(e).__name__, e) |
| 104 | +f.close() |
| 105 | +"# |
| 106 | + ); |
| 107 | + let out = Command::new("python3") |
| 108 | + .args(["-c", &script]) |
| 109 | + .output(); |
| 110 | + match out { |
| 111 | + Ok(o) => { |
| 112 | + let s = String::from_utf8_lossy(&o.stdout); |
| 113 | + let e = String::from_utf8_lossy(&o.stderr); |
| 114 | + if !s.trim().is_empty() { |
| 115 | + print!("{s}"); |
| 116 | + } |
| 117 | + if !e.trim().is_empty() && !o.status.success() { |
| 118 | + eprint!("{e}"); |
| 119 | + } |
| 120 | + } |
| 121 | + Err(e) => println!(" (skip h5py probe: {e})"), |
| 122 | + } |
| 123 | +} |
0 commit comments