[NEP-001 Feature] byte-identical port np.save/np.load and variants in NumPy 2.4.2

.claude/CLAUDE.md

@@ -297,7 +297,77 @@ Semantics (all probed against NumPy 2.4.2, pinned in `NDEvaluateTests.cs`):
 `bernoulli`, `beta`, `binomial`, `chisquare`, `choice`, `dirichlet`, `exponential`, `f`, `gamma`, `geometric`, `gumbel`, `hypergeometric`, `laplace`, `logistic`, `lognormal`, `logseries`, `multinomial`, `multivariate_normal`, `negative_binomial`, `noncentral_chisquare`, `noncentral_f`, `normal`, `pareto`, `permutation`, `poisson`, `power`, `rand`, `randint`, `randn`, `random_sample`, `rayleigh`, `seed`, `shuffle`, `standard_cauchy`, `standard_exponential`, `standard_gamma`, `standard_normal`, `standard_t`, `triangular`, `uniform`, `vonmises`, `wald`, `weibull`, `zipf`
 
 ### File I/O
-`fromfile`, `load`, `save`, `tofile`
+`fromfile`, `load`, `load_npy`, `load_npz`, `save`, `savez`, `savez_compressed`, `tofile`
+
+The `.npy`/`.npz` stack is a **port of NumPy 2.4.2's `numpy/lib/_format_impl.py`** (NEP-01) in
+`IO/{NpyFormat,NpzFile,PyLiteral}.cs`, and the writer is **byte-for-byte identical to NumPy's own
+`np.save`** — not merely readable by it. Format versions **1.0 / 2.0 / 3.0** (2- vs 4-byte header
+length; latin-1 vs UTF-8), **64-byte `ARRAY_ALIGN`** data alignment (so the data is mmap-ready),
+**`fortran_order`** both ways (write via the transpose; read via reshape-reversed + transpose), and
+**big-endian** files byte-swapped to native on read (per-component: `>c16` swaps in 8s, `>U1` in 4s).
+Headers parse through a real Python-literal parser (`PyLiteral`, standing in for `ast.literal_eval`)
+rather than `IndexOf`/regex, so key order, whitespace, quoting, trailing commas and the Python 2 `3L`
+suffix all work; `max_header_size` (default 10000) bounds hostile input and `allow_pickle` lifts it.
+
+`np.load` returns `object` — `NDArray` for `.npy`, `NpzFile` for `.npz` — mirroring NumPy's
+content-dependent return; **`np.load_npy` / `np.load_npz` are the typed forms** and need no cast.
+`NpzFile` is lazy + cached, takes `"w"` or `"w.npy"` as keys, exposes `.Files` stripped, supports
+`npz.f.weights` dot access (NumPy's `BagObj`) and **must be disposed**.
+
+| Dtype | Maps to | Notes |
+|-------|---------|-------|
+| Boolean/SByte/Byte | `\|b1` / `\|i1` / `\|u1` | single-byte types take the `\|` prefix, never `<` |
+| Int16…UInt64, Half, Single, Double | `<i2`…`<u8`, `<f2`, `<f4`, `<f8` | direct |
+| Complex | `<c16` | `<c8` widens to `Complex` on read |
+| Char | `<U1` | 2-byte UTF-16 ↔ 4-byte UCS-4; non-BMP rejected (needs a surrogate pair) |
+| Decimal | — | `NotSupportedException`: no NumPy dtype (cast to Double) |
+
+Object arrays, structured/subarray dtypes, `datetime64`/`timedelta64`, `\|S`/`<U`n>1/`\|V`, `<f16`
+and `<c32` are parsed then rejected with a precise message — including the zero-width `<U0`/`\|S0`/`\|V0`,
+which are *valid* NumPy dtypes and so report "unsupported", never "invalid descriptor". `mmap_mode`
+validates its value and throws `NotImplementedException` (sketch in `np.load.cs`'s `CheckMmapMode`).
+
+**Known intentional divergences** (differential-verified, everything else agrees): the 5 unsupported
+dtypes above are the *only* files NumPy loads that NumSharp refuses. NumPy reports `shape: (True,)`
+and overflowing dims via a later TypeError/OverflowError while NumSharp rejects them up front — both
+refuse, the text differs. Big-endian files byte-swap to native (NumPy keeps a byte-swapped dtype),
+and `Char` round-trips U+0000 as a NUL where NumPy's `<U` reports `''` (the bytes are identical —
+`<U` treats NUL as string padding).
+
+**Three traps this port has already fallen into — do not re-break:**
+- **`(1)` is not a tuple.** Only a comma makes one in Python, so `'shape': (1)` is the *int* 1 and
+  NumPy rejects it. `PyLiteral.ParseTuple` tracks `sawComma` for exactly this.
+- **A header-length field is attacker-controlled.** `ReadBytes` sizes its buffer from what the stream
+  actually holds (seekable ⇒ exact) and grows as bytes arrive — never from `count`. A 28-byte file
+  can claim a 4 GB header, and `new byte[claim]` turns that into a gigabyte spike; NumPy shrugs it
+  off in ~2 KB because `fp.read(n)` allocates only what it returns. Pinned by
+  `HostileHeaderLength_DoesNotAllocateTheClaim` (measures ~1.2 KB).
+- **`NpzFile`'s name map needs NumPy's TWO passes**, not one interleaved loop: `dict(zip(stripped,
+  full))` then `.update(zip(full, full))`. For an archive holding both `a` and `a.npy` — which both
+  strip to key `a` — pass 2 re-points `a` at the entry literally *named* `a`; one loop resolves it to
+  whichever came last. Within a pass, later entries win (a zip may repeat a name, and Python's dict
+  keeps the last; .NET's `GetEntry` would return the first). Pinned by the `npz_*_names` cases.
+- **Never buffer an npz member to sniff its magic.** `MemoryStream`/`byte[]` cap at 2 GB, so buffering
+  made NumSharp *write* 5 GiB archives it then failed to *read* ("Stream was too long") while NumPy
+  read them fine — and it doubled the cost of every ordinary load. `NpyFormat.ReadArray` is strictly
+  forward-only, so `StartsWithNpyMagic` sniffs on one `entry.Open()` and the reader streams from a
+  fresh one. Pinned by `Npz_StreamsMembers_RatherThanBufferingThem`.
+
+**Size limits (measured, not assumed):** `.npy` and `.npz` both round-trip **> 4 GiB** — verified at
+5 GiB against every 32-bit boundary (`int.MaxValue`, `uint.MaxValue`) with NumPy reading the result,
+and the 5 GiB header is byte-identical to NumPy's. Array data lives in unmanaged memory and both IO
+paths stream in 256 KB chunks, so no managed buffer ever scales with the array. Zip64 comes from
+.NET's `ZipArchive` on demand (NumPy forces `force_zip64=True` for the same reason, gh-10776). The
+only 2 GB ceilings left are the inherently `byte[]`-returning conveniences — `np.save(arr)`,
+`np.savez(...)` → `byte[]`, and `NpzFile.GetRawBytes` — which is a `byte[]` limit, not a format one;
+use the path/stream overloads above 2 GB. Tests: `NpyLargeFileTests` (`[HighMemory][LongIndexing]`,
+excluded from CI; needs ~8 GB RAM + ~6 GB disk).
+
+**Gates:** `test/oracle/gen_npy_oracle.py` → `IO/corpus/npy_oracle.zip` (286 committed cases of real
+NumPy output) replayed by `NpyOracleTests` under the **`NpyOracle`** category — read / byte-exact
+write / header-only write / verbatim error / round-trip / npz / live-view write / hostile-allocation.
+Reverse interop (NumPy reading NumSharp, the direction byte-equality can't cover for `.npz`) is the
+manual gate `python test/oracle/verify_npy_interop.py`.
 
 ### Printing / Formatting
 `array2string`, `array_repr`, `array_str`, `format_float_positional`, `format_float_scientific`, `get_printoptions`, `printoptions` (IDisposable context), `set_printoptions`
@@ -342,6 +412,7 @@ Semantics (all probed against NumPy 2.4.2, pinned in `NDEvaluateTests.cs`):
 | ILKernelGenerator | `Backends/Kernels/ILKernelGenerator*.cs` (per-chunk, NDIter-driven) |
 | DirectILKernelGenerator | `Backends/Kernels/Direct/DirectILKernelGenerator.*.cs` (whole-array, 63 partials) |
 | Kernel shared infra | `Backends/Kernels/{VectorMethodCache,ScalarMethodCache,KernelOp,StrideDetector,SimdMatMul.*}.cs` |
+| .npy/.npz format (NEP-01) | `IO/NpyFormat.cs` (port of NumPy's `_format_impl.py`), `IO/NpzFile.cs` (lazy archive), `IO/PyLiteral.cs` (header parser ≙ `ast.literal_eval`), `APIs/np.{save,load}.cs` |
 | Type info | `Utilities/InfoOf.cs` |
 | Generic NDArray | `Generics/NDArray\`1.cs` |
 
@@ -561,6 +632,10 @@ test/oracle/                          corpus generators (NumPy 2.4.2 — run by
                                       varstd,matmul,astype,stat,where,sort,manip}.jsonl
   gen_index_oracle.py                 getter/setter index oracle (token index over 15 base recipes)
   fuzz_random.py                      seeded random fuzzer (imports the other two)
+  gen_npy_oracle.py                   .npy/.npz format oracle: REAL np.save/savez output + a manifest ->
+                                      NumSharp.UnitTest/IO/corpus/npy_oracle.zip (286 cases)
+  interop_write.cs                    the reverse direction: NumSharp writes, verify_npy_interop.py reads
+  verify_npy_interop.py               manual gate — real NumPy reads NumSharp's output (needs Python + SDK)
 test/NumSharp.UnitTest/Fuzz/          C# replay harness (no Python)
   FuzzCorpus.cs                       rebuilds EXACT NDArray views from (dtype,shape,strides,offset,bytes)
   OpRegistry.cs                       op-name → NumSharp call (pairs 1:1 with gen_oracle.py)
@@ -571,6 +646,20 @@ test/NumSharp.UnitTest/Fuzz/          C# replay harness (no Python)
   MetamorphicTests.cs                 oracle-free invariants (round-trips / involutions / identities — no NumPy)
   HarnessSelfTests.cs                 proves the harness has teeth (BitDiff detects value/NaN/-0 diffs; non-vacuous)
   corpus/*.jsonl                      committed corpus (~68K cases / 40 tiers; op corpus ~53K incl. 3.7K Char woven + 579 Decimal), copied to test output via the csproj glob
+test/NumSharp.UnitTest/IO/            .npy/.npz format gate (no Python)
+  NpyOracleCorpus.cs                  opens npy_oracle.zip, rebuilds arrays from the manifest
+  NpyOracleTests.cs                   one [NpyOracle] test per claim: read / byte-exact write / header-only
+                                      write / verbatim error / round-trip / npz / live-view write /
+                                      hostile-allocation / corpus non-vacuity
+  NpyLargeFileTests.cs                the 2 GB and 4 GB walls; the 5 GiB pair is [HighMemory] (not in CI)
+  AllocationTests.cs                  MinAllocated: best-of-N, since GetTotalAllocatedBytes is process-wide
+  corpus/npy_oracle.zip               committed real-NumPy output (286 cases), copied via the csproj glob
+  NumpyCompatibilityTests.cs          hand-made NumPy fixtures (test_compat/) — v2.0/v3.0/fortran/scalar/
+  NpyFormatVersionTests.cs            empty/bool/multi/compressed. They resolve test_compat RELATIVE to the
+                                      working dir and must fail loudly when it is absent: they used to
+                                      `return` silently, and nothing copied the fixtures, so on a clean
+                                      checkout all 13 passed while asserting NOTHING. The csproj now copies
+                                      test_compat\* — do not drop that rule.
 ```
 
 - **Generators live in `test/oracle/`** and write the corpus into `test/NumSharp.UnitTest/Fuzz/corpus/` (path resolved relative to `test/oracle/`). CI replays the committed corpus, never the generators.
@@ -579,6 +668,27 @@ test/NumSharp.UnitTest/Fuzz/          C# replay harness (no Python)
 - **Regenerate** (deterministic; needs `numpy==2.4.2`): `python test/oracle/gen_oracle.py <mode>` (modes: `smoke astype_full binary divmod_power comparison unary reduce where place matmul bitwise unary_extra nanreduce scan stat logic modf manip sort tail params aliasing copyto errors`) + `python test/oracle/gen_index_oracle.py` (the `index_*` tiers) + `python test/oracle/fuzz_random.py 1234 2000 random_smoke.jsonl` + `dotnet run test/oracle/gen_decimal_oracle.cs` (the `decimal_*` tiers), then `dotnet build` (copies the corpus to test output).
 - **Run the gate**: `dotnet test --filter "TestCategory=FuzzMatrix"`. Each case is bit-exact (pass), a documented difference in `MisalignedRegistry` (excused, never silent), or a failure (red). Full divergence ledger: `test/NumSharp.UnitTest/Fuzz/README.md`.
 
+### The `.npy`/`.npz` format oracle (same philosophy, separate corpus)
+
+Same shape as above — NumPy is the oracle, the corpus is committed, no Python at test time — but the
+claim is stronger: NumSharp's writer must be **byte-identical** to `np.save`, not merely readable.
+
+- **Regenerate**: `python test/oracle/gen_npy_oracle.py` (deterministic; needs `numpy==2.4.2`), then `dotnet build`.
+- **Run the gate**: `dotnet test --filter "TestCategory=NpyOracle"` — 286 cases across every dtype ×
+  {0-d, empty, 1-D, 2-D, 3-D, unit, empty-2d/3d} × {C, F, strided, reversed, offset, broadcast,
+  transposed} × versions {1.0, 2.0, 3.0} × {little, big, native} endian, plus value fidelity (NaN
+  payloads incl. sNaN, subnormals, signed zero, integer extremes, BMP seams) and 42
+  malformed / unsupported / hostile files whose messages must match NumPy verbatim.
+- **`kind: "header"` cases exist because two writer branches are unreachable from a real array.** The
+  growth padding is normally INVISIBLE — shrink the body by 5 chars and the alignment padding grows by
+  5, leaving the file identical — so a wrong growth axis (`shape[0]` vs `shape[-1]`) passes every
+  ordinary test; it only shows on shapes that tip the header across a 64-byte bucket, and those need
+  10¹⁷ elements to allocate. The v1.0→v2.0 auto-selection boundary likewise sits at 21,817 dimensions.
+  Both are driven through the header dict, exactly as NumPy's `_write_array_header` was.
+- **Reverse interop** (manual; needs Python + the SDK): `python test/oracle/verify_npy_interop.py` has
+  NumSharp write 28 files and real NumPy read them. This is the only way to prove `.npz` output,
+  whose ZIP framing legitimately differs from Python's `zipfile` and so cannot be byte-compared.
+
 ## CI Pipeline
 
 `.github/workflows/build-and-release.yml` — test on 3 OSes (Windows/Ubuntu/macOS), build NuGet on tag push, create GitHub Release, publish to nuget.org. The `FuzzMatrix` gate runs here (replays the committed corpora; no Python).
@@ -759,7 +869,10 @@ A: No external runtime dependencies. `System.Memory` and `System.Runtime.Compile
 A: TensorFlow.NET, ML.NET integrations, Gym.NET, Pandas.NET, and various scientific computing projects.
 
 **Q: Can I save/load NumPy files?**
-A: Yes. `np.save()` writes `.npy` files, `np.load()` reads both `.npy` and `.npz` archives. Compatible with Python NumPy files.
+A: Yes, and byte-exactly. `np.save` writes `.npy`, `np.savez`/`np.savez_compressed` write `.npz`, and `np.load` reads both (returning `object`; `np.load_npy`/`np.load_npz` are the typed forms). `IO/NpyFormat.cs` is a port of NumPy 2.4.2's `_format_impl.py`, so a saved file is byte-for-byte what NumPy's own `np.save` would produce — versions 1.0/2.0/3.0, 64-byte alignment, `fortran_order`, and big-endian reads included. See the File I/O section above for the dtype map and the gates.
+
+**Q: Why does `np.load` return `object` instead of `NDArray`?**
+A: Because NumPy's does: `np.load` yields an `ndarray` for a `.npy` and an `NpzFile` for a `.npz`, dispatching on the file's magic bytes rather than its extension, and C# has no union type to express that. `np.load_npy` and `np.load_npz` are the typed escapes — prefer them when the kind is known; each also reports a directed error if handed the other kind.
 
 **Q: What random distributions are supported?**
 A: An extensive NumPy-matching set (40+ distributions and samplers — see the Supported np.* APIs → Random list), all on the `NumPyRandom` class in `RandomSampling/`, with 1-to-1 seed/state parity to NumPy.