From 48d0695190b5d72aecf4c8d5efadfeb15f5a0f1a Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Mon, 29 Jun 2026 21:18:55 -0400 Subject: [PATCH 1/7] feat: support native xarray reindex/align/concat over ManifestArrays MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Teach ManifestArray to handle xarray's reindex/alignment indexer at the array-protocol layer, so plain xr.reindex, xr.align, and xr.concat work over virtual arrays and stay lazy — missing positions become null-path manifest entries that read back as the array's fill_value, with no materialization and no dtype promotion. xarray's Variable._getitem_with_mask gathers with an integer indexer (-1 for missing) then masks with where(). The hooks that make this work lazily: - manifests/indexing.py: an integer-array indexer is routed to a chunk-grid remap (chunk_map_from_indexer -> ManifestArray._reindex_axis) instead of being rejected as fancy indexing. Each chunk-block must be all-missing (-> null chunk) or a chunk-aligned run (-> real chunk), else NotImplementedError; no chunk splitting. - manifests/reindex.py (new): chunk_map_from_indexer partitions the indexer by chunk and asserts the all-or-nothing condition. - array_api.py: np.where returns the gathered array unchanged when the requested fill mask matches its null chunks (the manifest's own fill_value governs); general masking still raises. np.result_type treats a lone ManifestArray's dtype as authoritative, so int arrays keep their dtype + declared sentinel fill rather than being promoted to float+NaN (also fixes a latent generator-consumption bug in the previous implementation). - array.py: __array_function__ now admits ndarray among the argument types (np.where's boolean condition arrives as one); adds ManifestArray._reindex_axis. Co-Authored-By: Claude Opus 4.8 (1M context) --- virtualizarr/manifests/array.py | 72 ++++++- virtualizarr/manifests/array_api.py | 58 +++++- virtualizarr/manifests/indexing.py | 70 +++++++ virtualizarr/manifests/reindex.py | 78 +++++++ .../tests/test_manifests/test_reindex.py | 132 ++++++++++++ virtualizarr/tests/test_native_align.py | 193 ++++++++++++++++++ 6 files changed, 598 insertions(+), 5 deletions(-) create mode 100644 virtualizarr/manifests/reindex.py create mode 100644 virtualizarr/tests/test_manifests/test_reindex.py create mode 100644 virtualizarr/tests/test_native_align.py diff --git a/virtualizarr/manifests/array.py b/virtualizarr/manifests/array.py index d7909b1d..2f519693 100644 --- a/virtualizarr/manifests/array.py +++ b/virtualizarr/manifests/array.py @@ -144,8 +144,11 @@ def __array_function__(self, func, types, args, kwargs) -> Any: return NotImplemented # Note: this allows subclasses that don't override - # __array_function__ to handle ManifestArray objects - if not all(issubclass(t, ManifestArray) for t in types): + # __array_function__ to handle ManifestArray objects. Plain ndarrays are + # also permitted among the argument types because some handled functions + # legitimately mix them with ManifestArrays — e.g. np.where, whose boolean + # condition arrives as an ndarray during xarray's reindex/alignment fill. + if not all(issubclass(t, (ManifestArray, np.ndarray)) for t in types): return NotImplemented return MANIFESTARRAY_HANDLED_ARRAY_FUNCTIONS[func](*args, **kwargs) @@ -333,6 +336,71 @@ def with_fill_value_only(self, fill_value: Any) -> "ManifestArray": ) return ManifestArray(metadata=new_metadata, chunkmanifest=empty_manifest) + def _reindex_axis( + self, axis: int, chunk_map: list[int | None], new_size: int + ) -> "ManifestArray": + """ + Return a new ManifestArray with the chunk grid along ``axis`` remapped. + + Each entry of ``chunk_map`` is the source chunk index to copy into that + target chunk slot, or ``None`` for a missing (null-path) chunk that reads + back as ``fill_value``. ``new_size`` is the new length of ``axis`` in + elements; the chunk size is unchanged. + """ + from virtualizarr.manifests.manifest import ( + MISSING_CHUNK_PATH, + ChunkManifest, + ) + from virtualizarr.manifests.utils import copy_and_replace_metadata + + manifest = self.manifest + src_paths = manifest._paths + src_offsets = manifest._offsets + src_lengths = manifest._lengths + + new_grid_shape = list(src_paths.shape) + new_grid_shape[axis] = len(chunk_map) + + new_paths = np.full( + new_grid_shape, MISSING_CHUNK_PATH, dtype=np.dtypes.StringDType() + ) + new_offsets = np.zeros(new_grid_shape, dtype=np.uint64) + new_lengths = np.zeros(new_grid_shape, dtype=np.uint64) + + new_inlined: dict[tuple[int, ...], bytes] = {} + for new_idx, src_chunk in enumerate(chunk_map): + if src_chunk is None: + continue # leave this slab as missing/null + src_slice: list[Any] = [slice(None)] * src_paths.ndim + src_slice[axis] = src_chunk + dst_slice: list[Any] = [slice(None)] * src_paths.ndim + dst_slice[axis] = new_idx + new_paths[tuple(dst_slice)] = src_paths[tuple(src_slice)] + new_offsets[tuple(dst_slice)] = src_offsets[tuple(src_slice)] + new_lengths[tuple(dst_slice)] = src_lengths[tuple(src_slice)] + # re-key any inlined chunks that lived in this source slab + for key, data in manifest._inlined.items(): + if key[axis] == src_chunk: + shifted = list(key) + shifted[axis] = new_idx + new_inlined[tuple(shifted)] = data + + new_manifest = ChunkManifest.from_arrays( + paths=new_paths, + offsets=new_offsets, + lengths=new_lengths, + validate_paths=False, + inlined=new_inlined if new_inlined else None, + ) + + new_shape = list(self.shape) + new_shape[axis] = new_size + new_metadata = copy_and_replace_metadata( + old_metadata=self.metadata, new_shape=new_shape + ) + + return ManifestArray(chunkmanifest=new_manifest, metadata=new_metadata) + def to_virtual_variable(self) -> xr.Variable: """ Create a "virtual" xarray.Variable containing the contents of one zarr array. diff --git a/virtualizarr/manifests/array_api.py b/virtualizarr/manifests/array_api.py index 164e6147..55e255d0 100644 --- a/virtualizarr/manifests/array_api.py +++ b/virtualizarr/manifests/array_api.py @@ -36,13 +36,28 @@ def decorator(func): @implements(np.result_type) def result_type(*arrays_and_dtypes: Union["ManifestArray", np.dtype]) -> np.dtype: - """Called by xarray to ensure all arguments to concat have the same dtype.""" + """ + Resolve a result dtype for ManifestArray arguments. + + Called by xarray both to check that concat/stack inputs share a dtype and, + during reindex/alignment, to combine a ManifestArray with a scalar fill + value. A ManifestArray's dtype is fixed by its metadata (the same metadata + that defines its fill value), so when exactly one ManifestArray is combined + with scalars/dtypes we return its dtype rather than promoting — keeping the + array's native dtype and declared fill value through a reindex. + """ from virtualizarr.manifests.array import ManifestArray - dtypes = ( + manifest_dtypes = [ + obj.dtype for obj in arrays_and_dtypes if isinstance(obj, ManifestArray) + ] + if len(manifest_dtypes) == 1: + return manifest_dtypes[0] + + dtypes = [ obj.dtype if isinstance(obj, ManifestArray) else np.dtype(obj) for obj in arrays_and_dtypes - ) + ] first_dtype, *other_dtypes = dtypes unique_dtypes = set(dtypes) for other_dtype in other_dtypes: @@ -54,6 +69,43 @@ def result_type(*arrays_and_dtypes: Union["ManifestArray", np.dtype]) -> np.dtyp return first_dtype +@implements(np.where) +def where(condition, x, y, /): + """ + Support xarray's reindex/alignment fill, which calls + ``where(~mask, gathered_array, fill_value)`` after gathering chunks. + + The gathered ManifestArray already carries null-path chunks (which read back + as ``fill_value``) at exactly the missing positions, so this is an identity + whenever the requested fill positions coincide with the array's missing + chunks. In that case ``x`` is returned unchanged and the manifest's own fill + value governs. Any other ``where`` usage (e.g. general boolean masking) would + require materializing values and is not supported. + """ + from virtualizarr.manifests.array import ManifestArray + + if isinstance(x, ManifestArray) and np.isscalar(y): + cond = np.asarray(condition, dtype=bool) + if cond.shape == x.shape and np.array_equal(~cond, _missing_element_mask(x)): + return x + + raise NotImplementedError( + "np.where on a ManifestArray is only supported for the reindex/alignment " + "fill pattern (filling an array's own missing chunks); general masking " + "would require materializing values." + ) + + +def _missing_element_mask(marr: "ManifestArray") -> np.ndarray: + """Boolean element-mask (shape == marr.shape), True at missing (null) chunks.""" + from virtualizarr.manifests.manifest import MISSING_CHUNK_PATH + + mask = marr.manifest._paths == MISSING_CHUNK_PATH + for axis, chunk_size in enumerate(marr.metadata.chunks): + mask = np.repeat(mask, chunk_size, axis=axis) + return mask[tuple(slice(0, length) for length in marr.shape)] + + @implements(np.concatenate) def concatenate( arrays: tuple["ManifestArray", ...] | list["ManifestArray"], diff --git a/virtualizarr/manifests/indexing.py b/virtualizarr/manifests/indexing.py index 18dc18db..ebe2e914 100644 --- a/virtualizarr/manifests/indexing.py +++ b/virtualizarr/manifests/indexing.py @@ -145,6 +145,68 @@ def apply_indexer( return output_arr +def _is_reindex_indexer(indexer_1d: T_BasicIndexer_1d) -> bool: + """True if this 1D indexer is an integer array (xarray's reindex indexer).""" + return isinstance(indexer_1d, np.ndarray) and indexer_1d.dtype.kind in ("i", "u") + + +def _is_full_slice(indexer_1d: T_BasicIndexer_1d, length: int) -> bool: + """True if this slice selects the whole axis (a no-op).""" + if not isinstance(indexer_1d, slice): + return False + return ( + indexer_1d == slice(None) + or indexer_1d == slice(0, length, None) + or (indexer_1d.indices(length) == (0, length, 1)) + ) + + +def _apply_reindex( + marr: "ManifestArray", indexers: tuple[T_BasicIndexer_1d, ...] +) -> "ManifestArray": + """ + Remap the chunk grid for one or more integer-array (reindex) indexers. + + Each integer-array axis is reindexed by translating its indexer into a + chunk-grid map (null chunks for missing positions). Any other axis must be a + full-axis slice (the shape xarray's reindex/alignment produces); combining + reindexing with chunk subsetting in a single call is not supported. + """ + from virtualizarr.manifests.reindex import chunk_map_from_indexer + + result = marr + for axis, indexer_1d in enumerate(indexers): + if _is_reindex_indexer(indexer_1d): + chunk_size = result.metadata.chunks[axis] + try: + chunk_map = chunk_map_from_indexer( + np.asarray(indexer_1d), chunk_size, result.shape[axis] + ) + except NotImplementedError as err: + # chunk_map_from_indexer is reached from deep inside xarray's + # reindex/alignment machinery, so re-raise with call-site context + # naming the operation, the axis, and what the user can do. + raise NotImplementedError( + f"Cannot align/reindex this virtual array along axis {axis} " + f"(chunk size {chunk_size}) without materializing data: the " + "target coordinate labels do not line up with chunk boundaries, " + "so filling them would require splitting a chunk. VirtualiZarr " + "only supports whole-chunk appends, inserts, and reorders. " + "Re-chunk along this dimension, or load the variable into memory, " + "if you need this alignment. See " + "https://github.com/zarr-developers/VirtualiZarr/issues/51." + ) from err + result = result._reindex_axis(axis, chunk_map, new_size=len(indexer_1d)) + elif _is_full_slice(indexer_1d, marr.shape[axis]): + continue + else: + raise NotImplementedError( + "Combining reindexing (an integer-array indexer) with chunk " + "subsetting in a single operation is not supported." + ) + return result + + def apply_selection( marr: "ManifestArray", indexer_without_newaxes: tuple[T_BasicIndexer_1d, ...] ) -> "ManifestArray": @@ -163,6 +225,14 @@ def apply_selection( # at this point there should be no ellipsis, no Nones, and one 1D indexer for each axis. assert len(indexer_without_newaxes) == marr.ndim + # An integer array indexer is xarray's reindex/alignment indexer (with -1 + # marking missing labels). Rather than reject it as fancy indexing, remap the + # chunk grid: missing positions become null-path chunks that read back as + # fill_value. This is what lets xarray's reindex/align/concat machinery work + # over ManifestArrays without materializing data. + if any(_is_reindex_indexer(ind) for ind in indexer_without_newaxes): + return _apply_reindex(marr, indexer_without_newaxes) + # validate types and reject anything we can never support per-axis narrowed_indexers: list[int | slice] = [] for indexer_1d in indexer_without_newaxes: diff --git a/virtualizarr/manifests/reindex.py b/virtualizarr/manifests/reindex.py new file mode 100644 index 00000000..74f015a3 --- /dev/null +++ b/virtualizarr/manifests/reindex.py @@ -0,0 +1,78 @@ +"""Translate xarray's reindex/alignment indexer into a chunk-grid map.""" + +from __future__ import annotations + +import numpy as np + +_SPLIT_MSG = ( + "Cannot reindex/align lazily: the requested labels would require splitting " + "or sub-chunk reordering of a source chunk, which VirtualiZarr does not do " + "(it would require reading chunk bytes). Only whole-chunk appends, inserts, " + "and reorders are supported. See https://github.com/zarr-developers/VirtualiZarr/issues/51." +) + + +def chunk_map_from_indexer( + indexer: np.ndarray, + chunk_size: int, + source_len: int, +) -> list[int | None]: + """ + Partition a reindex/alignment indexer into a chunk-grid map. + + xarray's reindex hands the backing array an integer indexer along the + reindexed axis, with ``-1`` marking positions absent from the source (the + would-be-fill positions). This partitions that indexer by chunk and returns, + for each target chunk slot, either the index of the source chunk to copy into + it, or ``None`` for an all-missing (null-path) chunk that reads back as + ``fill_value``. + + Each chunk-sized block of the indexer must be either entirely ``-1`` (→ null + chunk) or a contiguous, ascending, chunk-aligned run of source positions + (→ that source chunk). Anything else — a block mixing present and missing + positions, a sub-chunk reorder, or an unaligned start — raises, because it + cannot be expressed without splitting a chunk. + + Parameters + ---------- + indexer + Integer array (xarray's positional indexer), ``-1`` where missing. + chunk_size + Chunk size along this axis (number of elements per chunk). + source_len + Length of the source axis in elements (to size the trailing partial chunk). + + Returns + ------- + list of (int or None) + One entry per target chunk slot. + + Raises + ------ + NotImplementedError + If the indexer cannot be expressed at chunk granularity. + """ + pos = np.asarray(indexer) + chunk_map: list[int | None] = [] + for start in range(0, len(pos), chunk_size): + block = pos[start : start + chunk_size] + + if np.all(block == -1): + chunk_map.append(None) + continue + if np.any(block == -1): + raise NotImplementedError(_SPLIT_MSG) + + s = int(block[0]) + if not np.array_equal(block, np.arange(s, s + len(block))): + raise NotImplementedError(_SPLIT_MSG) + if s % chunk_size != 0: + raise NotImplementedError(_SPLIT_MSG) + src_chunk = s // chunk_size + src_chunk_size = min(chunk_size, source_len - src_chunk * chunk_size) + if len(block) != src_chunk_size: + raise NotImplementedError(_SPLIT_MSG) + + chunk_map.append(src_chunk) + + return chunk_map diff --git a/virtualizarr/tests/test_manifests/test_reindex.py b/virtualizarr/tests/test_manifests/test_reindex.py new file mode 100644 index 00000000..7e0aa90f --- /dev/null +++ b/virtualizarr/tests/test_manifests/test_reindex.py @@ -0,0 +1,132 @@ +import numpy as np +import pytest + +from virtualizarr.manifests import ChunkManifest, ManifestArray +from virtualizarr.manifests.reindex import chunk_map_from_indexer +from virtualizarr.manifests.utils import create_v3_array_metadata + + +def _idx(*vals): + return np.array(vals, dtype="int64") + + +def _marr(shape, chunks, entries, dims): + metadata = create_v3_array_metadata( + shape=shape, + chunk_shape=chunks, + data_type=np.dtype("float32"), + codecs=[{"configuration": {"endian": "little"}, "name": "bytes"}], + fill_value=float("nan"), + dimension_names=dims, + ) + return ManifestArray( + metadata=metadata, chunkmanifest=ChunkManifest(entries=entries) + ) + + +class TestChunkMapFromIndexer: + def test_append_chunk_size_1(self): + assert chunk_map_from_indexer(_idx(0, 1, 2, -1, -1), 1, 3) == [ + 0, + 1, + 2, + None, + None, + ] + + def test_chunked_append(self): + # source chunks [0,1],[2,3]; one new all-missing chunk appended + assert chunk_map_from_indexer(_idx(0, 1, 2, 3, -1, -1), 2, 4) == [0, 1, None] + + def test_prepend(self): + assert chunk_map_from_indexer(_idx(-1, -1, 0, 1), 1, 2) == [None, None, 0, 1] + + def test_gap_fill_insert_whole_chunk(self): + # [0,1] [missing] [2,3] + assert chunk_map_from_indexer(_idx(0, 1, -1, -1, 2, 3), 2, 4) == [0, None, 1] + + def test_whole_chunk_reverse_c1(self): + assert chunk_map_from_indexer(_idx(2, 1, 0), 1, 3) == [2, 1, 0] + + def test_whole_chunk_reverse_c2(self): + assert chunk_map_from_indexer(_idx(2, 3, 0, 1), 2, 4) == [1, 0] + + def test_trailing_partial_source_chunk(self): + # source len 5, C=2 -> chunk sizes 2,2,1; selecting just the size-1 tail + assert chunk_map_from_indexer(_idx(4), 2, 5) == [2] + + def test_raise_mixed_present_and_missing(self): + with pytest.raises(NotImplementedError, match="split"): + chunk_map_from_indexer(_idx(0, 1, 2, -1), 2, 4) + + def test_raise_sub_chunk_reorder(self): + with pytest.raises(NotImplementedError, match="split"): + chunk_map_from_indexer(_idx(1, 0, 2, 3), 2, 4) + + def test_raise_unaligned_start(self): + with pytest.raises(NotImplementedError, match="split"): + chunk_map_from_indexer(_idx(1, 2), 2, 4) + + def test_raise_partial_target_of_full_source(self): + # taking only half of source chunk [2,3] would split it + with pytest.raises(NotImplementedError, match="split"): + chunk_map_from_indexer(_idx(0, 1, 2), 2, 4) + + +class TestReindexAxis: + """ManifestArray._reindex_axis: the chunk-grid remap the indexer hook drives.""" + + def test_pad_with_null_chunk(self): + marr = _marr( + (4,), + (2,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 8}, + "1": {"path": "/b.nc", "offset": 8, "length": 8}, + }, + ["x"], + ) + result = marr._reindex_axis(axis=0, chunk_map=[0, None, 1], new_size=6) + + assert result.shape == (6,) + assert result.metadata.chunks == (2,) + assert result.manifest.dict() == { + "0": {"path": "file:///a.nc", "offset": 0, "length": 8}, + "2": {"path": "file:///b.nc", "offset": 8, "length": 8}, + } + + def test_whole_chunk_reorder(self): + marr = _marr( + (4,), + (2,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 8}, + "1": {"path": "/b.nc", "offset": 8, "length": 8}, + }, + ["x"], + ) + result = marr._reindex_axis(axis=0, chunk_map=[1, 0], new_size=4) + + assert result.manifest.dict() == { + "0": {"path": "file:///b.nc", "offset": 8, "length": 8}, + "1": {"path": "file:///a.nc", "offset": 0, "length": 8}, + } + + def test_only_target_axis_remapped_2d(self): + marr = _marr( + (2, 4), + (2, 2), + { + "0.0": {"path": "/a.nc", "offset": 0, "length": 16}, + "0.1": {"path": "/b.nc", "offset": 16, "length": 16}, + }, + ["y", "x"], + ) + result = marr._reindex_axis(axis=1, chunk_map=[0, 1, None], new_size=6) + + assert result.shape == (2, 6) + assert result.metadata.chunks == (2, 2) + assert result.manifest.dict() == { + "0.0": {"path": "file:///a.nc", "offset": 0, "length": 16}, + "0.1": {"path": "file:///b.nc", "offset": 16, "length": 16}, + } diff --git a/virtualizarr/tests/test_native_align.py b/virtualizarr/tests/test_native_align.py new file mode 100644 index 00000000..79ed7d5f --- /dev/null +++ b/virtualizarr/tests/test_native_align.py @@ -0,0 +1,193 @@ +""" +End-to-end tests that native xarray alignment machinery works over ManifestArrays. + +The point of the indexer/where/result_type hooks is that users call plain +`xr.reindex`, `xr.align`, and `xr.concat` — no VirtualiZarr-specific API — and the +fill stays virtual (null-path chunks that read back as the array's fill_value). +""" + +import numpy as np +import pytest +import xarray as xr +from obspec_utils.registry import ObjectStoreRegistry + +from virtualizarr.manifests import ( + ChunkManifest, + ManifestArray, + ManifestGroup, + ManifestStore, +) +from virtualizarr.manifests.utils import create_v3_array_metadata + +CODECS = [{"configuration": {"endian": "little"}, "name": "bytes"}] +F32 = { + 1.0: b"\x00\x00\x80\x3f", + 2.0: b"\x00\x00\x00\x40", + 3.0: b"\x00\x00\x40\x40", +} +I32 = {10: b"\x0a\x00\x00\x00", 20: b"\x14\x00\x00\x00", 30: b"\x1e\x00\x00\x00"} + + +def _ds_1d(x, path, dtype="float32", fill_value=float("nan"), chunk=1, name="foo"): + n = len(x) + n_chunks = -(-n // chunk) + entries = { + str(i): {"path": path, "offset": i * chunk * 4, "length": chunk * 4} + for i in range(n_chunks) + } + metadata = create_v3_array_metadata( + shape=(n,), + chunk_shape=(chunk,), + data_type=np.dtype(dtype), + codecs=CODECS, + fill_value=fill_value, + dimension_names=["x"], + ) + marr = ManifestArray( + metadata=metadata, chunkmanifest=ChunkManifest(entries=entries) + ) + return xr.Dataset({name: xr.Variable(["x"], marr)}, coords={"x": list(x)}) + + +def _ds_2d(x, t, path, name="foo"): + nt, nx = len(t), len(x) + entries = {} + k = 0 + for it in range(nt): + for ix in range(nx): + entries[f"{it}.{ix}"] = {"path": path, "offset": k * 4, "length": 4} + k += 1 + metadata = create_v3_array_metadata( + shape=(nt, nx), + chunk_shape=(1, 1), + data_type=np.dtype("float32"), + codecs=CODECS, + fill_value=float("nan"), + dimension_names=["time", "x"], + ) + marr = ManifestArray( + metadata=metadata, chunkmanifest=ChunkManifest(entries=entries) + ) + return xr.Dataset( + {name: xr.Variable(["time", "x"], marr)}, coords={"x": list(x), "time": list(t)} + ) + + +class TestNativeReindex: + def test_reindex_stays_lazy(self): + ds = _ds_1d([0, 1, 2], "/a.nc") + result = ds.reindex(x=[0, 1, 2, 3, 4]) # plain xarray, no vz API + + assert isinstance(result["foo"].data, ManifestArray) + assert list(result.x.values) == [0, 1, 2, 3, 4] + assert result["foo"].data.manifest.dict() == { + "0": {"path": "file:///a.nc", "offset": 0, "length": 4}, + "1": {"path": "file:///a.nc", "offset": 4, "length": 4}, + "2": {"path": "file:///a.nc", "offset": 8, "length": 4}, + } + + def test_reindex_float_reads_back_nan(self): + import obstore as obs + + store = obs.store.MemoryStore() + obs.put(store, "a.bin", F32[1.0] + F32[2.0] + F32[3.0]) + ds = _ds_1d([0, 1, 2], "memory:///a.bin") + + reindexed = ds.reindex(x=[0, 1, 2, 3, 4]) + + ms = ManifestStore( + group=ManifestGroup(arrays={"foo": reindexed["foo"].data}), + registry=ObjectStoreRegistry({"memory://": store}), + ) + vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values + np.testing.assert_array_equal(vals[:3], [1.0, 2.0, 3.0]) + assert np.isnan(vals[3:]).all() + + # xarray internally builds a NaN fill scalar and casts it toward the int + # dtype (a value we never use, since our np.where returns the array as-is), + # which emits a harmless "invalid value encountered in cast" RuntimeWarning. + @pytest.mark.filterwarnings( + "ignore:invalid value encountered in cast:RuntimeWarning" + ) + def test_reindex_int_keeps_dtype_and_sentinel(self): + # xarray would normally promote int->float+NaN; the manifest keeps int32 + # and its declared sentinel fill instead. + import obstore as obs + + store = obs.store.MemoryStore() + obs.put(store, "a.bin", I32[10] + I32[20] + I32[30]) + ds = _ds_1d([0, 1, 2], "memory:///a.bin", dtype="int32", fill_value=-9999) + + reindexed = ds.reindex(x=[0, 1, 2, 3, 4]) + assert reindexed["foo"].dtype == np.dtype("int32") + assert isinstance(reindexed["foo"].data, ManifestArray) + + ms = ManifestStore( + group=ManifestGroup(arrays={"foo": reindexed["foo"].data}), + registry=ObjectStoreRegistry({"memory://": store}), + ) + vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values + np.testing.assert_array_equal(vals, [10, 20, 30, -9999, -9999]) + assert vals.dtype == np.dtype("int32") + + def test_reindex_non_chunk_aligned_raises(self): + ds = _ds_1d([0, 1, 2], "/a.nc", chunk=2) + # the deep error is wrapped with call-site context: the operation, the + # axis, and what to do — and the low-level reason is chained. + with pytest.raises( + NotImplementedError, match="align/reindex this virtual array along axis 0" + ) as excinfo: + ds.reindex(x=[0, 1, 99, 2]) + assert "chunk boundaries" in str(excinfo.value) + assert isinstance(excinfo.value.__cause__, NotImplementedError) + + +class TestNativeAlign: + def test_align_outer_union(self): + a = _ds_1d([0, 1, 2], "/a.nc") + b = _ds_1d([3, 4, 5], "/b.nc") + + ra, rb = xr.align(a, b, join="outer") + + assert isinstance(ra["foo"].data, ManifestArray) + assert isinstance(rb["foo"].data, ManifestArray) + assert list(ra.x.values) == [0, 1, 2, 3, 4, 5] + + def test_align_exclude_then_concat_spatial(self): + # the ITS_LIVE pattern: align on the spatial dim (exclude time), concat on time + import obstore as obs + + store = obs.store.MemoryStore() + obs.put(store, "a.bin", F32[1.0] + F32[2.0]) # time=0, x=0,1 + obs.put(store, "b.bin", F32[3.0] + F32[1.0]) # time=1, x=2,3 + a = _ds_2d([0, 1], [0], "memory:///a.bin") + b = _ds_2d([2, 3], [1], "memory:///b.bin") + + ra, rb = xr.align(a, b, join="outer", exclude=["time"]) + # time is left alone; x is unioned + assert list(ra.time.values) == [0] + assert list(rb.time.values) == [1] + assert list(ra.x.values) == [0, 1, 2, 3] + + cube = xr.concat([ra, rb], dim="time", join="exact") + assert cube.sizes == {"time": 2, "x": 4} + assert isinstance(cube["foo"].data, ManifestArray) + + ms = ManifestStore( + group=ManifestGroup(arrays={"foo": cube["foo"].data}), + registry=ObjectStoreRegistry({"memory://": store}), + ) + vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values + np.testing.assert_array_equal(vals[0, :2], [1.0, 2.0]) + assert np.isnan(vals[0, 2:]).all() + assert np.isnan(vals[1, :2]).all() + np.testing.assert_array_equal(vals[1, 2:], [3.0, 1.0]) + + +class TestGeneralWhereStillRejected: + def test_boolean_masking_raises_clearly(self): + # we must NOT silently no-op general where(); it should raise, since it + # would require materializing values. + ds = _ds_1d([0, 1, 2], "/a.nc") + with pytest.raises(NotImplementedError, match="materializing|where"): + ds["foo"].where(ds["x"] > 0).compute() From 180aee13d98c152d630a19e05f9ab4dcae79e595 Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Tue, 30 Jun 2026 11:11:07 -0400 Subject: [PATCH 2/7] fix(manifests): support reindexing multiple axes at once MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A 2D+ outer-join reindexes more than one dimension in a single call, and xarray expresses that as a broadcast (vectorized) indexer: one integer array per axis, shaped to the array's rank with the real labels on its own axis and length 1 everywhere else (e.g. (N,1) and (1,M)). _apply_reindex assumed a 1D indexer per axis, so it misread the array's rank — raising "axis 0 (chunk size 1)" or, on newer numpy, TypeError from int() on a non-scalar block. Add _collapse_outer_indexer: since reindex indexing is orthogonal, each broadcast component reshapes losslessly to its 1D per-axis indexer. A genuine pointwise/vectorized indexer (non-axis dim != 1) is rejected with a clear error. This makes the real ITS_LIVE x+y mosaic work via native reindex/concat. Co-Authored-By: Claude Opus 4.8 (1M context) --- virtualizarr/manifests/indexing.py | 27 ++++++++++++++++++++-- virtualizarr/tests/test_native_align.py | 30 +++++++++++++++++++++++++ 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/virtualizarr/manifests/indexing.py b/virtualizarr/manifests/indexing.py index ebe2e914..5ad1d727 100644 --- a/virtualizarr/manifests/indexing.py +++ b/virtualizarr/manifests/indexing.py @@ -145,11 +145,33 @@ def apply_indexer( return output_arr -def _is_reindex_indexer(indexer_1d: T_BasicIndexer_1d) -> bool: +def _is_reindex_indexer(indexer_1d: T_BasicIndexer_1d) -> TypeGuard[np.ndarray]: """True if this 1D indexer is an integer array (xarray's reindex indexer).""" return isinstance(indexer_1d, np.ndarray) and indexer_1d.dtype.kind in ("i", "u") +def _collapse_outer_indexer( + indexer_1d: T_BasicIndexer_1d, axis: int +) -> T_BasicIndexer_1d: + """Collapse a broadcast reindex indexer back to its 1D per-axis form. + + When xarray reindexes more than one dimension at once it sends a single + *broadcast* (vectorized) indexer: one integer array per axis, shaped to the + array's rank with the real labels on its own axis and length 1 everywhere + else (e.g. ``(N, 1)`` for axis 0 and ``(1, M)`` for axis 1). Because reindex + indexing is orthogonal, each such array reshapes losslessly to the 1D + indexer for ``axis``. int/slice indexers and already-1D arrays pass through. + """ + if not (isinstance(indexer_1d, np.ndarray) and indexer_1d.ndim > 1): + return indexer_1d + if any(size != 1 for ax, size in enumerate(indexer_1d.shape) if ax != axis): + raise NotImplementedError( + "Pointwise (vectorized) fancy indexing is not supported on a " + f"ManifestArray; received a {indexer_1d.shape} indexer for axis {axis}." + ) + return indexer_1d.reshape(-1) + + def _is_full_slice(indexer_1d: T_BasicIndexer_1d, length: int) -> bool: """True if this slice selects the whole axis (a no-op).""" if not isinstance(indexer_1d, slice): @@ -175,7 +197,8 @@ def _apply_reindex( from virtualizarr.manifests.reindex import chunk_map_from_indexer result = marr - for axis, indexer_1d in enumerate(indexers): + for axis, raw_indexer_1d in enumerate(indexers): + indexer_1d = _collapse_outer_indexer(raw_indexer_1d, axis) if _is_reindex_indexer(indexer_1d): chunk_size = result.metadata.chunks[axis] try: diff --git a/virtualizarr/tests/test_native_align.py b/virtualizarr/tests/test_native_align.py index 79ed7d5f..11055812 100644 --- a/virtualizarr/tests/test_native_align.py +++ b/virtualizarr/tests/test_native_align.py @@ -142,6 +142,36 @@ def test_reindex_non_chunk_aligned_raises(self): assert isinstance(excinfo.value.__cause__, NotImplementedError) +class TestNativeReindexMultiAxis: + def test_reindex_two_axes_simultaneously(self): + # A 2D outer-join reindexes BOTH axes in one call. xarray then sends a + # single broadcast (vectorized) indexer -- one N-D array per axis, shaped + # (Nt, 1) and (1, Nx) -- rather than two 1D indexers. The chunk-grid remap + # must collapse each broadcast component back to its 1D per-axis indexer, + # or it misreads the array's rank (the real ITS_LIVE x+y mosaic case). + import obstore as obs + + store = obs.store.MemoryStore() + # 2x2 grid, row-major (time, x): (t0,x0)=1 (t0,x1)=2 (t1,x0)=3 (t1,x1)=1 + obs.put(store, "a.bin", F32[1.0] + F32[2.0] + F32[3.0] + F32[1.0]) + ds = _ds_2d([0, 1], [0, 1], "memory:///a.bin") # x=[0,1], time=[0,1] + + result = ds.reindex(time=[0, 1, 2], x=[0, 1, 2]) + + assert isinstance(result["foo"].data, ManifestArray) + assert list(result.time.values) == [0, 1, 2] + assert list(result.x.values) == [0, 1, 2] + + ms = ManifestStore( + group=ManifestGroup(arrays={"foo": result["foo"].data}), + registry=ObjectStoreRegistry({"memory://": store}), + ) + vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values + np.testing.assert_array_equal(vals[:2, :2], [[1.0, 2.0], [3.0, 1.0]]) + assert np.isnan(vals[2, :]).all() # appended time row reads back as fill + assert np.isnan(vals[:, 2]).all() # appended x column reads back as fill + + class TestNativeAlign: def test_align_outer_union(self): a = _ds_1d([0, 1, 2], "/a.nc") From a62f6e6eb8bcd524a697aad7e0f9216a40117d98 Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Tue, 30 Jun 2026 16:35:15 -0400 Subject: [PATCH 3/7] Move function level imports up. --- virtualizarr/manifests/indexing.py | 3 +-- virtualizarr/tests/test_native_align.py | 9 +-------- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/virtualizarr/manifests/indexing.py b/virtualizarr/manifests/indexing.py index 5ad1d727..435189cc 100644 --- a/virtualizarr/manifests/indexing.py +++ b/virtualizarr/manifests/indexing.py @@ -7,6 +7,7 @@ from virtualizarr.manifests.array_api import expand_dims from virtualizarr.manifests.manifest import ChunkManifest +from virtualizarr.manifests.reindex import chunk_map_from_indexer from virtualizarr.manifests.utils import copy_and_replace_metadata # indexer with only basic selectors, no new axes or ellipsis @@ -194,8 +195,6 @@ def _apply_reindex( full-axis slice (the shape xarray's reindex/alignment produces); combining reindexing with chunk subsetting in a single call is not supported. """ - from virtualizarr.manifests.reindex import chunk_map_from_indexer - result = marr for axis, raw_indexer_1d in enumerate(indexers): indexer_1d = _collapse_outer_indexer(raw_indexer_1d, axis) diff --git a/virtualizarr/tests/test_native_align.py b/virtualizarr/tests/test_native_align.py index 11055812..bc17be3d 100644 --- a/virtualizarr/tests/test_native_align.py +++ b/virtualizarr/tests/test_native_align.py @@ -7,6 +7,7 @@ """ import numpy as np +import obstore as obs import pytest import xarray as xr from obspec_utils.registry import ObjectStoreRegistry @@ -87,8 +88,6 @@ def test_reindex_stays_lazy(self): } def test_reindex_float_reads_back_nan(self): - import obstore as obs - store = obs.store.MemoryStore() obs.put(store, "a.bin", F32[1.0] + F32[2.0] + F32[3.0]) ds = _ds_1d([0, 1, 2], "memory:///a.bin") @@ -112,8 +111,6 @@ def test_reindex_float_reads_back_nan(self): def test_reindex_int_keeps_dtype_and_sentinel(self): # xarray would normally promote int->float+NaN; the manifest keeps int32 # and its declared sentinel fill instead. - import obstore as obs - store = obs.store.MemoryStore() obs.put(store, "a.bin", I32[10] + I32[20] + I32[30]) ds = _ds_1d([0, 1, 2], "memory:///a.bin", dtype="int32", fill_value=-9999) @@ -149,8 +146,6 @@ def test_reindex_two_axes_simultaneously(self): # (Nt, 1) and (1, Nx) -- rather than two 1D indexers. The chunk-grid remap # must collapse each broadcast component back to its 1D per-axis indexer, # or it misreads the array's rank (the real ITS_LIVE x+y mosaic case). - import obstore as obs - store = obs.store.MemoryStore() # 2x2 grid, row-major (time, x): (t0,x0)=1 (t0,x1)=2 (t1,x0)=3 (t1,x1)=1 obs.put(store, "a.bin", F32[1.0] + F32[2.0] + F32[3.0] + F32[1.0]) @@ -185,8 +180,6 @@ def test_align_outer_union(self): def test_align_exclude_then_concat_spatial(self): # the ITS_LIVE pattern: align on the spatial dim (exclude time), concat on time - import obstore as obs - store = obs.store.MemoryStore() obs.put(store, "a.bin", F32[1.0] + F32[2.0]) # time=0, x=0,1 obs.put(store, "b.bin", F32[3.0] + F32[1.0]) # time=1, x=2,3 From da8a8486bfc1b315cef326ec43d78be5c58a3939 Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Thu, 2 Jul 2026 11:41:19 -0400 Subject: [PATCH 4/7] Refactor reindexing tests to only touch public interfaces. --- .../tests/test_manifests/test_reindex.py | 177 +++++++++++++----- 1 file changed, 125 insertions(+), 52 deletions(-) diff --git a/virtualizarr/tests/test_manifests/test_reindex.py b/virtualizarr/tests/test_manifests/test_reindex.py index 7e0aa90f..a79f82ae 100644 --- a/virtualizarr/tests/test_manifests/test_reindex.py +++ b/virtualizarr/tests/test_manifests/test_reindex.py @@ -1,8 +1,22 @@ +""" +Reindex/alignment behaviour exercised through the *public* ManifestArray API. + +xarray's reindex/alignment machinery conforms an array to new labels by indexing +it with an integer "gather" indexer that has ``-1`` at every missing position +(e.g. ``[0, 1, 2, -1, -1]``). ManifestArray supports that indexer at the array +level -- ``marr[indexer]`` -- by remapping the chunk grid: a run of ``-1`` that +covers a whole chunk becomes a null-path chunk (reads back as ``fill_value``), +and a chunk-aligned run of real indices keeps its chunk reference. Nothing that +would split a chunk is allowed. + +These tests drive that public indexing path directly rather than the private +``chunk_map_from_indexer``/``_reindex_axis`` helpers that implement it. +""" + import numpy as np import pytest from virtualizarr.manifests import ChunkManifest, ManifestArray -from virtualizarr.manifests.reindex import chunk_map_from_indexer from virtualizarr.manifests.utils import create_v3_array_metadata @@ -24,59 +38,65 @@ def _marr(shape, chunks, entries, dims): ) -class TestChunkMapFromIndexer: - def test_append_chunk_size_1(self): - assert chunk_map_from_indexer(_idx(0, 1, 2, -1, -1), 1, 3) == [ - 0, - 1, - 2, - None, - None, - ] - - def test_chunked_append(self): - # source chunks [0,1],[2,3]; one new all-missing chunk appended - assert chunk_map_from_indexer(_idx(0, 1, 2, 3, -1, -1), 2, 4) == [0, 1, None] - - def test_prepend(self): - assert chunk_map_from_indexer(_idx(-1, -1, 0, 1), 1, 2) == [None, None, 0, 1] - - def test_gap_fill_insert_whole_chunk(self): - # [0,1] [missing] [2,3] - assert chunk_map_from_indexer(_idx(0, 1, -1, -1, 2, 3), 2, 4) == [0, None, 1] - - def test_whole_chunk_reverse_c1(self): - assert chunk_map_from_indexer(_idx(2, 1, 0), 1, 3) == [2, 1, 0] +class TestReindexGatherKeepsChunks: + """An integer gather indexer remaps the chunk grid, lazily, keeping refs.""" - def test_whole_chunk_reverse_c2(self): - assert chunk_map_from_indexer(_idx(2, 3, 0, 1), 2, 4) == [1, 0] - - def test_trailing_partial_source_chunk(self): - # source len 5, C=2 -> chunk sizes 2,2,1; selecting just the size-1 tail - assert chunk_map_from_indexer(_idx(4), 2, 5) == [2] - - def test_raise_mixed_present_and_missing(self): - with pytest.raises(NotImplementedError, match="split"): - chunk_map_from_indexer(_idx(0, 1, 2, -1), 2, 4) + def test_append_null_chunks(self): + # chunk size 1: every missing (-1) position becomes its own null chunk + marr = _marr( + (3,), + (1,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 4}, + "1": {"path": "/a.nc", "offset": 4, "length": 4}, + "2": {"path": "/a.nc", "offset": 8, "length": 4}, + }, + ["x"], + ) + result = marr[_idx(0, 1, 2, -1, -1)] - def test_raise_sub_chunk_reorder(self): - with pytest.raises(NotImplementedError, match="split"): - chunk_map_from_indexer(_idx(1, 0, 2, 3), 2, 4) + assert result.shape == (5,) + assert result.metadata.chunks == (1,) + # appended positions are absent from the manifest -> null chunks + assert sorted(result.manifest.dict()) == ["0", "1", "2"] - def test_raise_unaligned_start(self): - with pytest.raises(NotImplementedError, match="split"): - chunk_map_from_indexer(_idx(1, 2), 2, 4) + def test_chunked_append(self): + # source chunks [0,1],[2,3]; one all-missing chunk appended -> one null + marr = _marr( + (4,), + (2,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 8}, + "1": {"path": "/b.nc", "offset": 8, "length": 8}, + }, + ["x"], + ) + result = marr[_idx(0, 1, 2, 3, -1, -1)] - def test_raise_partial_target_of_full_source(self): - # taking only half of source chunk [2,3] would split it - with pytest.raises(NotImplementedError, match="split"): - chunk_map_from_indexer(_idx(0, 1, 2), 2, 4) + assert result.shape == (6,) + assert result.manifest.dict() == { + "0": {"path": "file:///a.nc", "offset": 0, "length": 8}, + "1": {"path": "file:///b.nc", "offset": 8, "length": 8}, + } + def test_prepend_null_chunks(self): + marr = _marr( + (2,), + (1,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 4}, + "1": {"path": "/a.nc", "offset": 4, "length": 4}, + }, + ["x"], + ) + result = marr[_idx(-1, -1, 0, 1)] -class TestReindexAxis: - """ManifestArray._reindex_axis: the chunk-grid remap the indexer hook drives.""" + assert result.shape == (4,) + # real chunks shifted to slots 2,3; slots 0,1 are null + assert sorted(result.manifest.dict()) == ["2", "3"] - def test_pad_with_null_chunk(self): + def test_insert_whole_chunk_gap(self): + # [0,1] [missing] [2,3]: a whole null chunk inserted between two real ones marr = _marr( (4,), (2,), @@ -86,16 +106,16 @@ def test_pad_with_null_chunk(self): }, ["x"], ) - result = marr._reindex_axis(axis=0, chunk_map=[0, None, 1], new_size=6) + result = marr[_idx(0, 1, -1, -1, 2, 3)] assert result.shape == (6,) - assert result.metadata.chunks == (2,) assert result.manifest.dict() == { "0": {"path": "file:///a.nc", "offset": 0, "length": 8}, "2": {"path": "file:///b.nc", "offset": 8, "length": 8}, } def test_whole_chunk_reorder(self): + # reversing two whole chunks just swaps their references; no bytes read marr = _marr( (4,), (2,), @@ -105,14 +125,34 @@ def test_whole_chunk_reorder(self): }, ["x"], ) - result = marr._reindex_axis(axis=0, chunk_map=[1, 0], new_size=4) + result = marr[_idx(2, 3, 0, 1)] assert result.manifest.dict() == { "0": {"path": "file:///b.nc", "offset": 8, "length": 8}, "1": {"path": "file:///a.nc", "offset": 0, "length": 8}, } - def test_only_target_axis_remapped_2d(self): + def test_select_trailing_partial_chunk(self): + # source len 5, chunk 2 -> chunk sizes 2,2,1; selecting just the size-1 tail + marr = _marr( + (5,), + (2,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 8}, + "1": {"path": "/a.nc", "offset": 8, "length": 8}, + "2": {"path": "/a.nc", "offset": 16, "length": 4}, + }, + ["x"], + ) + result = marr[_idx(4)] + + assert result.shape == (1,) + assert result.manifest.dict() == { + "0": {"path": "file:///a.nc", "offset": 16, "length": 4}, + } + + def test_only_indexed_axis_remapped_2d(self): + # a slice on axis 0 + gather on axis 1: only axis 1's grid changes marr = _marr( (2, 4), (2, 2), @@ -122,7 +162,7 @@ def test_only_target_axis_remapped_2d(self): }, ["y", "x"], ) - result = marr._reindex_axis(axis=1, chunk_map=[0, 1, None], new_size=6) + result = marr[:, _idx(0, 1, 2, 3, -1, -1)] assert result.shape == (2, 6) assert result.metadata.chunks == (2, 2) @@ -130,3 +170,36 @@ def test_only_target_axis_remapped_2d(self): "0.0": {"path": "file:///a.nc", "offset": 0, "length": 16}, "0.1": {"path": "file:///b.nc", "offset": 16, "length": 16}, } + + +class TestReindexGatherRejectsChunkSplits: + """Anything that would split a chunk raises, naming the offending axis.""" + + @pytest.fixture + def marr(self): + return _marr( + (4,), + (2,), + { + "0": {"path": "/a.nc", "offset": 0, "length": 8}, + "1": {"path": "/b.nc", "offset": 8, "length": 8}, + }, + ["x"], + ) + + def test_mixed_present_and_missing_in_a_chunk(self, marr): + with pytest.raises(NotImplementedError, match="chunk boundaries"): + marr[_idx(0, 1, 2, -1)] + + def test_sub_chunk_reorder(self, marr): + with pytest.raises(NotImplementedError, match="chunk boundaries"): + marr[_idx(1, 0, 2, 3)] + + def test_unaligned_start(self, marr): + with pytest.raises(NotImplementedError, match="chunk boundaries"): + marr[_idx(1, 2)] + + def test_partial_target_of_full_source(self, marr): + # taking only part of source chunk [2,3] would split it + with pytest.raises(NotImplementedError, match="chunk boundaries"): + marr[_idx(0, 1, 2)] From 059248399fcca46788ac72ff82e36e1c200f12d7 Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Thu, 2 Jul 2026 15:42:15 -0400 Subject: [PATCH 5/7] Refactor xarray align tests to use existing fixtures. --- virtualizarr/tests/test_native_align.py | 365 +++++++++++++----------- 1 file changed, 205 insertions(+), 160 deletions(-) diff --git a/virtualizarr/tests/test_native_align.py b/virtualizarr/tests/test_native_align.py index bc17be3d..6c543558 100644 --- a/virtualizarr/tests/test_native_align.py +++ b/virtualizarr/tests/test_native_align.py @@ -1,136 +1,128 @@ """ End-to-end tests that native xarray alignment machinery works over ManifestArrays. -The point of the indexer/where/result_type hooks is that users call plain -`xr.reindex`, `xr.align`, and `xr.concat` — no VirtualiZarr-specific API — and the -fill stays virtual (null-path chunks that read back as the array's fill_value). +Users call plain ``xr.reindex`` / ``xr.align`` / ``xr.concat`` -- no +VirtualiZarr-specific API -- and the fill stays virtual: null-path chunks that +read back as the array's ``fill_value``. + +Structure/laziness cases build ManifestArrays directly via the ``array_v3_metadata`` +fixture (as ``TestConcat`` in test_xarray.py does). Value-readback cases open real +tiled NetCDF files via the ``netcdf4_files_factory_2d`` fixture + a parser (as +``TestCombine`` does) and read the reindexed result back through a ManifestStore. """ import numpy as np -import obstore as obs import pytest import xarray as xr -from obspec_utils.registry import ObjectStoreRegistry +from virtualizarr import open_virtual_dataset from virtualizarr.manifests import ( ChunkManifest, ManifestArray, ManifestGroup, ManifestStore, ) -from virtualizarr.manifests.utils import create_v3_array_metadata - -CODECS = [{"configuration": {"endian": "little"}, "name": "bytes"}] -F32 = { - 1.0: b"\x00\x00\x80\x3f", - 2.0: b"\x00\x00\x00\x40", - 3.0: b"\x00\x00\x40\x40", -} -I32 = {10: b"\x0a\x00\x00\x00", 20: b"\x14\x00\x00\x00", 30: b"\x1e\x00\x00\x00"} - - -def _ds_1d(x, path, dtype="float32", fill_value=float("nan"), chunk=1, name="foo"): - n = len(x) - n_chunks = -(-n // chunk) - entries = { - str(i): {"path": path, "offset": i * chunk * 4, "length": chunk * 4} - for i in range(n_chunks) - } - metadata = create_v3_array_metadata( - shape=(n,), - chunk_shape=(chunk,), - data_type=np.dtype(dtype), - codecs=CODECS, - fill_value=fill_value, - dimension_names=["x"], - ) - marr = ManifestArray( - metadata=metadata, chunkmanifest=ChunkManifest(entries=entries) - ) - return xr.Dataset({name: xr.Variable(["x"], marr)}, coords={"x": list(x)}) - - -def _ds_2d(x, t, path, name="foo"): - nt, nx = len(t), len(x) - entries = {} - k = 0 - for it in range(nt): - for ix in range(nx): - entries[f"{it}.{ix}"] = {"path": path, "offset": k * 4, "length": 4} - k += 1 - metadata = create_v3_array_metadata( - shape=(nt, nx), - chunk_shape=(1, 1), - data_type=np.dtype("float32"), - codecs=CODECS, - fill_value=float("nan"), - dimension_names=["time", "x"], - ) - marr = ManifestArray( - metadata=metadata, chunkmanifest=ChunkManifest(entries=entries) - ) - return xr.Dataset( - {name: xr.Variable(["time", "x"], marr)}, coords={"x": list(x), "time": list(t)} +from virtualizarr.manifests.utils import copy_and_replace_metadata +from virtualizarr.parsers import HDFParser +from virtualizarr.tests import requires_hdf5plugin, requires_imagecodecs + +LOADABLE_COORDS = ["time", "lat", "lon"] + + +def _read_back(marr: ManifestArray, dims: list[str], registry) -> np.ndarray: + """Materialise a virtual ManifestArray's values by wrapping it in a + ManifestStore and opening it. + + Reindex preserves whatever ``dimension_names`` the source metadata carried; + parser-built arrays carry ``None``, which ``xr.open_zarr`` requires, so set + them here for the read. + """ + metadata = copy_and_replace_metadata(marr.metadata, new_dimension_names=dims) + ms = ManifestStore( + group=ManifestGroup( + arrays={"v": ManifestArray(metadata=metadata, chunkmanifest=marr.manifest)} + ), + registry=registry, ) + return xr.open_zarr(ms, consolidated=False, zarr_format=3)["v"].values class TestNativeReindex: - def test_reindex_stays_lazy(self): - ds = _ds_1d([0, 1, 2], "/a.nc") + def test_reindex_stays_lazy(self, array_v3_metadata): + metadata = array_v3_metadata(shape=(3,), chunks=(1,), dimension_names=["x"]) + marr = ManifestArray( + metadata=metadata, + chunkmanifest=ChunkManifest( + entries={ + "0": {"path": "/a.nc", "offset": 0, "length": 4}, + "1": {"path": "/a.nc", "offset": 4, "length": 4}, + "2": {"path": "/a.nc", "offset": 8, "length": 4}, + } + ), + ) + ds = xr.Dataset({"foo": ("x", marr)}, coords={"x": [0, 1, 2]}) + result = ds.reindex(x=[0, 1, 2, 3, 4]) # plain xarray, no vz API assert isinstance(result["foo"].data, ManifestArray) assert list(result.x.values) == [0, 1, 2, 3, 4] + # appended positions become null chunks (absent from the manifest) assert result["foo"].data.manifest.dict() == { "0": {"path": "file:///a.nc", "offset": 0, "length": 4}, "1": {"path": "file:///a.nc", "offset": 4, "length": 4}, "2": {"path": "file:///a.nc", "offset": 8, "length": 4}, } - def test_reindex_float_reads_back_nan(self): - store = obs.store.MemoryStore() - obs.put(store, "a.bin", F32[1.0] + F32[2.0] + F32[3.0]) - ds = _ds_1d([0, 1, 2], "memory:///a.bin") - - reindexed = ds.reindex(x=[0, 1, 2, 3, 4]) - - ms = ManifestStore( - group=ManifestGroup(arrays={"foo": reindexed["foo"].data}), - registry=ObjectStoreRegistry({"memory://": store}), + def test_reindex_two_axes_simultaneously(self, array_v3_metadata): + # A 2D outer-join reindexes BOTH axes in one call, which xarray sends as a + # single broadcast (vectorized) indexer -- one N-D array per axis, shaped + # (Nt, 1) and (1, Nx) -- rather than two 1D indexers. The chunk-grid remap + # must collapse each back to its 1D per-axis indexer, else it misreads the + # array's rank (the real ITS_LIVE x+y mosaic case). + metadata = array_v3_metadata( + shape=(2, 2), chunks=(1, 1), dimension_names=["time", "x"] ) - vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values - np.testing.assert_array_equal(vals[:3], [1.0, 2.0, 3.0]) - assert np.isnan(vals[3:]).all() - - # xarray internally builds a NaN fill scalar and casts it toward the int - # dtype (a value we never use, since our np.where returns the array as-is), - # which emits a harmless "invalid value encountered in cast" RuntimeWarning. - @pytest.mark.filterwarnings( - "ignore:invalid value encountered in cast:RuntimeWarning" - ) - def test_reindex_int_keeps_dtype_and_sentinel(self): - # xarray would normally promote int->float+NaN; the manifest keeps int32 - # and its declared sentinel fill instead. - store = obs.store.MemoryStore() - obs.put(store, "a.bin", I32[10] + I32[20] + I32[30]) - ds = _ds_1d([0, 1, 2], "memory:///a.bin", dtype="int32", fill_value=-9999) - - reindexed = ds.reindex(x=[0, 1, 2, 3, 4]) - assert reindexed["foo"].dtype == np.dtype("int32") - assert isinstance(reindexed["foo"].data, ManifestArray) - - ms = ManifestStore( - group=ManifestGroup(arrays={"foo": reindexed["foo"].data}), - registry=ObjectStoreRegistry({"memory://": store}), + marr = ManifestArray( + metadata=metadata, + chunkmanifest=ChunkManifest( + entries={ + "0.0": {"path": "/a.nc", "offset": 0, "length": 4}, + "0.1": {"path": "/a.nc", "offset": 4, "length": 4}, + "1.0": {"path": "/a.nc", "offset": 8, "length": 4}, + "1.1": {"path": "/a.nc", "offset": 12, "length": 4}, + } + ), + ) + ds = xr.Dataset( + {"foo": (["time", "x"], marr)}, coords={"time": [0, 1], "x": [0, 1]} ) - vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values - np.testing.assert_array_equal(vals, [10, 20, 30, -9999, -9999]) - assert vals.dtype == np.dtype("int32") - def test_reindex_non_chunk_aligned_raises(self): - ds = _ds_1d([0, 1, 2], "/a.nc", chunk=2) + result = ds.reindex(time=[0, 1, 2], x=[0, 1, 2]) + + assert isinstance(result["foo"].data, ManifestArray) + assert result["foo"].shape == (3, 3) + # the four real chunks are kept; the appended row and column are null + assert result["foo"].data.manifest.dict() == { + "0.0": {"path": "file:///a.nc", "offset": 0, "length": 4}, + "0.1": {"path": "file:///a.nc", "offset": 4, "length": 4}, + "1.0": {"path": "file:///a.nc", "offset": 8, "length": 4}, + "1.1": {"path": "file:///a.nc", "offset": 12, "length": 4}, + } + + def test_reindex_non_chunk_aligned_raises(self, array_v3_metadata): + metadata = array_v3_metadata(shape=(3,), chunks=(2,), dimension_names=["x"]) + marr = ManifestArray( + metadata=metadata, + chunkmanifest=ChunkManifest( + entries={ + "0": {"path": "/a.nc", "offset": 0, "length": 8}, + "1": {"path": "/a.nc", "offset": 8, "length": 4}, + } + ), + ) + ds = xr.Dataset({"foo": ("x", marr)}, coords={"x": [0, 1, 2]}) # the deep error is wrapped with call-site context: the operation, the - # axis, and what to do — and the low-level reason is chained. + # axis, and what to do -- and the low-level reason is chained. with pytest.raises( NotImplementedError, match="align/reindex this virtual array along axis 0" ) as excinfo: @@ -139,38 +131,21 @@ def test_reindex_non_chunk_aligned_raises(self): assert isinstance(excinfo.value.__cause__, NotImplementedError) -class TestNativeReindexMultiAxis: - def test_reindex_two_axes_simultaneously(self): - # A 2D outer-join reindexes BOTH axes in one call. xarray then sends a - # single broadcast (vectorized) indexer -- one N-D array per axis, shaped - # (Nt, 1) and (1, Nx) -- rather than two 1D indexers. The chunk-grid remap - # must collapse each broadcast component back to its 1D per-axis indexer, - # or it misreads the array's rank (the real ITS_LIVE x+y mosaic case). - store = obs.store.MemoryStore() - # 2x2 grid, row-major (time, x): (t0,x0)=1 (t0,x1)=2 (t1,x0)=3 (t1,x1)=1 - obs.put(store, "a.bin", F32[1.0] + F32[2.0] + F32[3.0] + F32[1.0]) - ds = _ds_2d([0, 1], [0, 1], "memory:///a.bin") # x=[0,1], time=[0,1] - - result = ds.reindex(time=[0, 1, 2], x=[0, 1, 2]) - - assert isinstance(result["foo"].data, ManifestArray) - assert list(result.time.values) == [0, 1, 2] - assert list(result.x.values) == [0, 1, 2] - - ms = ManifestStore( - group=ManifestGroup(arrays={"foo": result["foo"].data}), - registry=ObjectStoreRegistry({"memory://": store}), - ) - vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values - np.testing.assert_array_equal(vals[:2, :2], [[1.0, 2.0], [3.0, 1.0]]) - assert np.isnan(vals[2, :]).all() # appended time row reads back as fill - assert np.isnan(vals[:, 2]).all() # appended x column reads back as fill - - class TestNativeAlign: - def test_align_outer_union(self): - a = _ds_1d([0, 1, 2], "/a.nc") - b = _ds_1d([3, 4, 5], "/b.nc") + def test_align_outer_union(self, array_v3_metadata): + def _ds(path, xs): + metadata = array_v3_metadata(shape=(3,), chunks=(1,), dimension_names=["x"]) + manifest = ChunkManifest( + entries={ + str(i): {"path": path, "offset": i * 4, "length": 4} + for i in range(3) + } + ) + marr = ManifestArray(metadata=metadata, chunkmanifest=manifest) + return xr.Dataset({"foo": ("x", marr)}, coords={"x": xs}) + + a = _ds("/a.nc", [0, 1, 2]) + b = _ds("/b.nc", [3, 4, 5]) ra, rb = xr.align(a, b, join="outer") @@ -178,39 +153,109 @@ def test_align_outer_union(self): assert isinstance(rb["foo"].data, ManifestArray) assert list(ra.x.values) == [0, 1, 2, 3, 4, 5] - def test_align_exclude_then_concat_spatial(self): - # the ITS_LIVE pattern: align on the spatial dim (exclude time), concat on time - store = obs.store.MemoryStore() - obs.put(store, "a.bin", F32[1.0] + F32[2.0]) # time=0, x=0,1 - obs.put(store, "b.bin", F32[3.0] + F32[1.0]) # time=1, x=2,3 - a = _ds_2d([0, 1], [0], "memory:///a.bin") - b = _ds_2d([2, 3], [1], "memory:///b.bin") - - ra, rb = xr.align(a, b, join="outer", exclude=["time"]) - # time is left alone; x is unioned - assert list(ra.time.values) == [0] - assert list(rb.time.values) == [1] - assert list(ra.x.values) == [0, 1, 2, 3] - - cube = xr.concat([ra, rb], dim="time", join="exact") - assert cube.sizes == {"time": 2, "x": 4} - assert isinstance(cube["foo"].data, ManifestArray) - - ms = ManifestStore( - group=ManifestGroup(arrays={"foo": cube["foo"].data}), - registry=ObjectStoreRegistry({"memory://": store}), - ) - vals = xr.open_zarr(ms, consolidated=False, zarr_format=3)["foo"].values - np.testing.assert_array_equal(vals[0, :2], [1.0, 2.0]) - assert np.isnan(vals[0, 2:]).all() - assert np.isnan(vals[1, :2]).all() - np.testing.assert_array_equal(vals[1, 2:], [3.0, 1.0]) + +@requires_hdf5plugin +@requires_imagecodecs +class TestNativeReindexReadback: + """Reindex fill reads back as the array's fill_value, verified over real files. + + Uses the air_temperature tiles from ``netcdf4_files_factory_2d`` (as + ``TestCombine`` does): two time slices x two lat bands. + + NOTE: air's ``lat`` is descending, so ``xr.align(join="outer")`` -- which sorts + the union ascending -- would require reversing within a chunk (not lazy). The + chunk-aligned, lazy path is an explicit ``reindex`` onto a descending union, + which is also the real ITS_LIVE mosaic pattern. + """ + + def test_spatial_reindex_then_concat_on_time_keeps_int_sentinel( + self, netcdf4_files_factory_2d, local_registry + ): + # air_temperature is stored int16 with sentinel fill -32767. Two tiles that + # differ in BOTH time slice and lat band (the ITS_LIVE granule shape): + # reindex each onto the union lat, then concat on time. + early_north, _, _, late_south = netcdf4_files_factory_2d() + parser = HDFParser() + with ( + open_virtual_dataset( + url=early_north, + registry=local_registry, + parser=parser, + loadable_variables=LOADABLE_COORDS, + ) as en, + open_virtual_dataset( + url=late_south, + registry=local_registry, + parser=parser, + loadable_variables=LOADABLE_COORDS, + ) as ls, + ): + union_lat = np.concatenate([en.lat.values, ls.lat.values]) # descending + cube = xr.concat( + [en.reindex(lat=union_lat), ls.reindex(lat=union_lat)], + dim="time", + join="exact", + ) + assert isinstance(cube["air"].data, ManifestArray) + assert cube["air"].dtype == np.dtype("int16") # not promoted to float + + nt, nlat = en.sizes["time"], en.sizes["lat"] + fill = cube["air"].data.metadata.fill_value + vals = _read_back(cube["air"].data, ["time", "lat", "lon"], local_registry) + en_vals = _read_back(en["air"].data, ["time", "lat", "lon"], local_registry) + ls_vals = _read_back(ls["air"].data, ["time", "lat", "lon"], local_registry) + + np.testing.assert_array_equal(vals[:nt, :nlat, :], en_vals) # early/north + np.testing.assert_array_equal(vals[nt:, nlat:, :], ls_vals) # late/south + assert (vals[:nt, nlat:, :] == fill).all() # early/south -> sentinel + assert (vals[nt:, :nlat, :] == fill).all() # late/north -> sentinel + + def test_float_fill_reads_back_nan(self, netcdf4_files_factory_2d, local_registry): + # write the tiles as float32/NaN-fill so the padded band reads back as NaN + encoding = {"air": {"dtype": "float32", "_FillValue": np.float32("nan")}} + north, _, south, _ = netcdf4_files_factory_2d(encoding=encoding) + parser = HDFParser() + with ( + open_virtual_dataset( + url=north, + registry=local_registry, + parser=parser, + loadable_variables=LOADABLE_COORDS, + ) as n, + open_virtual_dataset( + url=south, + registry=local_registry, + parser=parser, + loadable_variables=LOADABLE_COORDS, + ) as s, + ): + assert n["air"].dtype == np.dtype("float32") + union_lat = np.concatenate([n.lat.values, s.lat.values]) + reindexed = n.reindex(lat=union_lat) + assert isinstance(reindexed["air"].data, ManifestArray) + + nlat = n.sizes["lat"] + vals = _read_back( + reindexed["air"].data, ["time", "lat", "lon"], local_registry + ) + assert np.isfinite(vals[:, :nlat, :]).all() # real band + assert np.isnan(vals[:, nlat:, :]).all() # padded band -> NaN class TestGeneralWhereStillRejected: - def test_boolean_masking_raises_clearly(self): + def test_boolean_masking_raises_clearly(self, array_v3_metadata): # we must NOT silently no-op general where(); it should raise, since it # would require materializing values. - ds = _ds_1d([0, 1, 2], "/a.nc") + metadata = array_v3_metadata(shape=(3,), chunks=(1,), dimension_names=["x"]) + marr = ManifestArray( + metadata=metadata, + chunkmanifest=ChunkManifest( + entries={ + str(i): {"path": "/a.nc", "offset": i * 4, "length": 4} + for i in range(3) + } + ), + ) + ds = xr.Dataset({"foo": ("x", marr)}, coords={"x": [0, 1, 2]}) with pytest.raises(NotImplementedError, match="materializing|where"): ds["foo"].where(ds["x"] > 0).compute() From 28e241071e1055ae4188db29a5886c9a472eee20 Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Thu, 2 Jul 2026 15:45:55 -0400 Subject: [PATCH 6/7] Remove native naming style that is no longer relevant. --- .../tests/{test_native_align.py => test_align.py} | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) rename virtualizarr/tests/{test_native_align.py => test_align.py} (98%) diff --git a/virtualizarr/tests/test_native_align.py b/virtualizarr/tests/test_align.py similarity index 98% rename from virtualizarr/tests/test_native_align.py rename to virtualizarr/tests/test_align.py index 6c543558..9656f7d3 100644 --- a/virtualizarr/tests/test_native_align.py +++ b/virtualizarr/tests/test_align.py @@ -1,5 +1,5 @@ """ -End-to-end tests that native xarray alignment machinery works over ManifestArrays. +End-to-end tests that xarray alignment machinery works over ManifestArrays. Users call plain ``xr.reindex`` / ``xr.align`` / ``xr.concat`` -- no VirtualiZarr-specific API -- and the fill stays virtual: null-path chunks that @@ -47,7 +47,7 @@ def _read_back(marr: ManifestArray, dims: list[str], registry) -> np.ndarray: return xr.open_zarr(ms, consolidated=False, zarr_format=3)["v"].values -class TestNativeReindex: +class TestReindex: def test_reindex_stays_lazy(self, array_v3_metadata): metadata = array_v3_metadata(shape=(3,), chunks=(1,), dimension_names=["x"]) marr = ManifestArray( @@ -131,7 +131,7 @@ def test_reindex_non_chunk_aligned_raises(self, array_v3_metadata): assert isinstance(excinfo.value.__cause__, NotImplementedError) -class TestNativeAlign: +class TestAlign: def test_align_outer_union(self, array_v3_metadata): def _ds(path, xs): metadata = array_v3_metadata(shape=(3,), chunks=(1,), dimension_names=["x"]) @@ -156,7 +156,7 @@ def _ds(path, xs): @requires_hdf5plugin @requires_imagecodecs -class TestNativeReindexReadback: +class TestReindexReadback: """Reindex fill reads back as the array's fill_value, verified over real files. Uses the air_temperature tiles from ``netcdf4_files_factory_2d`` (as From c82094ff808bbe588a9460eec1485a1937087b4b Mon Sep 17 00:00:00 2001 From: sharkinsspatial Date: Thu, 9 Jul 2026 11:56:42 -0500 Subject: [PATCH 7/7] Add more detailed usage documentation and a complete its_live example. --- docs/how_to/examples.md | 4 + docs/how_to/grid_alignment.png | Bin 0 -> 38148 bytes docs/how_to/validation.md | 11 +- examples/V2/README.md | 4 + examples/V2/its_live.ipynb | 288 +++++++++++++++++++++++++++++++++ 5 files changed, 306 insertions(+), 1 deletion(-) create mode 100644 docs/how_to/grid_alignment.png create mode 100644 examples/V2/its_live.ipynb diff --git a/docs/how_to/examples.md b/docs/how_to/examples.md index 9c319716..700e9822 100644 --- a/docs/how_to/examples.md +++ b/docs/how_to/examples.md @@ -15,6 +15,10 @@ This notebook accompanies the [_Old format, no problem!: Cloud-optimizing the GO 1. [Ingesting the GOES-16 archive](https://github.com/zarr-developers/VirtualiZarr/blob/main/examples/V2/goes-16-ingest.ipynb) - End-to-end notebook building the virtual store from the netCDF archive. +### ITS_LIVE glacier-velocity mosaic + +1. [Mosaicking ITS_LIVE granules into a virtual cube](https://github.com/zarr-developers/VirtualiZarr/blob/main/examples/V2/its_live.ipynb) - Aligns granules that share a global grid but cover different regions using native xarray `concat(..., join="outer")`, stacks them along `time`, and writes the sparse virtual cube to Icechunk. + ## V1 Examples !!! note diff --git a/docs/how_to/grid_alignment.png b/docs/how_to/grid_alignment.png new file mode 100644 index 0000000000000000000000000000000000000000..afee919e0c93af78037e886833aa3712190cfe9c GIT binary patch literal 38148 zcmZ^L2Rzl^|G!yEMuZZ%lo6rK>?@J2knK{oYh`AWPZ!BYZR`1 z-T%3ve1DI}|4|PQ?l|XlUg!0CzV{-+iKV;_cN{*3wK1w^A{&(mmP9e&IqdbCR9q zy2R^lF;XNDKRkkyKPb=pJ5CQt37o?xqJ-j|_~CYj(lpSi;r6K$ui-x8&q72}@)G)F zooyua|12ej|7>Y!m|d@9gYJL66yyossFM2XT;z=J_`d|aIS*+GeRVYhLdl4)D0^~$ zZLojsqH~W9QNmTI(n;D>R1fH$*F@0wN zZFDS;2oJZmqf5aYv6M2jsXtp|mI=EZSY|h_S5UhmdP&$WVQHxNYI+(;iOe(${}H)S z3Vs3S|1jL)LleTtm)m=HzG~dk`n%i4my7|-Ft77@5%IC5$E~)-M5K`r8VB1KfZ~eV7#|1xXT}c zI81tPeNR`2ocZuR^lX2oTdB))h(3+{Ke2c~oAMMe57I!?`nIi_B3~MwA_3YQl7gDSslNAKl=RQxh9pgtGt1! z+%f+sp(!>546pGBI%%F_FAyA ze9;%sJEZOL!fGpxSf6ZIeU1BeECMCimYN2Hl~c~=?~Q{~7wP&fQU*0NjoLTba^h~G z%yx9}Mj$_P3;Dg*CMN|751$57tq9@g*HBkiSLPi0mF#;kY||??oU4^jg|VxcOBQ!8 zF}1LW_x`9%QWEfhh2q4qNgz{v_;GvDey7^B9P-=Z4o2cj`->jQlsbGQrjth$7So|Ne|NQobs;HbHxyxe3K;Pd*T|FW_D~rzhGx!!WllR7( z;jJHmOs!vnsb_Nw8zQc2WXi`2*`A%{wIh^veZpsR?;l@rzd7XlokHxI{wgaI=sV3- zHHfDl-SJI45(rbvLRK-I)?4VV2YXb1XkUC8CYc za`JU@p{h{l$~0Qjdo`(1fqymU5J;d;2*9!ND|k?hR`u|Fh3J=JJiVjTw^~BWm}i;) z7i_-BLkX?P@&p%3PW994zZ&r^bms}7xBqH2o?>dw|OW7p>TK5 zz3~j5Z$vF!!5OYd2l-z@rzr8kqO)WLfZ1HILg6{?Yb~KW*HiG#Xr^#t{FgiO=kQU% zKLELIc`iCt)X(6kbCjm~RPE=3>hM_q=*AH^L5Ya+rbM;W&QWBL5)G$XzM=E4I-hd& z#GTzXg0AYIiX4QA5bxNeD_B1=;1eY+zS{po8Z${a%`qQ2bfJWB754ymR$XZB%saaJXh^LkUB#bMeyPR7$Q*LTKjy5C z+PTVSa7OZi&z}-EL<=1JNJ>6G_<8Y=acLaTp2K5*DfOc=et10dPb-2d@Mf<9aV8ch8jKRIa$QGNeG(x9`+j(G?s#86 z=(DrkwnUU1icKDmt7r^R4i1d8E1QTXGy{m z))sQvmvtuu1VUMYclt}R-$e1?qaQ?w=Uv=KHtV1MKp&HIw^0)3P3aNZ7*@P^j_P(W zBrmS2YF7y9IC{T_|L- z@4L=Kv6mIgX)k#qfR!<-uo+HMxxfxkFwue5^2q)jof@d&p>}=#i0&+>hBarn744x> z?BBnC=QC@)Q2qPcWi#~RK-$MfneFr#$msQs_lzeQ2CmYTMCW69WYWb37nIg;Ib1ts zPip^HfT@|;_p4#zEM)WoH)XzE7Ph-#J6b_D6O4Yg_iRb$U`5FlKBPov!#pEqdI~QPu?gVfgChI;}{LQt$WQPuPahB zF6K|Z_xbr+`0RiKGSRKosY`-AOmg=_-Fl}O-{UVN)CNXw<1TzY`%l-$;f%CC{SWR1 z`v`lsgfsWPzegeR^jAvVuTP|zbta)A-`Y(?m|H?KBwb#Tkep<+ao{oFpW6qQi`z=P_1ErQaZle^=({JmJw>>F;kQ$#$h`@lJ=@_* z`(*L$selrTt}E*kE2DO`9SlsUGt#bkU3loLlrQiCPxzITul&$o9DwEfG39|29r=jP zd6;Y3AmYL;xOem|LZpm zmsU%eRo|DhO+bjuA`Jsv#tYD18$I{ZFnnSRHpn-o&ckbHTAXZTNz9Q*(xjAFk0 zl6^@VJ5>co6s&!vGirIO4ZD6Nj2*Ieu(zs|rJC5Dh^%%(hMZ~{R1Issx;5$BXfsyz zV8U}%Vqq{pejZf~^ToJ`tW^D+0k6?`KaYNEgrwor?2dg86!*RAX<$A4=g-{dI~IIO z>>rtbKFVmI7AK&vfp10!dzy-Kr5rz)1fSi7IZwoH98Z0Kgvqyax;4((V?g%O+^Pv<8Ec^7%iYv6rNATfCz}3{ zySKF(U2M{33dlK~I{ z9N(Y{kh&9KN|uSP?V%TNTb(d;j;U<$$ThxR@3|%9Rs-+Fa) zmgYx`EzxkOIJdxAoFA42eJ}UgbaQEVtzScpW*FZ)D0t;IDzEfzV8Fsw`No;l8tyN5 zT@URU6%1#u=?%DEj2Hgw70ploSX7>;cDivogzK&mBBPN`XOy(FEt;2?LSiEvhlsd1 z!6L~$Q~dWAYR6ww(oU87gqzQrs~#mnM#(IKg>X#aq5Adoz@tjxqq)9Y3Ai=h-G#CS z#9A9~t?8vyY7ARgCG8bdFl*TO@9CGk ztp$d)?OY-owIivU_k63n18c&2uUQ5NtUdlJ*aZGj@*@0vhVrnki$xAxDFZH9*SzVs zJA@~1@LP07tu|s|XAgFO7|g!n@4RWm*eCXj-e%Hu)K(o}#xq^RO4ab?8V$KWfz;Mn z2KSja)9aJN^N7M^1D~B)KGPP;zF6Qb$G8_!qv%LG3b8jwV>k9}T?Z`&XjajBn6kOz z#X^`y@dq(lmmEQKW{Buru?3fm>*ZZLk1Z6fCqE7D<~2MhltAhhUdM>#0h1eJ*NAcL zyHOvXFVobaxN`Si*+L>p#Az-@>{mG)(K-I};F5I3UKg_j&Y&7x{RPyn^=FvW70W_2 zW#lvLuTy+cB`J+MPwQdpVWo&ncSpyPe6%hm4OYm#;_>~x_nT+Oxd57?ii`=J{slBP zF%@Rt-&J^!(2G(eR4bHfso-a=g!AA|AzP2UH*Q$SUSE)qv^bD!Ay>$)Tk@rDkjvil zvXBjj2U2f+X8a088__V2*vsL|?^=_M!lD=p$FJ>0cRJ{z#wW+58lTF#hF*er`aU%x z#}p4^7}HLZ+|_{wdZRq~1g(17-f>eK!b;6w7lcS&R55D)-O5!Mv#ixL4$Mty@)-h> z|5%_fAi_;M1iSlh0S#?+?&}GYLnP(h1JaxR15VfAO;#mbm|!o4Kv&zhwx4<1qBlD* zoJLoK3ziPl;!MXUdzfa%krZpav%e-1%@)5Utv0v4>)$NdYW9uJPPi+2+q2%M<##4l zCiI|Wc4g43`;9z%a>FUsA`i56PwG_M2A=n729Iy-o zbZz~-wpTe@ve!po;jLLid*5VxI~r8%j)f_ew_~lh2U2g7fd`uXB&s$7|ix?sD$F!Xvmvm7)j2 zV52Dz=0^Y%IAZkVm9?{SICmh$d>GOurKu$!eec70bR?NmJr`a80CWjUB7EfJIczS!EAR3uONeyQkE(G|jL_r^e57 zx-yTef^IbMz;GiVAwXaFLi|0qr<~f6iuqk{M7LDKbC9sC`1U8}{y@%%xTk(Ut@30n zqrG$<3Xvw7XCol_65M3;CR;{Svvjpo9s9icx%$BRF0J*W=O+3Jg7sIr7o*wXUwlGns&LFbQo$Z+t2oSJok87ypx;Y9^CJ} z=$YMLpkM5Mb~z|qE!EMa?xRNjJ57a}b_2DKLSMbAes<-(XHz^oh0)z+I*z#^W84vbc(y{C}G(~6C)e?q3co8rYLIUO-w;KhG z!#|Sy8TS+W50Ev`n;z&4MgT5V6+0Csh@HXe}>-P1F(W@lJB2X|J*CyY@C` zaV|l#;ZCUGthM7)Z1Oj?(b!Uc3y-JRE2W0XDa+Zq+F$H9np?!u4aVPhtR*McUkfbV zif^Y6hySuAj!g=AIFo|)vhmDx>;62bB7N{Zm9bH2Gjrf;W9FXY(+tm6Ru@@W)_ zT86c`p^AXQbH%}eA|wrcvQ4N&RrPl2JR2ePaZrlq{~uLKX)xDR6+}Ov+}o)!3c6d} zr17y$yXBH%l91Y>`1X_7t}%ddv3@09@8^1rzbIPWbj)YCHBYLsQV@o%A+!8&++aM2 zl8_p0=PF-ea2^e5(`h&;Hb;EyXp9m0#|3;M1|H~o+6};ly{^KE@Z~Pd&%%u*ng^^y zjpBV#_6@>@1qlqhCZ(FUD>7n?dwox9uUf~SV^-y4@J`c|L=$B9mqap1MsbTZkMsK0 zdAm)KJIH@!6pe`-_eDxtd`!yDvR2V}lXb`Jy7pJS0@DiJNJ1{hD>y+bqcFn!Em)!MBy>A6X-ayC-OQip1L zlcPE+1u3qPRg_Q5h?P3iI_CH+5UL{&C}x718zSUPgp|HgN8UC1Sza5Y{pPJ!l_a;} zIzwx{XI`7CLZGYNEtTP29fY6ZASNr`HgIc6VBMHqG{+eAb3DDOY&pA;SU<{khun5& zp#NcmvEh1$~(~jK(rT}pKCVG95kx0Cam%Yc4zJDN?a-R**5N>VC0Iv z%<$00HK)IBZsE@NvkofEKU7<_!*db~X)&d~zcmw|N3{xDo3B~e>zL1dRn2Qy*?B*7 z88?hWDeqA5=4w+&th&KSqp3%c>};1d50VeDI)U2@TtX|Qmq;}gy7ov!2)pMD z2<^4;iIlW%Gs>rPBt8{!?ZUTvkMk7vb>%c)qQ-nBCl@e2Ay}f{=~AGp=236cPT!8 zz1=UQLyZgmgdyyI7zvrGfq1UIqA=Yfz0QWYwRpwp>;~i4!HRmRgHfAOHX5gZd2fEH zw|RkG?`*paZp7PLu>@{mtV%SO9qZ?1y2zP^;4tm-*`$g&7uc@z(l`?Zzx$`yC5An} zd5yO6x{gW#AvcS!*sR^JsgtpV9us6~YRjFlTIeQ5pvOqD2OIB)@I~zZ@~+zNCbM^e}^lIZth9l#SPG5{7;}B9@(}nzd(e7)Wle;D8 zH%_L_T*5KTlY2e~yM?cA8RU44JkNAxluaqL>eFtp^{k0rirw%vuT^z694$~lfXr&9 zl$5WoHfcuUPbR)m&T`aT;|oW*%Q0liC-WMRN5l!#rbu!Z+TZdDzT|mHSNF_)ANBf= z1-@J0%FMGJZk%Z}G0dLTP4h2%7rRVBw8gs*BVXS6EIy@NsZhP9zO`^&9i0G|-TyfK zu>;aR>adG?!w^|tQWJ+*3yB?qUGN#X1#2B0&Y=2vSi7A$XrUb)uHsik?h;I!c~FI6Wch<=!3FoSiiJRa55xJ-JDaNU#zF{r_fPGy zC(;NYfC@cR0*49rpQ9esmVWE3m*@4BjEjeJ3cxYE*I^?wLEE>!Mw9+^8x@&HzDMiB z!mH(Wo-!K)56I;5ZtghS)wi~_(zVxBZ4UTarD`fPu*9&0v}<5iw-h>1Nw@Y~y*FN& zl8GyA?hsm*GIyn`kMgLywP;%X+T~o28j}4^HB5Hkvf9%xY z0vRg(rd!gkWoiMeq}=N5hSDMUQoE8ttE zcZT6<*(*zvIzNWA)bpFU)%Q9dsAJ!d<@&(A*E6LnEiKLP+HveH zg(IP&)g215b7sdCx++7z&7(#9SShu_y*p^?XIJ72qiIzx>(9=_Gl@UFH-gK?yi zEM1Ol@k9y-Lg|9fi90EX^345%|IW z%&(fl9qN~}NzD-Vy25V;Z40I_O?^Qc8F~PoafZjnUYAVn8jO1L<`LYxVeH^RVHd4i zeZk{Lj4wsK58$~m9`EfE<*Xs-Tkoi=T+x9jpm)5mM^^ymq7GEE@526$bxgz+9 z)2K|hmOMfYym!7v^#9hZPTL|@FSKLaF&p1yA)U@B zvz~|0ctgi<8H64ML=h1YUr^|gxQNI|+x)KhQK?1r$TZJ^RB=(LY$6O>KHK0?Z&pM8 zxW@=1GTWGUE!&&)_om3?VANPaH7?qJKQW(dQJ=I?!y=*qnU>xRm2z$U@O!mw+Kh`( z0%Te17$m)JPqJKke4Gm?g7x1)J9)kwxB7A3)}~NnPrklK3L$G~Vsli+!6!;%gMP#b zNC8wwdNbr9+;p?<`P$yx4UJbk@ZPs0<<_m6QU~F@>?*h18@A^>hwE6CR7D$}udw|U zL-1|C*94k?*PCa7WcjXRDAv$3RwY*Jw0(MpgQi;hk)$<=Tn{<5h1_X*blJG$1mCTT zvhQsD`rib+F;NVJGMeBK@vz%Qqw1;52^fd&9b3I5NE|yl)W;D!UfJ?CdLu$l5Ux`j zdnA%aT0%~nVny;FFAz^g8y64Z(VxGG;uF0N73i;uv7RaF`c3eFtb}pWlRoHJzYgC9 z)%>&HK}$jB6XTHo^G-5B#byX(RLhR@FV7D3Imp7l``OadcPaM&XCsoR2N1hwOMIdZ zw~v(gFAZk2kQN^PN7=;T|AUhtK`ITpBdBNY(v7Pe*4qAlh2*X}K%j`YC>jxb#XlAN zLzM!0WxxRvte+y~B2BoU%l|_JsIC$XfBjy>h`)8F`|;l(1E^w9F~E~oYIp-z!!HtS zQ2j$||1Xq~7S&hl%Dm^HDYx+}PH>n~8uU7lJXXJeLm&nG)_+tuB!gauz0OPAxDx$i zZs>7813;ax0O{xmKb!yAT=;DF;VEV>%jWbW;7Oow# z?JH@~=_C4B*rDN!B1*^S_is6SpZ@Se1)z1qb^Q8!oDkz~j?xof?+EjHz;d#+MSn}1 z4q15le_9EGaG0vRlMQOGhgT`jfbW63_iN9<8!QQCyZ^-Fegxhg4uFCH*x!2<@$B#m zf36^?3=dSQ{YHrxs7)<8jkp_rPY%WLJ-(l{`NreDZxTN;-*t3Sb_CQrACZ`vn}>8E zO9jnLO^5W7c`JSUQKerEc|YsuG^eR1V@J|pM*W3c7?BLb!@0VL6$eT1w$+5M`zxKR zHLm*wwTir45i-PJANF!M&x9m$G%mvAc%R4-VYhiBkvO1$Z!X9yDpD09;BZA%)rxXS zdpqn01pH;~HulRWQCb%d41yNz@!rR9q3b%6pge|~m zSMa-A(j0R5n6bC;&nn^pK_I&Q-T7~K2;htkY2DhywS~PFER7GIXE`c3;w=-&dm z#jYB8@v%+lKD0cJd5RH^CzuCK(eY1?KH!|hmX+P#6>=(mg>imy&gZ98FJ_KzCiyd4 z%ZkL_OYT6%GaOSTJk+MOOaJv4oDk$`-KoCsPG(a|OkbVo7W5g;To|Bvwv?thbRty4 zBX~glbboV#Z($jK?Q4+d%=jbWwtB=?_~JvZ(~Esj&mXO1{w}nX4v~H`&uoO>^k|*f zfRf`S^ymA<*59u$E5`o}ds8Zw0NFA?7&bf|kLT@%)mgv)bsAszEN}r_2Z#PAfOtS> z>ESFPQm+yTmALFPNvmUyk1q^3P`hA5W670!LQ%iGY-+#^m+VJ3Ch<&6SDpIZE4cf{ z7^elUpKKR-2L6Engfgn>$LBz^?~AV;M>T$${d~tb7vte#-BO*EZ2B%1G9`rn%%(E& z)`-3wKKY3Xc6pOjplW!hcc4PApiH_|*Xh}3qrIkUVM@xSHA%yY2k5bU^;)@yrJY3h zl{gZYU-e9E)1+7&6A67rvNTbuHfTFWU*%-$vX^>5*V^sL6G-k)ix-Tc?@;+fCr+7hSH=2KTByb`^YHLLpr3|j>Uvm}RnQqOa z3*%5tam+ZYUhXIxnWLS5tyxI)&PMIfs+4bhnbgAV9x1LtyNF83gq)YxhzwW(Xy)|E zngGwxpM>N-(V>)+(|ph@er7O=YKNlCc@J2N+6f5kfJJa^i%m-e+<+30bju#@CeU4R zfwo}3L(eqS4-R14&?(1~1GftF9PtNw9Hum#^8Xx2(Pp7BciLJs$7i{TkCdcJ)(6s^58G`OtG~ zvLz_PHJNe9rC?1=($SMlWa)Es;@6O9PgTi%bJvEe2R)@G4K7-g4guhcmpiF!k6b?g zptJWcf~+UZcq_^KaYmZn`ZVTxZ^`f*+2zNV2||HQa@|-+eR6c?;l>2cpl^kS!>D1V z`@P3NtD{GJX7AQ&;6L#K>n@*FpY{u7Yf&+7$o1SoWDh49vbYT7aHTK(J)cF9>@-UnWmmp$a?sV|3 zhYDEFIH6ml=3Wid6q$r>JpW}W(G zsMYHI&P$j{;Q<|u8$Oi?oGTff{vSLO)DO~KJHC!0z4SF4&n@;#s@rdLIjQf>o!N?x zCBGUTM4+Ic04iy=DxeRdVi4J}%mn^;LWfQMV=@YHCd2-Z$vE;9^?kc$-8$#74cXt~ z9{;YbhpBlOzG>y8O#I0RdV;b%qyCR$7q$`JG;*^2^Q;L#AF15Ck zfDO1btqS;Uo)s1Hzflz~l!1s%bXs29mAQ11;ehLp0ICFj{U%*Zm4YkxfJKgeG<^rq zp(-g>)}caB_SZ)L`2giJ8y*6pX_z1Rp|1S7JU>4}@C8=|4bVh)f!3JxulRt}7aQ;n zEZCz~fC`!taS;FA1JB<9=;u*a62z6DA5=Gv=$kSIfO$HZAr5CRs+9Z&|7}HS>O}Sj zdA~H`NgPK_L^f2AmWnfR)`#SCDEiqRtwd=A==$v!?i$2kfhPxcF}cE_ebTf7rbQ(5Rw0R*!)Ngft7` z@4C+&?b1SHj61}aw?MEB5>-VtBK%GbUB7-Ev_n$&%tox0P#iOa@*|G_B!+-0Xa3M2 zLVU7{Gi5)g?)DNNE=Bdw&j&O-nE;q(^t86rNhby&%QXoJiTgMFUHzUNeG0g6j|ndP z2TdcTG-cnV2{*eOA{_A1nsYvh@I#- zRakz{{qKXpY*62yIAu=iWb;0w*Q5{W-rVt$W>S(n78gJ`m7HJbiyPpo`7WLQ^3w0r zh^Uvr-uutr&;=Ddk^duFI|Qn>%+eJaUFVidWT8y|Wy&-Q4h@NUGP zxn@cdMd+#}yogcLobP{?ftBR;oT}^G{=oEQtjfWKS@f%m{n_5>AP_^gyXPuZEe!*= z>o-e^(R#N9=TIykAK#bG-MpF)2#c4s2jfA9px||tfA*$;Y`cmsX~R;XPk zlsPkFs*ft&m8*vV=pi_So?Yz3{X*eI#9LgwmlK4J01yTyAMPab?>0JpYt_XSd*wSX zI4FlND__yB2$B0gSo%dCuwU+1ur7D~iGEw0((8ja8oGGvNwrd9yntLyukVBmxBd0! z|EM&i5eT0T#ULdK$$WsKz=Bwr$tBCv^j`>$^$jP^w_!jMerE-p`lbkF(o9o8v3F}d zzY|RMPYy^R_o04}T0D|mk{O!rEemQ(B+ZjHz9S2GxRVf2km4s1Qw{+Ux7k^K&Ez1v!`2!)G{T}ZT_}fsfoAT2PS-WlI~N-;`#Ud z91ogM=#W3M8L&ddUy5r`4a&fCIpQ`{832;WC)Er)+F!ebwgfi9$kLaz~V5CVzN0VIa1@ff~L^gA_LoSPVTY!Cb$h;NhQ6SZk?L4JOf(3Iu%9hUo?vVTzq1bw)kNgz04r{*;;VF~ChQNNxnTj}Y=QCj8_ zojFE7U|W3#!!{rHvhlfi*ztG)9KJTNs1An6PLFVC!kpG;i{ii{%D6G~sf8KC;+Nt=o z|MvNzM&HI)jAB|I=tWkQIDt39lU}9wk#eC7*SYn|!ZlPyWJTRpZYaIGaXl*VuO|5E z6H^JaxS|%)#=&L#l@f2PY~Z@waGil<5?@H?m0EPY%S!UdX@Myx+I7y%qHe(bj% zP9iD;2p5-{I<;aK=~Y5fk;@uu6dB<-~-L!1S8K*`@95nO>Yf)q~!QxhiycM$G5Rq z>*gd%V^48Y+~~7axpk)OJwi}X4n)o0USDnloq6eCK*Dk0q4sysh^jL2K1_0?{d4R* zuahgzIuw^m=JG76lzbhpImPW(J9lEgdRXC8Tjk;p8o+a%_>gRY89p>XM8lopyk(sD zL>qzc*_Ju|D>(Px>NplQf%ZI8;7PDq@#5wX-XZdI-{+>POBIBYsT9 zp+TRM5DGL)W)!quvPfVK6JHC$Da`K1f=R_|1`+kvw?C7w27Dws(qjYOCVRmGsprm3 z9Rwci`i=1xXlW*gyD7F-Hz`z_E>){Zpfz=CU2RjH_m(S9(?4zFt_(W$bGN&6wP8aT|J@CoO)+c%fp`+-;G(uo~I38Df?X0QSx-PS97X8H#he{Bm z)OIoDhfxE_f+yXvatVauTGpY?#2k#rNCPrM(yr;U*D3sO64$vC*fR1JTyIai4GIfm zyikfG%hHqqK$}P6iGY|SaUN7ZFdQDz1?%MkaZn{*{t+DcM*wr7)cc5u;1*>!RBZKg z_v3D&9Gp)&b{(2~h`J7bfFwdy0+l1@+JBoi#qE!Y2)+vN{x&i$>921~Osw-~w@v1;@VIOCMF@L`l#phZ+BSCKiX9f69>1cQ& z1N|ix&(JXguz^eM@@rL7|Ba0eXiDO{cWgu>!Cr#hZhsFcM9{BS2af(*Xi2)Mn1E;` zR6_8?2pdRE-YQIlh|OjBodUHMKQzWhLT}X)mGGBv9QDOC9Fi@ewzVI?eJh^tUNuKC z+ID2NiL#$`9icYF2A52m-@O_+?+^;+YCuab4cCyXYMUkKN&0KLg0BNc4X>+O{d?)p zZ)bp;n|PBwIrQ(%9e&_0&KDc`6lGL{CvwE{KTZ*z|83krTqL!EaN1R9_W!H~uWO+WhOWPX8_m>dx4G=nk3dicUiA77kSdnY&QW9U!V881q3(nQN z2MPyWd%n7~Yl0W;Pg{oPD@>-^+NUtuM{FGDXmQdpMSSOMiEkXgk=!efdoMEU?dFJc zb&HpZGHZ^b$# zAHT9ZsmMw6iW`B4Au3SbH52rY^56N=DO%36h#sCU@Krdh-$*9kFAPJq`Z8{_Bvci6 zJ!mZ0n3Te}0}O@AOK+C7C#OIDdRD!`O9a=y-Ed`jxK>+o&$>i2SL>k%Hovw%*YiK+ z*jpMx_TO8h3#VGo2Q-A!M@uPxf=nZt&WHIhztwZ$dHS8Tz2?-g`($#yO!JQa(paBb z+b8O7)EIE|5Z|H_5^Vn1FI=J(7JO98nY_G+b~M~6B((pJ8u9jzTfUbiRR|n#lItxE zTcHPz>u7BSwDrJQ^c-^MWk*Hh9^sM@w_12E56Wy1y;7;m@kkxB=6~lCq!$5YIXyIJ zB#gmNh>Nq%51Hs8&nxnSIuOkKOT?2%L((C&>5iopw!*cS`wQ78y`&6*Ipn3`jEV*& zwwPt0RuCZq7&fa9_*poI*wsewPCm(9J-@ZTl7u{S_eRTsQ`eQ>PsK|t7aU^+-Tb|T zB{W(=|FGz*G57II*B<|PppkJOXSwGB;TeYgRa!GDc_b%pFv*{MTjPGZjHM~I*W*9x6qne zFbRxy-OjI17%6z|Y&of=Ctzz+@wM;npP^*L+aS)zV@PBrHpyuC5hU6QqR=+`GHg zJ*(avxjwoU6<JF1s2wXCsl;YALgse z@-(ExgDDg>N$jrJdQ+)Wm#D^uv?Z>KjNyjE-NX^cjW=YhJFmzO#D?=ba*=GiQa2qP zBCU=pB`g;Yp;tATVtGV<6iaRGxaK>JB06=R6+i+ z2wE$|ov}I)-Bz908ICLBF^-fAGF`O&KC(Ea)hBQWs|j3CT%O4PE)C!eIc>K2>&1`e ztV{%X<0uX96b+cKw+10fL7x+sFK)>l?{ul z8Sk_vXuB;KBwsIrZzX36c#n^%IU(E1UL~7)%-x8re1ZL(CP#*f#{~2J$o?7Zn7P1? zlQBPjVqD4M3rK#tiYLhK-p}!qEK^HnE;9bj7paY`X?XM#+wM{66p`^Lz+>%mSgA_5 zVI$0N!CP9ml=lG}&SlKFyTFe*PMSi9TR>Vx8YDSRtH{a;lWr_Nfu$nP%&!xJ0%-NZ zJIR#7$#QHszT5bK`+IeQ^zq61_w%wzwe{ggBQNT6RIFrSI8 za+Ye!#GINRAT*fknXiK$%(@#bzx`vbZi$}PJkDnRW`RmkX#&8PRudnWOW^=!%~?TF zk*}7djO<_MAcUNFUR?l0hVQeLhDoaKLkVqjhgpe7JhJo~Kj%AfNjV)`KOkqV*CT6`tfK)po*me2Iq^VZhj_R)`nLwQxW3OS~5? z^~Wm-Y~YC^CyW@$J}YGaWUsJQ7vK@{EC`YTOwn6PxEllpq7vbyx4P0l*t)xx{wu2a zK|g^|6Dcup158{?;>GP^33K=5PK~L+$~Y=#5T7<&qL zXG=@BZc8vk%#@ED=Y=U=7%F-cP`A67%Lm3Cu5bT9f}$0xG2C3OqJ~HV!?pN_xYapvsuYLXKf~{p1g056I7b; z+@TU7JkABIN?a1M(J`*&0SyB>=3PnL^uHRw_*aSXM}q$^0V#z`z1nWNtKZ2h#)CcC z#dMh0aQ@Cm(pOAU6Q&1XU|}YJLh`0tYR8MQ@2`a?I2INbz-(FA7&nVAn0_%Et9pV9 zSBNWgNoU)ladSYepwB)}rJ_U0I@Hq;B)Ply{a9jqr~5m&dapu^raY4;YeGD+yMdbgO2HAqQo~!lRJoOvuJ>s5QE8I0HB=*Zg z#qTeL*bz6-&-j3WE4B}VIpSBs#1FQAI)b=G>1zlLZhj6I+0PfK_a#YqDjShA+%{fY znQWW{bP6lS349;>E`{A+SK0A!Gg@s$Nm2o&oDu*m%t|~z65tm=ECA9Lz4skgOvSc| zzCG&O@aW6>aqqT*!t6rnX~hwyeD8h~*YCP*uZEx{L?c^qQKvpqhC<>RNHS^_!njH4 zcm&#NB*8Sy@14CJFajQOz%Wr}6)A)`cL}lV9(gd{o#t#X^N}_$B;|3g`+YG#3-`ap^}+_&M#>zI{n$I z{V0A!G4EDQ*tfCW+H?f%cXyzconq`|B>37SPce;ni@fd`R<1pEY?3goAa+vE2z5 zMMcF_)Sa+_hM5j~Ztxt0@fxL0u~qT!`i13|j~o9=_rdJ$jF^|_8rS0{)@*Bjytu6FlYdXxlBwJV1C^F??uNo52TN^7 zg{l`a~A*TWeMA4rF8*~0%%b4kw;_*ze+|9PO{cl$amtpA0|XnWl5?O zq1Z=7&uq!61gmoh1!(b!Sgt~^hf;k|&j_&!JN3z&x$7>IUuLf8U43NgHF!{s59P1r z%E?w5ZG-c-y_)>+)e8(nV!`tY+LNV>IwiKonk3;Lt+gjsT~Lw|s>s3aEil+PZw+B| z>3`5$^~7Yo-x^9TJ;W?j}Dn zlJWxn>U;W{lsi)a6s=~_hHSXZn21kzinO`+%2>6zA^Wc8+Do1q`4{AjoZ3d-A+DXk8wyEY$JMb5)5K;x?*X9L(qPr~}4nFH5@o%K8QfH|!fR=1E!w zS2H$N&;VO4cg4U=Z^s>+|IYxDy317>cLM&Jxmy$hxB=O-4=hu%jry*~)Hp9twd8&- zK^GKF%q9ujC%Mb(7-#uPC-+L==Bv?z`etA-FHMac(WGlc_AXL1UfWW)O$xISiT#j( zWqLY=PpuqwMMBgbcdhv_SkBV< z*VrbV0BgoUTcMd)Fj<}5;<}LAC}1_3qm^&n+Ts-RRQyUdsu7hhGbpGw;SJmkJOWi* zRlhYhg)#JeH$rt@%M&HJn=%%BbzPT)~l_29jnx@K?6ZRQ7e zmIQ`zi4c&D&?2THAST+N_Ee8>ze!Hs+I9w#_?#9mEVEx?V}xpU}r(q_K7(Cu5hVbKZ@ z_Nwj;GuAXbZ#jrED74@F0kTe6qwU@& zl1QXGGn4)pekrmul9abzF;lfG9G{RWx~;GwwJML=4HKOeI;LTLf9J5hSW!^nXrD*9r^;n5T1oLm|1H(vF$HIYPkK@h@{|TZ7$d#J(*~{}hM?rZ{bNo2M>ogZoNYvv`$;5nZO*txA}b z*6y1bXp;m%S^qE&-X*2V$o-3mU024;2$l?MT>}l}0(Q}++kIi*sl&u69WNo$4Vcx= zdEBpY6=9z-c_{`>zMZ?r`7nngF3)9LYw*IE!s(Iw#gAD_l}2;bhVylbpXpje7(oMU zRX`BnWiKX027kg*n=G3MPN-dC*%l z*RJc>#|W$54uQ|+pzLWGVfixnP4nFIlZT9hxC+ zy=41(6s)I$PMKPAdnhW^znljs*LNmUc+Q7h|1#sA$4)#Oc3UnOyT}E8pL0v=6WE~i z`dv9=OmJ4y3(am;Vou) zfH+Y*N9REvVs>6g{SiLnT=qwy$x$Ag&$Mo!*p@3Hucj8=rz*zHB1NiKVRIc68zRA= z|J|27QCB?XnM)gi@p-=ab^QhI88x`U<}UH|SISA^x*OvozsHPlV>=&7hN1%@qmppz z5E*|OF739d)ULKrM&FiA#T|~ZMFC}RvvrfQ44G6^wo5;HtN}pKGiqvv8A)` z)vzW&r7DeFu{&9H*q+9WQFo;L_|?$xO1ZFTM&JEMo6%aX&3$6>$?$?G{{8jUtOLjd zK1JI^{a$;qcR_=?hqF}I4WsV&@ncUhf^3X#OSjDRts<$_#3yDDKLzZy289yAWbqvL z+P6`YbjP=?Yr?+6a9=}S-n#vo;jBn4_CyZv!dk^nu`QdvM(@2XlzF5gYJJ#(^?P-q z7>y!#oG^bg+##P-w9ufs`6xFYiP zv5SG^&k^Z+itD{z^a6-nZS8{UHOHMk--%w>_cKHB3~qS_=&&|rBcrX_{%k{~&iNr} z0S>I~Pi)yXR1BS1KqJ`FNmkR*Ny$#{w>X~It8VpGu7;t0&p@+eL|4xm;Ux1xBX&>99XPkMV3ItVoKdu;v`dm}Z*!S{2j;mfSUR4! z%}&jx91b2Gu@{|&#exy=8Q@Ftq04cP#A-b7=wV)UhYvzF!9y(~SH^4G_EsC2h6aGS zxg_drX>7PU2Ja;Uv*vte4OixxW?DoA7-Mu{Ukchny8m<}eya6GtUH^m%m9}1Q z_ehV$-Kgst&T=?twmDBC(o&J}Q1cmhdQrpn;^n?xWRm-x-5tVVX^IKi1j?X)qJTJW*P%ff^&<-$r8$ z3mVq=7JNUulvu37bj!_OU*SZ$p;*3KxC*o-Mu4=6r>n3+pl|%^4CPR1E({=$A?>WXBZezV_JdHDtb*rt*rGAZ>R0vynB^Yc}UU76?i#4zRpzy?d(74if z5|wT&`m;|j&)2Fvk1{f{cPV7=kWsc{WF^P?zV3ta{yZMP-}l=e9z8hc zzOVbfuIqKZ=JQn>4neE?^U;n!pBpZK5_}d{YJT8@p7X?Mn`#_F3|O6423_K$Go7B- ztY_Ujl)&ne@9rk+T=*&hn$%AqF@qSeH1$GXSM$1m$&I^=n34TDAzB? z^Ypu|`SL$EI=$KP1+a3jJIvW>O4c8luVUGt+@m%5?mD|({|9Z2cQ}t9ZN=S1$M1vj z^QOH$f|4@fv1H=N^Gg@B3?IDx@H?P6up;Afz(MvG#6)+QZxcu5nvQJ_p0`wMTo^de z8S@(e48fBdl+#H{;Xw}h({;{tt8FThfk1J&viPa+m_>n<#L6&ri|7O}PMR4H*e+4z zll&-xZ+Wz9s!nujTLu=-6e%*gEnf96AIyTl`*}c&ED*|oGMaEBU5hn(R#&}P>MY3UmnexzgVu@zfvwVoQ2u5=;!+O zeg}u4_$-x4PV5W+8{8`VOowiM>haJFUvM+yEI5$#=O-?${K^C=6>(eTBQFK+yb*b} z5V7)Mcgb6`X4t3jjn>d{RddtvMEA>cSlj`_Vu3TYsk$?yu>;uaJZuzi*}u@zO$Gjf=QgL-}75;~(cdpDeTGofsA zK%BFA>3!Jdl5*DB_B2g3qeBYC0+%}*Wn%Mybn}dY zwk!GBSgC2oYFYxNq3C=KOS(GZibJ;}Zb;>+U*yKRp+rUYwJXZs{hP#08E*@MplR`X zYLb~k9BpHtZC2aE!Sl*dv7d3lcRaEwNA&bW#x*s018koScZ_S{JQScbBw6T0=iiQf zs}igt^25<4-{HH<#Qd~_+Hi8Wje%Q~SG-A-)f<8Ru83kIx4CgeVog78{MUR>p#xn| z`c0Fk#ym?+xoZY>Hl67)xAJL&8QsPuoW3ws9lnC=oBipVl^5;ky04|U%5r@`cx5qX z^XopIm%POyQBi#DYXSCbbss+~%6k1MpU3r%etJDEYX7XPNRcmSxPyCn@5F4P2RpX& zeekD^P5Tv(l~m36r$-kj33PSd^d?M`YFuMk$!{k>M8w-}C?qPI3X9${o zYo5w^njRi_l9f@Vh3RlGkDhK?myt>DSbeD4xbdsx@q6hnEN5>9M9fc@MBzU z&EdRd`O2k*s$VTb1O9SRM^5sTh<#-GGg=~@7%oX>*;|x2PMt5gvB+oA`sZfh)GS?S zhLf1==G{ByXWyC~u+?6jZ7$=#oA7;gFupa;NsU_7z6zAfA$*o zeskC~=p5s@ssGd|{3S8#jp%)56eb-*b2XB)C9?KwS!(^q1v-h!aXJnMN^)Zqo%Pt| znCucta9`LpXEq69`fAPfv3b+;9`^=1`jX5)JLQ`8YRwc59bQ>q<{te+N_rsZL{n1e z(BefBw~DX~SH20Cwyr+G{k>hdpWG(A=EPmPn+&UtSf+KUL2u3sKd@}jf~%6p56H#? z5;4xf4(tJSr`giwVMu@x@WR?q)j#xzQ#rcaRPSG9w_3irZmDzN6!#G5E5Vx^%jIzk zrxk;Ue$(d1TphGk&hUO(IyXm9>Q?S+5AmF@4xZz*OVXIaQxnT z!=D`F>aaZGmQQ`K5^lVmbkS>y$`)-}8W20d#bN0p1{%_^2?g0pGZhOLqXSen1UlQ3 zDVf80XSi3|2Yl;zYWMzJa&#jI?M$C`c33>tT0WawbYVPTqv<>A#P zui3AUvsjb{_{nn5zhGE-H{tdBwB@AhjclhQAr_o(n{HhAQQst9?0;o#kip4vF?@cd zn!28PAl*7_b@kV8vS`;|jqFYW<^^IaHG94&my4Q@fHWOj(NWSHV3PbAVU*i7G*h_H zDZQSqo|-gslK4zaWiF8QWULrY(RbfEWA*8PP>4M?HF5Dw_vb!4!O%78SJgw&-o}~E zEwXgAK(6EADPPO}ZmZTv8orKO;#;+gbI2y$e}oya0LLF4m<9Y7`q$f4o|0eE7uFZs zM)G#=}H)*%hb6Ruc$uo zK8E$3T4%l}>DTPm<8FGnxVsjKt3Atx*hK*5jKl*X{!ETPZl9hoq4LG zp96I?_3o)>#hp>WrrbPd)NNK1X^hU87Ku}4u*ny6uz8+StCken^}NrztCV%1JF<~~ zp@f<{-+qQ_uB7)~>C~KKf^uZf=ilm*f0_lg>=~Q%%4;Uu+pN_w~t=H9~Tn6GP9pGl6m>)l*Zl%XQmb>K#%7Q{h@>b;0KZ=_MV z(HCeL_y6?mNDjb0(DF?JBcxZ%=PtPT;uy7j-Br)I{5ESbU^sF6fMr+g1pRl)#|BnpSoFpGSL~Q3~f59C&RJxY|B--VXDKS{MC@Ocv;p zvswTv?Ykt+d=-eO??{?)H0BKC3<~8-7o~yUhC2~W5{8StT{*<5{!YJaHQkAPrp?Z+ zH7_W>F7_0+%A=bM(u`jpgmt}BdE{;$T|Ly&V4$F-{D|a$3{zE79r3yY=8%*(1J9!RB|D_o*^&e-u*k zrR>+D?>Rq1w-YT=q$Q=OBlWBJ$`1UcJkgKg#weP~x#3 z97tr&P>DTz*bcWgovh2$6}Tz#OIG-S&Pp6q)RlYS7QBbSi3k`;%C&ZW0Bq;xs^spcQ( zI3aa<|9e#yBG(O-yrqer&%|9lHSO)!BVTV#2H7P~1PzS-4&mx*2r54!d(0Z+@GgB& zgqxfD!=qXxaTCH=ro_fbO)AwF*x65l>g=E8d;Rg!X+#{yn_x)>F4RAJ-hm)Tt&3nwIRy z8zE3CFmqXSI<1*;dr3W0A#|plMP=;G-dq}Domi~rS~2SnavjtE@|D4Qy5`AB>r?a9 z`H=^`UnM-y4xds!TX&l-7>tjV*LcwB%XfGih^kY%XI?837m~ZNv60*>9M{?Z$cv(m zN7F@iOjiB8`HH>MO0QnxY_Cs3U5<@@8`G-K*TX`+NS1KBXd8V*^#KBgzsZddN#2l? zi7|Vv9ik*W+ud!wSlZ^T#MiZUHAR0w2s<3i?V?}hfL)l@E;oPgns2LJroXH+;G>JX z&>Z}^e7do#%<4vls~g6(lwFZSgPqMpBa&T`Vv;>@tKO?TbG5e*n^@F>*&2*%=ZBW2<3$Y{jn~&d*VW3=3qTJqquyq;_Td861 z9@1y3f9c_SAugNWQZ-lgoo1oyiBD2=LjG#>T>+$xMI7J2D4Cqdoy$K4mOg_L# zZ`+``su7dw5LtIEouol>MN+?+V6~W3#Xo?LzV?A#DhU z_OO>!exuzUl)5%T%faO<=$gx`*ltj6_nNjPUvxl|gW99tOrlB{`$vSenl#5fc(qj{ z#K81gzO!uM;IQ+2gIjW#sqJ#vqj}}gOQsIr8$rW192c&1{H}|#{YP(LcEm_uxOYfj zJN?f`7w;=fj?S^|ljD;&f-2VHelN)`M7y*e#F#7J3@LIAH~C_aEWp38dA3F}qR(vq zou8}Dk+p80H09)h?)Atthwnx0MLNO2S2JcB`W1|q9$5de5|;bC5t61uiqgCA$8z$_ zd?f>0oyW|D&-_#C(SpzWmxMT-SUc94aAMz6-**389o67_HTGQM&Mh^u<=17qQU$5F zF(D+~DHUy8m;>bl9aU zy3v-=qqR44F9P+ezp`W6nQr+phA$yUQK@Xl7clba=;&F?a5=s`b(s2%C!6)Xcgv1> z)VsCS^J44t!g+9|hMO6F8=GS;b;A*Na(t91UHvvQtmU7W4yL&^Cpd%)SdPW!%rEuz zY{<@dF&hcK=jML(b%i&1{eeJ2UzWfkFW;2`E}D=*}%744d zsq%X9ed4#8iPE=oB&d|S(ywQ^7FLug&~D7Pp|4MqyAxz~F!4|Rrd)FX<*RPvo@L(m zrN;AoVp_ohs`=#YqNVx1IwvPn92}M~cj!Mv3%Taqlu8TC|F9Cq?V@x+>qd)-rHL%j z?^Pb^dvPS4U+CRgbn@Tix~+Fcc>RpDb2ZPYEwbr5i9S$5ii_pGoo zp3$T2LiG_*0Ws`jk)pTX+E?32=35xp!q><0-*Q?GG=JsV{7Lz=;i1R;H@4Juwbws3 zg>{2g*Ms)U3l1!FCda;7Kk;FfV{Ekt4x$Ezb@HP+XPV_KuP@!d9$>Utd_&pU@u}0errTzb#b)Vawr5yABL1fL+K)B~&DC&Xc)v>=u}m*dX>rdZQpuZVXumqi%uY4xm8U*t@#*i#G(4U2iT=FB zjpOH6%v~CF?5DL?t3CMb@BdMqA!P|N&unlT9#XfSt6T6^QgJU5`4%Vr)j@w+;p7I3 zZAb~^*JkxM3`d6d)^zHSu4&8qrm1(=`-GW&1E)$!N^&V(-~e~;Jjpka{Ur$mTFa)^ ziUY>@+`a9*aQl78%#-J}26E39V8=x+ zNPX!4-08WRceXlg`F2yS+dx5X`QQVb+e)tHg-=zlQp43QRf@g7`AYUE{KMQ|I%{|P z36cweS22+vWSY!LQua0J+#{YO0-RyT_X*Ls(G_0w54*71H5;@M`RZ_~jYXnFtz4*E zrSjtHaMR5Wz^xz@vFI=4f88&5y`&cfmE!?f*E(LHwP9601Q=|Rfy?lvk{~V zJ{{MK?Mh+qbojiaTuhQ26xo>WRYt5mk+3Qn;t;MqE6gLj$+euTRT*S7VtSp^h`%D) z#>c*8S=Y4RlFs}<$-P@l?n0pZ^l{8*R4RsFJr?Nk?7DsZwC7dgy>});J)0sybzER(m#~0Rd9*gUnQPhjH ze0mc5-E2`Ms_4U$qmAiOj%}pwe*1JIIP+g>-=BYild-~f=r}L0S6%9;QV{QH=X>w0 zwXzx=SD8}Y<!*WTyRaPbWn`2a?tumMn8CQtk8=+9^i~Zy^Qb+8uGV#m4FLrkB zh6?r~3cKI9SipoaQo9%8h;(_WJUH8h_%&bp1ZBOfK96d5ALCH0wcbmEk(V8=o|_lR z-7+I3VwB+hxpmx`z`5I|D)o@A^U5eKmj1+vFZ3t#`w6|bG%z@4vmRCGp3*k-q;M41 zR{Sm2*O8a_oct|9=YQ=7{H#!gKQckT-(<8tPZPofml8KqQy)BP$kwiq=*w5em{=GKm*6RA^q6?LMGQAFpkWF(h^gXmd5 zwvZ=@O{c;o#+RM9=JVm;a>S-UZH@d-u8 zZR;6oqpz**C|B-U!gMpS>O%;vbCan$m{@)3WPJ@)HV(hH>_)^XzwpX@YO z->#&!m)hyv!`0;&jg&A>z4pz%|4p-XqafnW%#E^}4XHOc*m{CTE{cZ!*6gYbWH?3~ zmOxk##p3oPq-Pzs);<4ai92zSy2GVJSi@p|;BD8kqmZSldJ=`0(uXy)6c_X0UBa+` zuZoHR6GGJ8GVxU}@xEt@)Kz(@5=jzTw-61Vr|OByx#S)UQ=X^YO50P8Y1dUdS^FAT;ltFo|JKf0<)UwqbZ1DgIMGM9$K~yv=O64W zUUah0GJb^aQD+OIywJEjP})UFYnJkc11#tw>0z}0dt`QcM$DsjRCnS6YFr@@+a4#WfORUZ>3_?w@5W!@DJmnK=+o_nkylXtxqWm}j>P%`@J2VP8CL!TT(ruIU?8azp;4 zA(-$&(c(S6=)EYXJyQt@r%qHH-U77+sCy3Ki=h-C3+Qa%v;AZRDrhe+ZHLUYNA@``CTz`ccn=#Ung$<~mN>Nd2H8M=y=RBgR_ z3}4`~RkepNZ&CSI-qPfYRp3btxoP7*4MpQtogSO@Pd-<@M%^LphS<4(yswDDrB{p$ zmm%lofj*lF4eI$Rzh$oscS>z9Z?P&I)9cGB9=fk1&QRI=m3nNtBW(o|AWhr4-)StgFheva z`~x|r>R>Zd5&bPJgW0-Ngs^yW^0qxmC4{miBBa|F0JRd3?0odBdyh$SErcQrpc&NZ zYvT-p(M@a;he)m#kSLV=Ib@dHwJD!W(!00Qyl5x;olHi@lbniOiFt~7b=Y4eDBPk* z{?dyhC%SIgybOJHlrO0Yx9Adw<7K@TcALuS(IFm*ZL-RiM`-}H*mweYM?&AN9*N<` z<6BcTc#iVk_6+raY+EV2RF~-MS)u}#(fTf1Ri1?ulW{E6yOpb|%gkBt4sg~qe*Jo$ z^7$W#OYYHv2#((#Qi>X2@$BjK)*kPy=Cn^T+#o0pm&}~KRC?9!RiMp-uA{y9WwN6N zU=9(-W%8vl7egjQHl)WM6?4pl^j-OWXDs&$r1y>Vn7*fdRZacG@Y<=}13*-~Mmm7} zby|o6t^ZoNyW0dX>R-?cC?Gg_YJGgluZ$j=iTNACdLSjgR4#b`pY8!(xZ6Y@DYBPp zKSP`%5P}XhFOKpJ$9gfofdpj|wd4k;6bj5A>Sh9Mrl@3b7JO#gE{^~enP^R3m>&28 z&uC=Sf4=$obZ12WNvfylegbXDsn#yLCs+6H?w{d_rXhOf8}CU^DceeFL_V%#ilxZQ z;}sXzF#?qOzj0YTgEXj%%0vxhTFkM0NeKx=xahDjdW4g->toOs$#43_0-E!uG;>iT z`;psxSHz0?1O5p@PqZ-XX?!yt!F}V)nRY1!H&PukFH&SKp^;Q*{6au}z)E5kovga< ziS8x_Lc^64b4PB6`=0yj%z-$VGI*$lhGkk`{yE%C8e#I2A}@U`^y{uMlObjKNP0gM zQQ?X*zO`RTZ={Zq{dX5-gy=CD8mj}p@;pXtywk9a*>4a=x!0Yr*$o~JJi+6CHQw7^ z@E#bxwj*_yVCBxCrnZNK5eC?ZZ@j;}#&@3ErAc+`;<2RQwPQjpZdD->5=q0|@co-) z<u8QEf0ZwdPF zj?l|yAU&%jHW&L(poVHn6Lx$7&ZNEf6`+S!Zp{`Z4Re`X^)NC{PlRqM$gD<;I=z=> ze`cw!rF?hy76k}>oV@uo@|5n#>m0k)N=-+;H8FxZv$v^lrxMjpqv&o~VVuD6tXH9+ z7<|^wCEpl)iBP?qs{A1AdS_yyo9jzc#lVWe&~bnMBE&xKGVkcUD9_YaD^Uu?ca&`5&?$-#%LM5M2p z(mTlNJig?YfnOj2Co3#KBE#gysz^q@$^~?aGo5JXG~gIcp2-Bc=QQ`3E)iVhaTlxQ zUecc~;(lW{%&jQ*Gdh(R}>6D2d_)Maw(wg4+IIC4c4oREWS-@{9}& zG_ZDpEL_M#hmJoLX>!hS7^8d%-7G!<%iyibNvQrT4BgyAgQf4#3!AE7_Rvf&>M%WWRu`xZ{uz1|tRk@@6XiNDA1Y$^F4VNb z?K7=LLC(bwcn}xf0=vr0RzootlWiEy-_K^BB(QK^7JqfEWLBL#K9k6>sWdfgs|Y}l zBj2WVu~b9K`Y;i^7^Fp$KV%rnM@zUN-B`{;N12IlS@(uHV*iS>ib~YEz1iupmX-A^ zujRZvOuS_@V*_IC?vsy%TpsQa5s-~bQ~la3BJjle%BgvgKeNG3ci%RR@!7ED`{$jY zlE7`v z^S%dB#)*nVVAbO`TFCk;?d4}ppvA&J>5IoHQJug*;c?c?ZPVAzbP#ER-2dI583u9I zm|A%rWk3aj9veAzo~AnXf7VKh0yI`0ZAupn~%J6f5fU^PtuP3yImCm z3t5ZOb~em`@DKfM(-?D`j_Y3N1S0moZwVdq0%Vv0(XbW%yb&cz0jkNCu(v{FvS=f9 zI=&3h>K(!KR=5aO(mYbLR(X4ia))s|6a$r^_%r4 zLWP^bC&&|zj*v}t{feiC{6oU10=aK0d*y-c!@b5fKp){;5zMCj>C#oG!791o8 zrBow=m41SVkh7vQ55PB*FLr9uPonTRg9-#g<)NAR7@f7Px6V12eq`EttP+9)$P<}C zL({y-2=YEwfJlW}Itn{v_1Tk{@N7+nGEIyU93SuDN^{PJ6B<=S(AW}XVB!e{!+Vb5 zsm1uWw1qB4hr;yNRaIBcomJ(nFbr`liGG0F!NDmXjxi@IkOJ9*G!^|VOK?-H`Q%w6d}9;0^)su^@3oOmlySIy1YC8_+FmJh>YX=U`X~pYSE@hgGF)i zJ^}Nj*j;O^;1aX)F{mR=W(%_TyK}bw{+LZQC?3WcSG`r{caHX^IFB0Z%jNE{8%<4_ zJDlE)7gt>G5`Tkz z@Ot>>_AP*|}Ud6^3<4VuI8U>(RGVQ|8;XR z{dMqD5>qTXmmxuMo>oKhgHY{TF} zGYy?yNy=0~10ZTbo2$?ngzi&|Hy5~oDO_-lBZ9e{Y zim}Mk_g8j4HT47|)9&?-Sc3370T7ZfN02x9he*|Jm?-VX4?s^ydvdS+RlSe1qHU<4 zzbdGp!cBWS>VVwwc#K1@HA;7uvd73XCh(J1{JI|Nrdd=K#ox zylfOzhVVTeux2{zCS9SQ%@IhCiUvsejZFlG0g?RZQNTdf0Zs^K z_MmH9a*UPdrF79PVHG@-C+q-v0c(n{N%3Q@Hw9AZ^yB|WB#E?lNjaD5&h>^9;8Xx|%}Xr2 zA28B?)?Oy-*R#0oF`0LHe@dZi3km?;tPuT)p^Z@uDOn%9La6lwT}AYJrvFTvq`P*| z`l1`Nl9GoG1w*~>O$M5s#VZSomTvkZo(GKaB-m#`eUnczN7a9Mx|a&zg z5?nq>n?kB?P41M$L3B-jL=f)m)?$&7gFP*hvCZuSaGYSeCP#6eQE1Vzd=Yj2jtq|g z%}w(W>F&ZMhn)Hyt^$xG{1_NCwtJ#8v!d!2_V1O<8gfcH!0lfYTBv-ZF*JkEixt5| z!+ktgkrW>fiizzVMTj|lbbusCJEFCgi^*fM+c3jFQBt|}=$Y$;^IJ$UssxyEKUMv$ zV)%*)l#BhepWfld{rn3BLCGN_hSQDwHj!N6E0M?RN7Z%jALMzC2R)DpBZ*U*V~;ww z#0^zm#k)dDB46{scAc*pDaM_tT^XEmOaq-<*whkt9g!|BJ}Ks}C>li8okH#@9c3?q z7m-x2P)>S%b{{9kMzMG6_=-t7HR{6-f%gDYMHABn4~1R?*%WW2SiAKizp+Ary_6>J zNDkVKyN8W5$S(Ht8Ka5%geQ`&CV=QbS71>5klJ+%s72E_wacViOOdL(;C|M)s#V?H zNm9$NUQLu_2~3j5ET`WLbY=K4gnR5u%Rx3g5ephpRLBrAsxX`Nmhv%Zj2R(*>e;k` z+fy4pI#yVh&Y)SDS`vVP^)+__&4YAu-f75*ND^i-Okyas1mO%hB;CMtC3F}sAKBR} zoKM18!1hx!xe5O$4^C&g$OWZopBHLQE>n3Kq^;{%uBzZPo`{Z;m2^GbP5PU(Qxh+j zF7px=9TR}1nbLBNcWXiYCQGYdn9N`a&p+h$#xcjgDd6GTb*Iuq>yAmxrF0#B&c91f z`Jk1H=+!utzkLLDvcQmc`*9!4txm;Yg0Bv$j+jm?%DNSo7KLfxG*meN=O}(v5ug&u zi2^aqqU@~ax8KMW%thz}H%F;QOyBi0Kp8UzA7Coty zOtmdvo4X#2xY=p=3nm9YYb|lRc6&2uDb)a~pZ8$7c8;xXSi<{`WN-0pJ;R9>MZaAv z*X6wN?s}ru^M?02n<;2r@g&2(r>94ac6={ewyFk~joRNH*;_;R`~+8HaDVu{MNDg# zIfOh`4D8=`Xm342XXKRJzB~Ny8#}b_U=dV5;iuvGAI0UFBuu~`(d0!kOg6cd7%Jo7 z|F`P-Ir{(7r^4qyswa4L7D|x%Qub?}FUf}Kn!B#rVTWi%)fVi57HqWi5+TrskCgRu zojD_4qWeDm%!Lcxow$<1^?0OZH@`n!AP=5>u<$|fMMSC=?XtQ-jXW4b zL-Y~Km?`95-rE@%UL(PL+pC)M2wB0d*2X(wl&zIIc$QF0(XT;ZXNPXz3t9jRF8ML{ zvKqhptCVeB@|g$53YM=%tk>>gAzy!dg33tB{s`*$e`vbMJK4&;)VzOA|F3qmeK8l< z3I&p%6CM1H|Ig)KBjmCdEQxU5UE6IK>`FfzBWniu_O46*O*kC1Ak^B?hVj=uhfLT} zaatBXN7UjG{7cz-Gmc>S`Y!Wb*#vQWgmC`f_GO0VAE2c3UL%8~oy>DQVXuqySC)(Zdu*44m= zSR^)cW&c8e5fg7qq1rMXsZb5QZ1VXM{nBT`dc*P+p>snj4?-UQUJpj8L2y6-AUV1^ zw@73CNC%s)!wz$s&+&9&9D*Fbgh!5V)Ol|*a6lK7(tFhuNIU9_#WK)*8QB(AX#@wh zmm@})4nR#vUH0N+ha%J!u|dIAi5O(BLp$I~zvpSOmcJX?WV}4PHKGL)UOG?nQm&we z3GmzDB+<7W73nns!eWGg3etl5i!TEwrPwde5`YQJcN;|$_feSM{FFb`R@Im=(4B)% zcFi3VBXuKugLj`)M51;abaZ=(7@&U-(mlEm$~F_hCu|KJ%oU&|^9Yb3DT<20(eoxO zOaWoCdFpCQ={?iO`)1J0(kPm(eSwz3#5|fBa9?dW-kYR!Yz0cnxoyT80};9^kW9?; zNJYUJa@B*%Gg#zM<2U-#2dC^VOKh#EpSL_`@Mjc+d&Zfyne3CyzeCe-q>J@DxEB?{ z!*&W-qe8n0(`3?~T*z{F-3Mc5f*$6V;-HN)(r4WX`a2K5f{rc&%>!Fi}xwZ78c>BxTToAfS!h@kNhN9F&5zH|xXT$hSPCBV!JJTbZ z`Fs9u*~i=!7*Jvvk*7rAmluPAzLZvMtgk9D z0rb-1sJQb7dCQWyr`=z6)I`_>KT$N`up%aU2mzu3$g|h=75{y+aPG^o5h%qB(6%g9 zKU*|@vQJjJ;vkBY;yCX#(t^5P?H81lag0US`2Oc&uZhMvQCSCJNl1)_+(SBD$6frQ zk*;+bsn`tQ(UTD5PQF3pO_e#LzdSpj7balJ3aRFSNOi@t2(ujt4l}}h33dHg7wQCY zchF-A^%ChPQVmynXZ}{2S zqYb;6!h}L(!W=dBdybgf=AnTMQ}2)oOlpdwU83F4?E4M4n?;L%eE1t4N7l_xY0xPP z+Ag0BCL7F7->}$&Qd}1PiDgZo;`)X@LxYejrm0^ysKLfXv=`8oc`%A z*LE;d!?S1C6rf{}czo+)q$x@xUCfhRLlg!#Mehl=hla@BAwOmPK07;m7^1?@0erDz zB@u#Xr2duAhG{Z+CTDquotAf@7)R=4~jXLAt%y z$g@SB0C=%YqtBcT14cqlzDD)NVn4>}#5mVhca`ZwLT;&_M)NaO? zwFz?)S$x~J{D4Ag*H^j;r@I~DhYMufB%=SAUBTP(h?3@(SAcy&i&w=KHcMlXCP1m+ zlqpTSOo4q*Lhp^dEk0 zBu*nyRM4N?SDiDLJK5p}(QNORq4q<)C1H=jeM8|!qZtx{j1~t@C{zw&809{^?S4dJo z0$%*Zi(_OPZ3&NeF8FmLP1QVmPbGv`3d({{cj>yMUx^e^*L519&6=_DQl-&^4`Eb^Ymikh z47;lLU@`{p^_pX8eyhd?ByHdNa3K|t$EN_zHbKPC-z{1t$_VWl~b8vn=4sFL;+ zImixVm|_{m-B}AD=q(P$L@Gn6p8ZG7m>nO(mi$MH?K`D90Gnf6zfB#vH6zjX_L*y@ z#ldIxa`{C@a>Kc)c$=`u*VBgRWZ$8OQn)R*NNQ!{2-lZj(sF=-t*!0b85Cr ziHzh*T$p`l3(#=hyy4;2DdNwf8aJ0A{yjUugQN^V!a^u^63Vr>u5Z-5W=WLwIc)v) zOfo`UGpK!Y`pJAIA*KYCsG}ijKgHXn-5+g~uASegO;jgHF2h@sm(_u^g#(z5DhhB~ zL8cjnmy|^&PRRq+!=kUK7Enms&~E-RJE_PqCk&hi@=)2Wzc}84v7O4QzEpNn;mQ%h zQqlERMV?SfqTDuOx&Hcy3u;gIA3Z75XX}3j<{N1}nE&0{=eUUNcH$F(sFUByA0CAq!v+~I-&c_FVXIv*^58g%1K{hnH}qRQi(r1T9}HbOD(POrZ{E{B z3}B4`WQ%+T;s04rqq8X;6oBw(qn&fMAnqHZ0Z3-{J;eo%_k3f{oI8h%3#(mSxYaxM{o~_Jy-C-Da^bETjPx$S`>$f^XjyDzb}%H7u_X#XKXmj;H*) z>#Uwn=y?P0#>;}RXM6W6+dVjoUWVPoftOupG`uu*?K;%2W!;tKYrp38Dx3j-4N#H! z3G>`<$9uVmF}4&-bvhigJGK{ddEdB(4nhPva7>V9-6_v{L_zH{yh?L$E0M_!!U7M3 zZ;Dx^Y>2+)B&B3*KU$0AmxWiprNG?WtIC72fOL~O~RwM?s^hNh)B*KVU z!OF4%SW+EZwChO-aVLO8>!=qPEc@xr`H}tQMXH*E%|8nLS(Oec6&ZAA8yGqP{jrDx zAklYj!(e5V)=Q|i(T1kWBtmOk*h z!|NLD>qsRsqSKkqOZWNoHM|~S37I_BvNZDU&V43PuUX_23a!_) zG7cxi5i1mDEFRd4V=ya?Pw$HAuk9qgFj*cg8(h{oTks9{T-?5KsGe~+B@_cK;Y9U4AmqM&xxVKmkS6UoJ_`WQB<`V(^6dhFmLVXO}wPY@)RD1r{*!~R^! z9-_Uyu~!%n(V0Mes6fn3*suB`v|$TIUx%_%7hi=+8gKMD6auq^8`0$iH;_GMgWXxb z=lCblu=tbEXN=*D{8zkH&5QbL3SeX)il;oBqk?EWn5siD`FN>>+-3`(aD{m57Xhpl z?jYuov*Z~a>|4hR{)<~-YK$?N>Q0FL@hP~?A-d4*lLr`lfZxM^jSWftnM-R1a3=mD zRyhRbc#ZTjv_6_5l320b;ng34PPBN?_zqtDpLVEPga?iP!;Ak@^tPGB|M245_ko*j z+3Qv@;0M{GLaww+cdqC|ePF~?(f;uuARyS>{3T@E3jZNN1j`-*XXoXn;us4>_zi)) Mtg=jo^zDcL2NYTRR{#J2 literal 0 HcmV?d00001 diff --git a/docs/how_to/validation.md b/docs/how_to/validation.md index 4a9803d9..60de83e1 100644 --- a/docs/how_to/validation.md +++ b/docs/how_to/validation.md @@ -231,7 +231,6 @@ VirtualiZarr provides some utilities for creating virtual datasets with such spa The `ManifestArray.with_fill_value_only` method returns a new ManifestArray with the same schema (shape, chunks, codecs, dimension names, attributes) as a given ManifestArray, but with an empty chunk manifest and the given `fill_value`. This is useful for filling in missing files/variables, by creating virtual datasets containing manifestarrays which have no chunk references. -These empty manifestarrays can then be concatenated with the manifestarrays containing chunk references, to create a virtual dataset which a logical grid spanning regions with both present and missing data. ```python exec="on" session="usage" source="material-block" result="code" # a virtual variable from a file where the data is present @@ -255,6 +254,16 @@ print(combined) The combined variable is still a single `ManifestArray` of shape `(2, 3, 4)`, but only the first timestep references real chunks. Reading data from the second timestep would return the array's `fill_value`. +In the previously mentioned case when a dataset is "inherently sparse at a global level, whilst still being dense at a regional level", when users want to concatenate these sparse datasets they can use normal xarray concatenation patterns with outer join. + +``` +xr.concat([ds_1, ds_2], dim="time", join="outer") +``` +This will internally align the sparse datasets on a larger logical grid spanning both regions with null chunks for the sparse regions where no source chunks exist prior to concatenating them. + +![Aligning two regionally-dense but globally-sparse datasets onto a shared logical grid before concatenation along time](grid_alignment.png) + + ### Inhomogeneous Codecs Zarr arrays are by definition made up of a grid of chunks with identical codecs (e.g. for compression). diff --git a/examples/V2/README.md b/examples/V2/README.md index 269a6ad3..da461c49 100644 --- a/examples/V2/README.md +++ b/examples/V2/README.md @@ -20,6 +20,10 @@ Virtualizes the same GOES-16 file using `CachingReadableStore` and `SplittingRea uv run --script examples/V2/goes_with_caching_stores.py ``` +### [its_live.ipynb](its_live.ipynb) + +Mosaics two [ITS_LIVE](https://its-live.jpl.nasa.gov/) glacier-velocity granules that share a global grid but each cover only part of it. Uses native xarray `concat(..., join="outer")` to align them onto the shared grid (sparse, with fill where a granule has no data) and stack them along `time`, keeping the data variables virtual (`ManifestArray`) throughout, then writes the result to Icechunk without copying any pixel data. + ## Performance Comparison Run both scripts and compare the "Virtualization time" printed at the end of each: diff --git a/examples/V2/its_live.ipynb b/examples/V2/its_live.ipynb new file mode 100644 index 00000000..e53068a3 --- /dev/null +++ b/examples/V2/its_live.ipynb @@ -0,0 +1,288 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "its-live-intro", + "metadata": {}, + "source": [ + "# Mosaicking ITS\\_LIVE granules into a virtual cube\n", + "\n", + "This notebook virtualizes two [ITS\\_LIVE](https://its-live.jpl.nasa.gov/) glacier-velocity\n", + "granules that share a global grid but each cover only part of it, then uses native xarray\n", + "`concat(..., join=\"outer\")` to align them onto the shared grid (sparse, with fill where a\n", + "granule has no data) and stack them along `time`. The data variables stay virtual\n", + "(`ManifestArray`) throughout — only the coordinates are loaded — and the result is written to\n", + "Icechunk without copying any pixel data.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "adc1e1cd-7d0f-4345-a6b0-54a70ab044de", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "import numpy as np\n", + "import obstore\n", + "import xarray as xr\n", + "from obspec_utils.registry import ObjectStoreRegistry\n", + "\n", + "import virtualizarr as vz\n", + "from virtualizarr.parsers import HDFParser" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5cf7b30-577a-46c3-ac48-a6f62ba1ea86", + "metadata": {}, + "outputs": [], + "source": [ + "bucket = \"s3://its-live-data\"\n", + "key = \"test-space/virtual-cubes\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "356d1624-041a-4928-a45c-ae0070ce0b56", + "metadata": {}, + "outputs": [], + "source": [ + "granule_1 = \"LC08_L1GT_020121_20231013_20231102_02_T2_X_LC09_L1GT_020121_20231106_20231106_02_T2_G0120V02_P084.nc\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42bca045-986d-4d17-834b-7526ac8afabb", + "metadata": {}, + "outputs": [], + "source": [ + "store = obstore.store.from_url(bucket, region=\"us-west-2\", skip_signature=True)\n", + "registry = ObjectStoreRegistry({bucket: store})\n", + "parser = HDFParser(drop_variables=[\"mapping\", \"img_pair_info\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce09e8ef-8959-459f-8c1f-4676cb21a735", + "metadata": {}, + "outputs": [], + "source": [ + "vds_1 = vz.open_virtual_dataset(\n", + " url=os.path.join(bucket, key, granule_1),\n", + " parser=parser,\n", + " registry=registry,\n", + " # load x/y/time as real, indexed coordinates so xarray has an index on every\n", + " # dimension for alignment. Data variables stay virtual (ManifestArray).\n", + " loadable_variables=[\"time\", \"y\", \"x\"],\n", + " decode_times=True,\n", + ")\n", + "vds_1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "884dffd7-f125-4c1f-add0-25075fc7ce14", + "metadata": {}, + "outputs": [], + "source": [ + "granule_2 = \"LC08_L1GT_020120_20201121_20210315_02_T2_X_LC08_L1GT_020120_20210124_20210305_02_T2_G0120V02_P051.nc\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2198e8ee-62fc-431d-986b-b4f59a9a4667", + "metadata": {}, + "outputs": [], + "source": [ + "vds_2 = vz.open_virtual_dataset(\n", + " url=os.path.join(bucket, key, granule_2),\n", + " parser=parser,\n", + " registry=registry,\n", + " loadable_variables=[\"time\", \"y\", \"x\"],\n", + " decode_times=True,\n", + ")\n", + "vds_2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e7d4099-bbfc-416b-b0e4-d8eb1bfb0b57", + "metadata": {}, + "outputs": [], + "source": [ + "# The ITS_LIVE data uses a descending y coordinate. Internally, xarray's alignment machinery requires ascending\n", + "# ordering so the y coordinate needs to be flipped prior to concatenation.\n", + "flip = lambda ds: ds.assign_coords(y=-ds.y) # descending y -> ascending -y\n", + "flip_1 = flip(vds_1)\n", + "flip_2 = flip(vds_2)\n", + "cube = xr.concat(\n", + " [flip_1, flip_2], dim=\"time\", join=\"outer\", combine_attrs=\"drop_conflicts\"\n", + ")\n", + "cube = flip(cube)\n", + "cube" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d2b78bc-9560-4a0a-b55c-21933e11ff21", + "metadata": {}, + "outputs": [], + "source": [ + "import shutil\n", + "\n", + "import icechunk as ic\n", + "\n", + "# the referenced chunk data lives on this public, anonymous S3 bucket\n", + "url_prefix = \"s3://its-live-data/\"\n", + "store_path = \"its_live_cube.icechunk\"\n", + "\n", + "# start fresh so this cell is re-runnable\n", + "shutil.rmtree(store_path, ignore_errors=True)\n", + "\n", + "# register a virtual chunk container so icechunk knows how to read the s3 refs\n", + "config = ic.RepositoryConfig.default()\n", + "config.set_virtual_chunk_container(\n", + " ic.VirtualChunkContainer(\n", + " url_prefix, ic.s3_store(region=\"us-west-2\", anonymous=True)\n", + " )\n", + ")\n", + "repo = ic.Repository.create(\n", + " storage=ic.local_filesystem_storage(store_path),\n", + " config=config,\n", + " authorize_virtual_chunk_access=ic.containers_credentials(\n", + " {url_prefix: ic.s3_credentials(anonymous=True)}\n", + " ),\n", + ")\n", + "\n", + "\n", + "# icechunk metadata is strict JSON and cannot represent NaN/inf, but ITS_LIVE\n", + "# stores NaN in some attributes (e.g. 'stable_shift_stationary'); drop those.\n", + "def _drop_nonfinite_attrs(ds):\n", + " ds = ds.copy()\n", + "\n", + " def clean(attrs):\n", + " keep = {}\n", + " for k, v in attrs.items():\n", + " arr = np.asarray(v)\n", + " if arr.dtype.kind == \"f\" and not np.isfinite(arr).all():\n", + " continue\n", + " keep[k] = v\n", + " return keep\n", + "\n", + " ds.attrs = clean(ds.attrs)\n", + " for var in ds.variables.values():\n", + " var.attrs = clean(var.attrs)\n", + " return ds\n", + "\n", + "\n", + "# write the cube: only the virtual references are stored, no pixel data is copied\n", + "session = repo.writable_session(\"main\")\n", + "_drop_nonfinite_attrs(cube).vz.to_icechunk(session.store)\n", + "snapshot_id = session.commit(\n", + " \"its_live virtual cube: 2 granules mosaicked + stacked on time\"\n", + ")\n", + "print(\"committed snapshot\", snapshot_id)\n", + "\n", + "# reopen from the committed store\n", + "cube_roundtrip = xr.open_zarr(\n", + " repo.readonly_session(\"main\").store, consolidated=False, zarr_format=3\n", + ")\n", + "cube_roundtrip" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6aa95554-5a0c-481c-bc2a-945576737e00", + "metadata": {}, + "outputs": [], + "source": [ + "v = cube_roundtrip[\"v\"]\n", + "dv = v.isel(time=1) - v.isel(time=0)\n", + "dv" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd10354b", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "# surface speed `v` at each time step, on a shared color scale\n", + "v = cube_roundtrip[\"v\"].load()\n", + "vmax = float(np.nanpercentile(v, 98))\n", + "\n", + "nt = v.sizes[\"time\"]\n", + "fig, axes = plt.subplots(1, nt, figsize=(7 * nt, 6), constrained_layout=True)\n", + "axes = np.atleast_1d(axes)\n", + "for ax, t in zip(axes, range(nt)):\n", + " im = v.isel(time=t).plot.imshow(\n", + " ax=ax, cmap=\"viridis\", vmin=0, vmax=vmax, add_colorbar=False\n", + " )\n", + " ax.set_aspect(\"equal\")\n", + " ax.set_title(str(v[\"time\"].values[t])[:10])\n", + "fig.colorbar(im, ax=list(axes), label=\"speed (m/yr)\", shrink=0.8)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1cc2f9b4-ce0b-4614-8b91-e2416390434e", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "lim = np.nanpercentile(np.abs(dv), 98)\n", + "\n", + "# plot\n", + "fig, ax = plt.subplots(figsize=(9, 7))\n", + "dv.plot.imshow(\n", + " ax=ax,\n", + " cmap=\"RdBu_r\", # red = speedup, blue = slowdown\n", + " vmin=-lim,\n", + " vmax=lim, # symmetric, centered on zero\n", + " cbar_kwargs={\"label\": \"Δ speed (m/yr)\"},\n", + ")\n", + "ax.set_aspect(\"equal\")\n", + "plt.tight_layout()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}