Validate get_subset_from_range against coordinates, not variant count - #178
Conversation
low and high are genomic coordinates. The guard compared high to positions.size, which is the number of variants, so it rejected any range whose end coordinate was larger than the variant count. That is nearly every real range: a matrix of 100 variants over 1 Mb refused the window [200000, 300000). WindowedAnalyzer.compute_region is the only caller and failed the same way. The guard also hid a second fault in the other direction. A range that passed it but held no variant reached the constructor, which rejects zero-size arrays, and raised "genotypes cannot be empty". get_subset already had an empty-result path; this method had none. Keep the low < 0 and low >= high checks and drop the third. A range that reaches past the last variant is a valid query and returns the variants that do fall inside it. Two more faults in the same method: * chrom_end was set to high, which is exclusive. chrom_end is inclusive: the constructor counts accessible bases up to chrom_end + 1, from_ts stores sequence_length - 1, and the bp window iterator writes end - 1. The subset claimed one base more than it covered. * The sample names were dropped, so load_pop_file did not work on the result. Move the empty-result builder out of get_subset into _empty_subset, which both methods now call. It also sets fields, which the old code left unset -- reading .fields on an empty subset raised AttributeError. Fold the three copies of the accessible-base count into _accessible_bases_in_range. An empty region matrix now reaches the window iterator, which read positions[0] without a guard. Return no windows instead.
|
I was about to add that the subset was dropping the acc mask, but it seems this PR is fixing that as well. |
|
@mufernando -- would you review this? |
mufernando
left a comment
mufernando
left a comment
There was a problem hiding this comment.
it looks good to me. I tested on my use case and it works now.
I just left a question about whether we need to warn the user that the high is greater than the accessible mask they already have. the current code does the right thing and clips the high anyway.
| # only bounds that mean anything are that the range starts at or after | ||
| # zero and is not empty. A range that reaches past the last variant is | ||
| # a valid query: it returns the variants that do fall inside it. | ||
| if low < 0 or low >= high: |
There was a problem hiding this comment.
should we check high is not greater than the size of the accessible mask array?
There was a problem hiding this comment.
the current code does the right thing if that's the case, because accessiblemask.slice get's the min out of the high or len of the array.
There was a problem hiding this comment.
good question -- i don't think we need the extra check. slice and count_accessible both clip, and bases past the mask end are just inaccessible, which is the answer we want
|
i'm gonna merge this |
2 participants
Fixes #172.
lowandhighare genomic coordinates. The guard comparedhightopositions.size,which is the variant count:
Everything else in the method treats them as coordinates: it filters
self.positions,calls
accessible_mask.slice(low, high), and stores them aschrom_start/chrom_end.The sibling
get_subsetvalidates againsthaplotypes.shape[1]correctly, because itsargument really is indices. This check looks inherited from there. It has been in the
file since the first commit.
What it broke
Almost every real range was refused. Only ranges inside
[0, n_variants)passed. Amatrix of 100 variants over 1 Mb rejected the window
[200000, 300000)withValueError: Invalid range specified.WindowedAnalyzer.compute_region(
windowed_analysis.py:631) is the only caller and failed the same way.A range with no variants crashed. A range that passed the guard but held no variant
reached the constructor, which rejects zero-size arrays, and raised
ValueError: genotypes cannot be empty.get_subsethas an empty-result path; thismethod had none, so the guard was hiding the gap.
The fix
Keep
low < 0andlow >= high, and drop the third check. A range that reaches pastthe last variant is a valid query -- it returns the variants that do fall inside it.
Validating against
chrom_start/chrom_endwould recreate the same problem for windowsat the edge of a chromosome.
Two further faults in the same method:
chrom_endwas set tohigh, which is exclusive.chrom_endis inclusive: theconstructor counts accessible bases up to
chrom_end + 1,from_tsstoressequence_length - 1with a comment saying so, and_iter_bp_windowswritesend - 1. The subset claimed one base more than it covered.sampleslist was dropped, soload_pop_filedid not work on the result.Two cleanups that the fix needs:
_empty_subsetnow holds the empty-result builder that onlyget_subsethad. It alsosets
fields, which the old code left unset -- reading.fieldson an empty subsetraised
AttributeError. That object is reachable today:test_windowed_analysis.pybuilds one, and every variant-free window produces one.
_accessible_bases_in_rangereplaces three copies of the same accessible-base count.Outside the method
An empty region matrix now reaches
WindowIterator, which readpositions_np[0]with noguard.
_iter_bp_windowsandcount_windowsreturn nothing for a matrix with novariants. Without this,
compute_regionwould trade one error for another. Drop thispart if you would rather keep the PR inside
haplotype_matrix.py.Tests
get_subset_from_rangehad two tests and both were written around the bug.test_accessible_mask.py:370says so outright: "get_subset_from_range checkshigh <= positions.size, so use a range within the number of positions".
test_haplotype_matrix.py:221passes only because that fixture's 10th position is belowits variant count. Both still pass unchanged.
New tests cover a coordinate range far past the variant count, the inclusive
chrom_end,the kept sample names, an empty window, a range past the last variant, malformed ranges,
and
fieldson an empty subset.compute_regiongets its first two tests.Not fixed here
_iter_bp_windowsanchors windows atpositions[0]and ignoreschrom_start. Afterthis fix,
compute_region(hm, "1", 200000, 300000)returns windows that start at thefirst variant in the region and can end past 300000, so two neighbouring region calls
give windows that do not line up. That is a decision about what
compute_regionpromises, so it belongs in its own issue.