11#!/usr/bin/env python3
2- """Reproducible-features demo: train a model on `rosalind features` output and
2+ """Reproducible-features demo: train a model on native Arrow feature output and
33prove the inputs (and therefore the trained model) are bit-reproducible.
44
55The headline is REPRODUCIBILITY, not biological novelty: identical inputs ->
6- byte-identical features -> bit-identical trained weights, verifiable by the
7- BLAKE3 receipt. Runs with numpy only (no pandas/ sklearn/pyarrow needed).
6+ byte-identical Arrow streams -> bit-identical trained weights, verifiable by the
7+ BLAKE3 receipt. Runs with NumPy and PyArrow (no pandas or sklearn needed).
88
99 python examples/reproducible_features_demo.py [path-to-rosalind-binary]
1010"""
2020import tempfile
2121
2222import numpy as np
23+ import pyarrow as pa
2324
2425REPO = os .path .dirname (os .path .dirname (os .path .abspath (__file__ )))
25- sys .path .insert (0 , os .path .join (REPO , "python" ))
26- from rosalind import features # noqa: E402
26+
27+ NUMERIC_FEATURES = (
28+ "depth" , "raw_depth" , "a" , "c" , "g" , "t" ,
29+ "a_fwd" , "a_rev" , "c_fwd" , "c_rev" ,
30+ "g_fwd" , "g_rev" , "t_fwd" , "t_rev" ,
31+ "mean_bq" , "mean_mapq" ,
32+ )
2733
2834
2935def find_binary (argv : list [str ]) -> str :
@@ -44,12 +50,37 @@ def sh(cmd: list[str]) -> None:
4450
4551
4652def output_blake3 (manifest_path : str ) -> str :
47- """Read the recorded BLAKE3 of the feature TSV from the run receipt."""
53+ """Read the recorded BLAKE3 of the Arrow stream from the run receipt."""
4854 with open (manifest_path ) as fh :
4955 m = json .load (fh )
5056 return m ["outputs" ][0 ]["blake3" ]
5157
5258
59+ def extract_features (binary : str , reference : str , bam : str , directory : str ):
60+ """Write one native Arrow artifact and return its table and evidence paths."""
61+ os .makedirs (directory , exist_ok = True )
62+ artifact = os .path .join (directory , "features.arrow" )
63+ sh ([
64+ binary , "features" , "--reference-pack" , reference ,
65+ "--alignments" , bam , "--format" , "arrow-ipc" , "--output" , artifact ,
66+ ])
67+ with pa .memory_map (artifact , "r" ) as source :
68+ table = pa .ipc .open_stream (source ).read_all ()
69+ return table , artifact , f"{ artifact } .manifest.json"
70+
71+
72+ def model_arrays (table : pa .Table ):
73+ """Convert only the explicit training columns after bounded extraction."""
74+ columns = [
75+ table [name ].to_numpy (zero_copy_only = False ).astype (np .float64 , copy = False )
76+ for name in NUMERIC_FEATURES
77+ ]
78+ values = np .column_stack (columns )
79+ positions = table ["pos" ].to_numpy (zero_copy_only = False )
80+ references = np .asarray (table ["ref" ].to_pylist (), dtype = "U1" )
81+ return values , positions , references
82+
83+
5384def train_logreg (x : np .ndarray , y : np .ndarray , iters : int = 300 , lr : float = 0.2 ):
5485 """Deterministic pure-numpy logistic regression (zero init, fixed schedule)."""
5586 mu = x .mean (0 )
@@ -78,54 +109,61 @@ def main() -> int:
78109 work = tempfile .mkdtemp (prefix = "rosalind-demo-" )
79110 data = os .path .join (work , "toy" )
80111
81- print ("== building a toy dataset + index ==" )
112+ print ("== building a toy dataset + analysis reference pack ==" )
82113 sh ([sys .executable , os .path .join (REPO , "scripts" , "generate_toy_data.py" ), data ])
83114 ref = os .path .join (data , "reference.fa" )
84115 reads = os .path .join (data , "reads_R1.fastq" )
85116 raw = os .path .join (work , "raw.bam" )
86117 bam = os .path .join (work , "sorted.bam" )
87- idx = os .path .join (work , "ref.idx " )
118+ reference_pack = os .path .join (work , "ref.rref " )
88119 sh ([binary , "align" , "--reference" , ref , "--reads" , reads , "--format" , "bam" , "--output" , raw ])
89120 sh ([binary , "sort" , "--input" , raw , "--output" , bam ])
90- sh ([binary , "index" , "--reference" , ref , "--output" , idx ])
91-
92- print ("== extracting features TWICE (independent runs) ==" )
93- a = features (idx , bam , binary = binary , workdir = os .path .join (work , "a" ))
94- b = features (idx , bam , binary = binary , workdir = os .path .join (work , "b" ))
121+ sh ([binary , "reference" , "build" , "--fasta" , ref , "--output" , reference_pack ])
122+
123+ print ("== extracting native Arrow features TWICE (independent runs) ==" )
124+ a , arrow_a , manifest_a = extract_features (
125+ binary , reference_pack , bam , os .path .join (work , "a" )
126+ )
127+ b , arrow_b , manifest_b = extract_features (
128+ binary , reference_pack , bam , os .path .join (work , "b" )
129+ )
130+ data_a , pos_a , ref_a = model_arrays (a )
131+ data_b , pos_b , ref_b = model_arrays (b )
95132
96133 # --- Reproducibility proof -------------------------------------------------
97- tsv_a = os . path . join ( work , "a" , "features.tsv" )
98- tsv_b = os . path . join ( work , "b" , "features.tsv" )
99- bytes_identical = open ( tsv_a , "rb" ). read () == open ( tsv_b , "rb" ). read ()
100- h_a , h_b = output_blake3 (a . manifest_path ), output_blake3 (b . manifest_path )
101- local_hash = hashlib .sha256 (open ( tsv_a , "rb" ). read () ).hexdigest ()[:16 ]
134+ with open ( arrow_a , "rb" ) as first , open ( arrow_b , "rb" ) as second :
135+ bytes_a , bytes_b = first . read (), second . read ( )
136+ bytes_identical = bytes_a == bytes_b
137+ h_a , h_b = output_blake3 (manifest_a ), output_blake3 (manifest_b )
138+ local_hash = hashlib .sha256 (bytes_a ).hexdigest ()[:16 ]
102139
103- print (f" rows: { len ( a ) } ; numeric features: { a . data .shape [1 ]} " )
140+ print (f" rows: { a . num_rows } ; numeric features: { data_a .shape [1 ]} " )
104141 print (f" features byte-identical across runs: { bytes_identical } " )
105142 print (f" receipt BLAKE3 (run A): { h_a } " )
106143 print (f" receipt BLAKE3 (run B): { h_b } " )
107144 print (f" receipts match: { h_a == h_b } " )
108- assert bytes_identical , "feature TSVs differ across runs"
145+ assert bytes_identical , "native Arrow streams differ across runs"
109146 assert h_a == h_b , "receipt hashes differ across runs"
110147
111148 # --- A genuine supervised task: is the reference base a purine (A/G)? -------
112149 # Read counts peak at the ref base, so this is learnable from the features.
113- def label (ft ) -> np .ndarray :
114- return np .isin (ft . ref , [b "A" , b "G" ]).astype (np .float64 )
150+ def label (references : np . ndarray ) -> np .ndarray :
151+ return np .isin (references , ["A" , "G" ]).astype (np .float64 )
115152
116153 # Deterministic train/test split by position parity.
117- def split (ft ):
118- test = (ft . pos % 2 ) == 1
154+ def split (positions : np . ndarray ):
155+ test = (positions % 2 ) == 1
119156 return ~ test , test
120157
121- ya , yb = label (a ), label (b )
122- tr_a , te_a = split (a )
123- w_a , mu_a , sd_a = train_logreg (a . data [tr_a ], ya [tr_a ])
124- acc_a = float ((predict (w_a , mu_a , sd_a , a . data [te_a ]) == (ya [te_a ] >= 0.5 )).mean ())
158+ ya , yb = label (ref_a ), label (ref_b )
159+ tr_a , te_a = split (pos_a )
160+ w_a , mu_a , sd_a = train_logreg (data_a [tr_a ], ya [tr_a ])
161+ acc_a = float ((predict (w_a , mu_a , sd_a , data_a [te_a ]) == (ya [te_a ] >= 0.5 )).mean ())
125162
126- w_b , mu_b , sd_b = train_logreg (b .data [split (b )[0 ]], yb [split (b )[0 ]])
127- pred_a = predict (w_a , mu_a , sd_a , a .data [te_a ])
128- pred_b = predict (w_b , mu_b , sd_b , b .data [split (b )[1 ]])
163+ tr_b , te_b = split (pos_b )
164+ w_b , mu_b , sd_b = train_logreg (data_b [tr_b ], yb [tr_b ])
165+ pred_a = predict (w_a , mu_a , sd_a , data_a [te_a ])
166+ pred_b = predict (w_b , mu_b , sd_b , data_b [te_b ])
129167
130168 weights_identical = np .array_equal (w_a , w_b )
131169 preds_identical = np .array_equal (pred_a , pred_b )
0 commit comments