Speed up voxel lookup, and reject an unordered lookup instead of trusting it

[gouwens]

Two related changes to _matching_voxel_indices, the shared voxel-matching
helper: it stops paying a cost proportional to the size of the lookup, and it
stops trusting the ordering it depends on.

---

1. Resolve membership from the binary search, not a full-column scan

The helper answered two questions with two operations:

has_match = np.isin(voxel_inds, matching_voxel_lookup[:, lookup_ind])
lookup_results_ind = np.searchsorted(matching_voxel_lookup[:, lookup_ind], ...)

np.isin walks and internally sorts the entire key column, so the call cost
O(rows) regardless of how few voxels were queried. The real
closest_surface_voxel_lookup.h5 holds 61,911,881 rows, and
angle.find_closest_streamline calls the helper with a single voxel — so one
coordinate paid for a full pass over the whole table.

np.searchsorted already returns the left-hand insertion point, so a query is
present exactly when the key sitting at that row equals it. Membership is now
read off the search the helper already performs, in O(queries · log rows).

Same matches, same missing_value, same tie-breaking (still the left insertion
point, so the stable sorter IsocortexCoordinateProjector passes for issue #12
keeps doing what it did).

call	before	after
1 voxel vs the 61.9M-row lookup	12.5 s	0.046 ms
40k mixed hit/miss queries	13.5 s	105 ms
find_closest_streamline, 10 coordinates	127 s	1.45 s
project_coordinates, 300 coordinates	14.3 s	1.17 s

Outputs are identical, old vs new, on:

• the full 61.9M-row closest surface voxel lookup (40k queries mixing present
  keys, absent keys, and both out-of-range ends);
• all nine real view lookups through the sorter path, including
  flatmap_butterfly with its 350,956 tied keys, which is where tie-breaking
  would show up if it had changed;
• angle.find_closest_streamline and
  IsocortexCoordinateProjector.project_coordinates end to end on real files.

Plus 1,600 randomized small-array comparisons against the previous
implementation (duplicate keys, out-of-range queries, both missing_values,
with and without a sorter).

---

2. Reject an unordered lookup

searchsorted requires non-decreasing keys and does not verify it, so a
reference file whose lookup column was out of order produced wrong voxel matches
with no error — the failure mode CLAUDE.md carried as a gotcha, and which
test_without_a_sorter_an_unsorted_lookup_returns_wrong_answers_silently
characterized.

_check_lookup_is_ordered now raises ValueError naming the first out-of-order
entry, so a bad file can be located rather than merely rejected:

Lookup column 0 must increase monotonically; entry 5000001 (268579732) is
smaller than entry 5000000 (268579733). Matches are found by binary search,
which returns wrong voxels for an unordered column.

It covers both orderings the helper searches — the plain column, and the column
as reordered by sorter.

Where the check runs is the design problem

The scan is O(rows): 74 ms on the real 61.9M-row table, against 0.03 ms
for the lookup itself. find_closest_streamline looks up a single voxel per
call, so checking on every call would reinstate precisely the per-coordinate
O(rows) cost part 1 removes — the check would undo the fix.

So the result is memoised per array in a module-level WeakValueDictionary.
Values are weak, so an entry dies with its array and an id cannot be reused
while its entry is live.

	after part 1	with the check
first lookup vs the 61.9M-row table	0.046 ms	74 ms (scans once)
every later lookup	0.046 ms	0.030 ms
find_closest_streamline, 10 coords	1.45 s	1.76 s
project_coordinates, 300 coords	1.17 s	~1.2 s

The sorter path is not memoised and does not need to be:
_calculate_2d_coordinates builds its sorter with np.argsort on every call, so
an O(rows) check is strictly cheaper than what is already there.

The trade-off, stated plainly

A lookup array reordered in place after its first successful use is not
re-checked. Reference files are read once and left alone, so this is a fair
trade — but test_the_ordering_check_is_remembered_per_array pins it, so the
cache cannot be dropped without someone noticing they have restored a
per-coordinate O(rows) cost. If you would rather have unconditional checking,
deleting the cache lookup is a one-line change that costs ~74 ms per coordinate
in find_closest_streamline.

---

Notes for review

• The helper always required the lookup column to be ordered — that is what
  searchsorted needs — but only the caller-facing gotcha said so. It is now
  the docstring and an enforced precondition.
• test_the_key_column_is_never_scanned_end_to_end forbids np.isin/np.in1d
  via monkeypatch as a stand-in for the O(rows) cost; a timing assertion would
  make the same point far less reliably.
• The CLAUDE.md gotcha claimed the opposite behaviour and is updated.

uv run pytest: 366 passed. Real-data tier: 36 passed.

🤖 Generated with Claude Code

[gouwens]

Follow-up: the unordered lookup now raises

Pushed 21684ad, which closes the gotcha this PR's description left open — an
out-of-order lookup column used to yield wrong voxels silently.

_check_lookup_is_ordered runs inside _matching_voxel_indices and raises
ValueError naming the first out-of-order entry, so a bad reference file can be
located, not just rejected:

Lookup column 0 must increase monotonically; entry 5000001 (268579732) is
smaller than entry 5000000 (268579733). Matches are found by binary search,
which returns wrong voxels for an unordered column.

It covers both orderings the helper searches — the plain column, and the column
as reordered by sorter.

Where the check runs is the design problem

The scan is O(rows): 74 ms on the real 61.9M-row table, against 0.03 ms
for the lookup itself. angle.find_closest_streamline looks up a single voxel
per call, so checking on every call would reinstate precisely the per-coordinate
O(rows) cost this PR removes — worse than useless.

So the result is memoised per array in a module-level WeakValueDictionary.
Values are weak, so an entry dies with its array and an id cannot be reused
while its entry is live.

	before this commit	with the check
first lookup vs the 61.9M-row table	0.046 ms	74 ms (scans once)
every later lookup	0.046 ms	0.030 ms
find_closest_streamline, 10 coords	1.45 s	1.76 s
project_coordinates, 300 coords	1.17 s	~1.2 s

The sorter path is not memoised and does not need to be: _calculate_2d_coordinates
builds its sorter with np.argsort on every call, so an O(rows) check is strictly
cheaper than what is already there.

The trade-off, stated plainly

A lookup array reordered in place after its first successful use is not
re-checked. Reference files are read once and left alone, so this is a fair
trade — but test_the_ordering_check_is_remembered_per_array pins it, so the
cache cannot be dropped without someone noticing they have restored a
per-coordinate O(rows) cost.

test_without_a_sorter_an_unsorted_lookup_returns_wrong_answers_silently became
test_without_a_sorter_an_unsorted_lookup_raises — the characterization it
documented is now a defect that errors.

uv run pytest: 366 passed. Real-data tier: 36 passed.