diff --git a/README.md b/README.md index d3cf58d..d8467b0 100644 --- a/README.md +++ b/README.md @@ -132,6 +132,13 @@ is `Outer`, which has no error result: it panics when its result is too big, which takes operands of more than 2^30 elements each. An empty array whose other axis is huge, `(2^58, 0)`, costs nothing to reduce, concatenate or scan. +A shape that is representable can still be far larger than memory, and Go +cannot recover from running out of it: the process dies, where NumPy raises +`MemoryError`. Several operations produce more than they are given: +broadcasting (`(n, 1) + (1, n)` is n² elements), `Outer`, `MatMul`, and summing +the empty axis of `(0, n)` (n zeros, as in NumPy). When shapes come from +someone else, bound the result's size before computing it. + The assembly kernels take raw pointers. Their Go wrappers check every operand the kernel reads or writes against the length it is given, so a mistake is an index panic in Go, never an access outside a slice; a test places each operand @@ -156,6 +163,7 @@ cols, _ := m.SumAxis(0, false) // sum down each column -> [3 5 7] // NumPy basic-indexing views (share data): col0, _ := m.Slice(ndarray.All(), ndarray.A(0)) // m[:,0] -> [0 3] sub, _ := m.Slice(ndarray.R(0, 2), ndarray.Step(2)) // m[0:2, ::2] +row1, _ := m.Slice(ndarray.A(1)) // m[1]: unindexed axes are whole // ufuncs and masks roots := m.Sqrt() diff --git a/SECURITY.md b/SECURITY.md index 7cad1a7..ff81a2a 100644 --- a/SECURITY.md +++ b/SECURITY.md @@ -39,11 +39,40 @@ bounds checks do not cover, so they are held to the following: `TestEmptyArraysWithHugeAxesReturnPromptly` covers concatenation, the axis reductions and the scans. +- **Memory is the caller's to bound.** Validation guarantees a shape is + representable, not that it fits in memory, and a Go program that runs out of + memory dies (NumPy raises `MemoryError`). Broadcasting, `Outer`, `MatMul` and + summing the empty axis of `(0, n)` produce more than their operands hold; + with untrusted shapes, check the result's size first. +- **Fuzzed against an oracle.** `FuzzOps` (`fuzz_test.go`) drives slices, + reshapes, transposes, squeezes, broadcasting arithmetic, axis reductions, + `Take`, `MaskSelect`, concatenation, `MatMul` and scans from arbitrary bytes, + and compares every result element by element with an oracle written in the + test from NumPy's definitions (Python's `slice.indices`, the broadcasting + rule, ...). The library must reject exactly the calls the oracle rejects and + never panic. + The minimum Go version (`go` in `go.mod`) is kept on a release that fixes the standard library's known vulnerabilities; `govulncheck ./...` reports none. ## Audit log +**2026-10-04 (v0.4.0).** + +- `FuzzOps` ran 47 million executions without a divergence. A first version + that only checked an array against its own `Copy` was **blind**: with + `Squeeze` deliberately given wrong strides it found nothing in 31 million + executions, because a wrong view reads the same wrong values both ways. The + oracle version finds the same sabotage in 11 seconds; its failing input is + kept in `testdata/fuzz/FuzzOps`. +- It found two places where the library and NumPy differ. `Slice` with fewer + indices than axes was an error; it now takes the rest whole, as NumPy does. + `MaskSelect` accepts a mask that broadcasts to the array, which NumPy + rejects; this is kept as a documented extension. +- `SumAxis`/`ProdAxis` along a zero-length axis now return 0s/1s as NumPy does + (they were errors, unlike `Sum`/`Prod` of an empty array). This makes the + memory caveat above concrete, and it is documented there. + **2026-10-04 (v0.3.0).** - `govulncheck` found no reachable vulnerability. At `go 1.26.4`, 10 standard diff --git a/fuzz_test.go b/fuzz_test.go new file mode 100644 index 0000000..30d9815 --- /dev/null +++ b/fuzz_test.go @@ -0,0 +1,610 @@ +package ndarray + +import ( + "math" + "testing" +) + +// FuzzOps is a differential fuzzer over the public API. Arbitrary bytes build +// arrays of small (possibly empty, possibly invalid) shapes and choose a +// sequence of operations and their arguments: slices with any bounds and +// steps, reshapes with -1, transposes, squeezes, broadcasting arithmetic, axis +// reductions, indexing, concatenation and matrix products. +// +// Every result is compared with an oracle written here from NumPy's +// definitions, independently of the library: the oracle reads the operands +// element by element with At and computes what each output element must be. +// A self-consistency check alone is blind: a view with wrong strides reads +// back the same wrong values through At and Copy. Values are small integers, +// so every sum and product is exact and the comparison is equality. +// +// The oracle also says whether the call is valid; the library must agree, +// returning an error exactly when the oracle rejects the arguments, and must +// never panic. The seed corpus runs as an ordinary test; `go test -fuzz +// FuzzOps` explores. +func FuzzOps(f *testing.F) { + for _, s := range [][]byte{ + {}, + {0, 2, 3, 1, 1, 2, 0, 0, 3, 5}, + {1, 0, 3, 7, 0, 255, 2, 1, 9, 9, 9}, + {2, 4, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}, + {3, 1, 1, 3, 0, 2, 254, 1, 4, 0, 6, 2}, + {0, 3, 3, 2, 11, 1, 12, 0, 13, 2, 14}, + } { + f.Add(s) + } + f.Fuzz(func(t *testing.T, data []byte) { + r := &byteReader{b: data} + var pool []*Array + for range 2 { + if a := r.array(t); a != nil { + pool = append(pool, a) + } + } + for step := 0; step < 12 && r.more() && len(pool) > 0; step++ { + a := pool[r.n(len(pool))] + b := pool[r.n(len(pool))] + if res := applyChecked(t, r, a, b); res != nil && len(pool) < 6 { + pool = append(pool, res) + } + } + }) +} + +// applyChecked runs one operation chosen by r and checks it against the +// oracle, returning the result when there is one. +func applyChecked(t *testing.T, r *byteReader, a, b *Array) *Array { + t.Helper() + var ( + got *Array + err error + want *ref // nil: the oracle rejects the call + name string + ) + switch r.n(14) { + case 0: + specs := make([]sliceSpec, r.n(4)) + idx := make([]Index, len(specs)) + for i := range specs { + specs[i] = sliceSpec{kind: r.n(4), x: r.i(), y: r.i(), z: r.i()} + idx[i] = specs[i].index() + } + name = "Slice" + got, err = a.Slice(idx...) + want = refSlice(a, specs) + case 1: + shape := make([]int, r.n(4)) + for i := range shape { + shape[i] = r.i() + } + name = "Reshape" + got, err = a.Reshape(shape...) + want = refReshape(a, shape) + case 2: + name = "Transpose" + got, err = a.Transpose(), nil + want = refTranspose(a) + case 3: + name = "Squeeze" + got, err = a.Squeeze() + want = refSqueeze(a) + case 4: + ax := r.i() + name = "ExpandDims" + got, err = a.ExpandDims(ax) + want = refExpandDims(a, ax) + case 5: + name = "Add" + got, err = a.Add(b) + want = refBroadcast(a, b, func(x, y float64) float64 { return x + y }) + case 6: + name = "Mul" + got, err = a.Mul(b) + want = refBroadcast(a, b, func(x, y float64) float64 { return x * y }) + case 7: + ax, keep := r.i(), r.n(2) == 1 + name = "SumAxis" + got, err = a.SumAxis(ax, keep) + want = refReduce(a, ax, keep, 0, func(acc, x float64) float64 { return acc + x }, true) + case 8: + ax, keep := r.i(), r.n(2) == 1 + name = "MaxAxis" + got, err = a.MaxAxis(ax, keep) + want = refReduce(a, ax, keep, math.Inf(-1), math.Max, false) + case 9: + idx := make([]int, r.n(4)) + for i := range idx { + idx[i] = r.i() + } + name = "Take" + got, err = a.Take(idx...) + want = refTake(a, idx) + case 10: + ax := r.i() + name = "Concatenate" + got, err = Concatenate([]*Array{a, b}, ax) + want = refConcat(a, b, ax) + case 11: + name = "MatMul" + got, err = a.MatMul(b) + want = refMatMul(a, b) + case 12: + name = "MaskSelect" + got, err = a.MaskSelect(b) + want = refMask(a, b) + default: + ax := r.i() + name = "CumSum" + got, err = a.CumSum(ax) + want = refCumSum(a, ax) + } + switch { + case want == nil && err == nil: + t.Fatalf("%s(%v, %v): accepted, oracle rejects; got shape %v", name, a.Shape(), b.Shape(), got.Shape()) + case want != nil && err != nil: + t.Fatalf("%s(%v, %v): %v; oracle gives shape %v", name, a.Shape(), b.Shape(), err, want.shape) + case want == nil: + return nil + } + checkAgainst(t, name, got, want) + return got +} + +// ref is an oracle result: a shape and its row-major values. +type ref struct { + shape []int + vals []float64 +} + +// checkAgainst asserts got has want's shape and, read element by element +// through its own strides, want's values; and that Copy agrees. +func checkAgainst(t *testing.T, name string, got *Array, want *ref) { + t.Helper() + if !equalInts(got.Shape(), want.shape) { + t.Fatalf("%s: shape %v, oracle %v", name, got.Shape(), want.shape) + } + if got.Size() != len(want.vals) { + t.Fatalf("%s: Size %d, oracle %d", name, got.Size(), len(want.vals)) + } + flat := got.Copy().contiguousData() + k := 0 + eachIndex(want.shape, func(idx []int) { + v := got.At(idx...) + if !same(v, want.vals[k]) || !same(flat[k], want.vals[k]) { + t.Fatalf("%s: At%v = %v, Copy = %v, oracle %v (shape %v)", + name, idx, v, flat[k], want.vals[k], want.shape) + } + k++ + }) +} + +func same(x, y float64) bool { return x == y || (math.IsNaN(x) && math.IsNaN(y)) } + +func equalInts(x, y []int) bool { + if len(x) != len(y) { + return false + } + for i := range x { + if x[i] != y[i] { + return false + } + } + return true +} + +// eachIndex calls f with every index of shape in row-major order. +func eachIndex(shape []int, f func([]int)) { + n := 1 + for _, d := range shape { + n *= d + } + idx := make([]int, len(shape)) + for k := 0; k < n; k++ { + f(idx) + for d := len(idx) - 1; d >= 0; d-- { + idx[d]++ + if idx[d] < shape[d] { + break + } + idx[d] = 0 + } + } +} + +// build evaluates at(idx) for every index of shape. +func build(shape []int, at func([]int) float64) *ref { + out := &ref{shape: append([]int{}, shape...)} + eachIndex(shape, func(idx []int) { out.vals = append(out.vals, at(idx)) }) + return out +} + +func values(a *Array) []float64 { + return build(a.Shape(), func(i []int) float64 { return a.At(i...) }).vals +} + +// normAxis maps a possibly negative axis into [0, n), or -1 if out of range. +func normAxis(ax, n int) int { + if ax < 0 { + ax += n + } + if ax < 0 || ax >= n { + return -1 + } + return ax +} + +// sliceSpec is one fuzzed index: kind 0 = A(x), 1 = Rng(x, y, z), +// 2 = Step(z), 3 = All(). +type sliceSpec struct{ kind, x, y, z int } + +func (s sliceSpec) index() Index { + switch s.kind { + case 0: + return A(s.x) + case 1: + return Rng(s.x, s.y, s.z) + case 2: + return Step(s.z) + } + return All() +} + +// positions is Python's slice.indices: the positions selected along an axis +// of length n, or ok=false for a zero step. +func (s sliceSpec) positions(n int) (pos []int, ok bool) { + step := 1 + hasBounds := false + switch s.kind { + case 1: + step, hasBounds = s.z, true + case 2: + step = s.z + } + if step == 0 { + return nil, false + } + var start, stop int + if step > 0 { + start, stop = 0, n + } else { + start, stop = n-1, -1 + } + if hasBounds { + clamp := func(v int) int { + if v < 0 { + v += n + } + lo, hi := 0, n + if step < 0 { + lo, hi = -1, n-1 + } + return min(max(v, lo), hi) + } + start, stop = clamp(s.x), clamp(s.y) + } + for i := start; (step > 0 && i < stop) || (step < 0 && i > stop); i += step { + pos = append(pos, i) + } + return pos, true +} + +func refSlice(a *Array, specs []sliceSpec) *ref { + shape := a.Shape() + if len(specs) > len(shape) { + return nil + } + // Per input axis: either a fixed position (an integer index, axis dropped) + // or the list of positions it ranges over. + fixed := make([]int, len(shape)) + ranges := make([][]int, len(shape)) + var outShape []int + var outAxes []int // input axis of each output axis + for d, n := range shape { + s := sliceSpec{kind: 3} + if d < len(specs) { + s = specs[d] + } + if s.kind == 0 { + i := s.x + if i < 0 { + i += n + } + if i < 0 || i >= n { + return nil + } + fixed[d], ranges[d] = i, nil + continue + } + pos, ok := s.positions(n) + if !ok { + return nil + } + ranges[d] = pos + outShape = append(outShape, len(pos)) + outAxes = append(outAxes, d) + } + if outShape == nil { + outShape = []int{} + } + return build(outShape, func(idx []int) float64 { + in := append([]int{}, fixed...) + for o, d := range outAxes { + in[d] = ranges[d][idx[o]] + } + return a.At(in...) + }) +} + +func refReshape(a *Array, shape []int) *ref { + n := a.Size() + known, infer := 1, -1 + for i, d := range shape { + switch { + case d == -1 && infer < 0: + infer = i + case d < 0: + return nil + default: + known *= d + } + } + out := append([]int{}, shape...) + if infer >= 0 { + if known == 0 || n%known != 0 { + return nil + } + out[infer] = n / known + } else if known != n { + return nil + } + return &ref{shape: out, vals: values(a)} +} + +func refTranspose(a *Array) *ref { + s := a.Shape() + out := make([]int, len(s)) + for i := range s { + out[i] = s[len(s)-1-i] + } + return build(out, func(idx []int) float64 { + in := make([]int, len(idx)) + for i := range idx { + in[i] = idx[len(idx)-1-i] + } + return a.At(in...) + }) +} + +func refSqueeze(a *Array) *ref { + out := []int{} + for _, d := range a.Shape() { + if d != 1 { + out = append(out, d) + } + } + return &ref{shape: out, vals: values(a)} +} + +func refExpandDims(a *Array, ax int) *ref { + s := a.Shape() + if ax = normAxis(ax, len(s)+1); ax < 0 { + return nil + } + out := append(append(append([]int{}, s[:ax]...), 1), s[ax:]...) + return &ref{shape: out, vals: values(a)} +} + +// refBroadcastShape is NumPy's rule: align on the right, each pair equal or 1. +func refBroadcastShape(x, y []int) ([]int, bool) { + n := max(len(x), len(y)) + out := make([]int, n) + for i := 0; i < n; i++ { + dx, dy := 1, 1 + if j := len(x) - n + i; j >= 0 { + dx = x[j] + } + if j := len(y) - n + i; j >= 0 { + dy = y[j] + } + switch { + case dx == dy, dy == 1: + out[i] = dx + case dx == 1: + out[i] = dy + default: + return nil, false + } + } + return out, true +} + +// atBroadcast reads a at an index of a broadcast shape of rank len(idx). +func atBroadcast(a *Array, idx []int) float64 { + s := a.Shape() + in := make([]int, len(s)) + for i := range s { + if s[i] != 1 { + in[i] = idx[len(idx)-len(s)+i] + } + } + return a.At(in...) +} + +func refBroadcast(a, b *Array, op func(x, y float64) float64) *ref { + out, ok := refBroadcastShape(a.Shape(), b.Shape()) + if !ok { + return nil + } + return build(out, func(idx []int) float64 { return op(atBroadcast(a, idx), atBroadcast(b, idx)) }) +} + +func refReduce(a *Array, ax int, keep bool, init float64, op func(acc, x float64) float64, hasIdentity bool) *ref { + s := a.Shape() + if ax = normAxis(ax, len(s)); ax < 0 { + return nil + } + if s[ax] == 0 && !hasIdentity { + return nil + } + var out []int + for i, d := range s { + switch { + case i != ax: + out = append(out, d) + case keep: + out = append(out, 1) + } + } + if out == nil { + out = []int{} + } + return build(out, func(idx []int) float64 { + in := make([]int, 0, len(s)) + in = append(in, idx[:ax]...) + in = append(in, 0) + rest := idx[ax:] + if keep { + rest = idx[ax+1:] + } + in = append(in, rest...) + acc := init + for k := 0; k < s[ax]; k++ { + in[ax] = k + acc = op(acc, a.At(in...)) + } + return acc + }) +} + +func refTake(a *Array, idx []int) *ref { + flat := values(a) + out := &ref{shape: []int{len(idx)}} + for _, i := range idx { + if i < 0 { + i += len(flat) + } + if i < 0 || i >= len(flat) { + return nil + } + out.vals = append(out.vals, flat[i]) + } + return out +} + +func refConcat(a, b *Array, ax int) *ref { + sa, sb := a.Shape(), b.Shape() + if len(sa) != len(sb) || len(sa) == 0 { + return nil + } + if ax = normAxis(ax, len(sa)); ax < 0 { + return nil + } + for i := range sa { + if i != ax && sa[i] != sb[i] { + return nil + } + } + out := append([]int{}, sa...) + out[ax] += sb[ax] + return build(out, func(idx []int) float64 { + if idx[ax] < sa[ax] { + return a.At(idx...) + } + in := append([]int{}, idx...) + in[ax] -= sa[ax] + return b.At(in...) + }) +} + +func refMatMul(a, b *Array) *ref { + sa, sb := a.Shape(), b.Shape() + if len(sa) != 2 || len(sb) != 2 || sa[1] != sb[0] { + return nil + } + return build([]int{sa[0], sb[1]}, func(idx []int) float64 { + s := 0.0 + for k := 0; k < sa[1]; k++ { + s += a.At(idx[0], k) * b.At(k, idx[1]) + } + return s + }) +} + +// refMask follows MaskSelect's contract: the mask must broadcast to a's shape +// without enlarging it (an extension of NumPy, which wants the same shape). +func refMask(a, m *Array) *ref { + if s, ok := refBroadcastShape(a.Shape(), m.Shape()); !ok || !equalInts(s, a.Shape()) { + return nil + } + out := &ref{} + eachIndex(a.Shape(), func(idx []int) { + if atBroadcast(m, idx) != 0 { + out.vals = append(out.vals, a.At(idx...)) + } + }) + out.shape = []int{len(out.vals)} + return out +} + +func refCumSum(a *Array, ax int) *ref { + s := a.Shape() + if ax = normAxis(ax, len(s)); ax < 0 { + return nil + } + return build(s, func(idx []int) float64 { + in := append([]int{}, idx...) + acc := 0.0 + for k := 0; k <= idx[ax]; k++ { + in[ax] = k + acc += a.At(in...) + } + return acc + }) +} + +// byteReader hands out small integers from the fuzz input, and zeros once it +// is exhausted. +type byteReader struct { + b []byte + p int +} + +func (r *byteReader) more() bool { return r.p < len(r.b) } + +func (r *byteReader) next() byte { + if r.p >= len(r.b) { + return 0 + } + v := r.b[r.p] + r.p++ + return v +} + +// n returns a value in [0, m). +func (r *byteReader) n(m int) int { return int(r.next()) % m } + +// i returns a small signed value in [-8, 8), the range that exercises +// negative indices, out-of-range bounds and reversed steps. +func (r *byteReader) i() int { return int(r.next()%16) - 8 } + +// array builds an array of rank 0-3 with dimensions in [-1, 4] and small +// integer values, or nil when the shape is (correctly) rejected. +func (r *byteReader) array(t *testing.T) *Array { + shape := make([]int, r.n(4)) + for i := range shape { + shape[i] = r.n(6) - 1 + } + a, err := New(shape...) + if err != nil { + for _, d := range shape { + if d < 0 { + return nil + } + } + t.Fatalf("New(%v): %v", shape, err) + } + for i := range a.data { + a.data[i] = float64(int(r.next())%7 - 3) + } + return a +} diff --git a/indexing.go b/indexing.go index 9eac0c0..52aa880 100644 --- a/indexing.go +++ b/indexing.go @@ -11,8 +11,13 @@ import "fmt" // Take is integer fancy indexing into the flat array. // MaskSelect returns a 1-D array of the elements of a where mask is truthy -// (non-zero), in row-major order — NumPy's a[mask]. mask must broadcast to a's -// shape (typically it has exactly a's shape, e.g. a mask from a.Greater(...)). +// (non-zero), in row-major order: NumPy's a[mask] for a mask of a's shape, the +// usual case (a mask from a.Greater(...)). mask may also broadcast to a's shape, +// which NumPy rejects (IndexError): a (3,) mask selects the same columns of +// every row of a (2, 3). Two NumPy forms differ: a mask covering only the +// leading axes, a[rowmask], is rejected here unless it broadcasts; and a 0-d +// mask selects all or none of a, flattened, where NumPy returns a with a new +// leading axis of length 1 or 0. func (a *Array) MaskSelect(mask *Array) (*Array, error) { shape, err := broadcastShape(a.shape, mask.shape) if err != nil { diff --git a/ndarray.go b/ndarray.go index f2423d5..142c6be 100644 --- a/ndarray.go +++ b/ndarray.go @@ -853,8 +853,8 @@ func (a *Array) reduceLayout(axis int) (outer, axisLen, inner int) { // reduceAxis is the shared driver for the axis reductions. It validates the // axis, materialises the data, runs the supplied kernel, and wraps the result // in an array of the reduced shape. It rejects reduction along a zero-length -// axis (the empty-reduction case), matching this package's whole-array -// reductions which error on empty input. +// axis, as numpy does for the reductions with no identity (max, min, argmax, +// argmin); Sum and Prod, which have one, go through reduceAxisOr instead. func (a *Array) reduceAxis( axis int, keepdims bool, kernel func(dst, src []float64, outer, axisLen, inner, lo, hi int), @@ -886,13 +886,37 @@ func (a *Array) reduceAxis( // is kept with length 1 (e.g. (2,3) summed over axis 0 -> (1,3)); otherwise it // is removed (-> (3,)). A negative axis counts from the end. func (a *Array) SumAxis(axis int, keepdims bool) (*Array, error) { - return a.reduceAxis(axis, keepdims, kernels.SumAxis) + return a.reduceAxisOr(axis, keepdims, kernels.SumAxis, 0) +} + +// reduceAxisOr is reduceAxis for a reduction with an identity: along a +// zero-length axis every output element is the identity, as in numpy +// (np.zeros((3, 0)).sum(axis=1) is [0, 0, 0], .prod(axis=1) is [1, 1, 1]), +// and as Sum and Prod already return for an empty array. +func (a *Array) reduceAxisOr( + axis int, keepdims bool, + kernel func(dst, src []float64, outer, axisLen, inner, lo, hi int), + identity float64, +) (*Array, error) { + ax, err := a.normalizeAxis(axis) + if err != nil { + return nil, err + } + if a.shape[ax] != 0 { + return a.reduceAxis(ax, keepdims, kernel) + } + shape := a.reduceShape(ax, keepdims) + dst := a.alloc(prod(shape), false) + for i := range dst { + dst[i] = identity + } + return &Array{data: dst, shape: shape, strides: rowMajorStrides(shape), ws: a.ws}, nil } // ProdAxis returns the product along the given axis. See SumAxis for the // axis/keepdims semantics. func (a *Array) ProdAxis(axis int, keepdims bool) (*Array, error) { - return a.reduceAxis(axis, keepdims, kernels.ProdAxis) + return a.reduceAxisOr(axis, keepdims, kernels.ProdAxis, 1) } // MaxAxis returns the maximum along the given axis. See SumAxis for the @@ -910,13 +934,17 @@ func (a *Array) MinAxis(axis int, keepdims bool) (*Array, error) { // MeanAxis returns the arithmetic mean along the given axis. See SumAxis for the // axis/keepdims semantics. func (a *Array) MeanAxis(axis int, keepdims bool) (*Array, error) { - r, err := a.SumAxis(axis, keepdims) + ax, err := a.normalizeAxis(axis) if err != nil { return nil, err } - // The reduced axis length is positive (SumAxis rejected zero-length axes), - // so this scaling is well-defined. - ax, _ := a.normalizeAxis(axis) + // numpy returns NaN with a RuntimeWarning for the mean of an empty slice; + // Go has no warnings, so this is an error, as it is for Mean. + if a.shape[ax] == 0 { + return nil, fmt.Errorf("%w: mean along zero-length axis %d", + ErrShapeMismatch, ax) + } + r, _ := a.SumAxis(ax, keepdims) // ax is valid and non-empty: cannot fail n := float64(a.shape[ax]) for i := range r.data { r.data[i] /= n diff --git a/ndarray_test.go b/ndarray_test.go index 3e7191f..fd0eaad 100644 --- a/ndarray_test.go +++ b/ndarray_test.go @@ -623,14 +623,35 @@ func TestAxisReductionErrors(t *testing.T) { } } - // Reduction along a zero-length axis is rejected (empty-reduction case), - // for both the kernel-backed reductions and the Sum-derived Mean. + // Along a zero-length axis, as numpy 2.5.3 does: Sum and Prod give their + // identity (np.zeros((0, 3)).sum(axis=0) is [0, 0, 0]); Max, Min, ArgMax + // and ArgMin have none and are rejected (numpy raises ValueError); Mean is + // rejected too, where numpy returns NaN with a RuntimeWarning. z := mustArr(t, ok(New(0, 3))) - if _, err := z.SumAxis(0, false); !errors.Is(err, ErrShapeMismatch) { - t.Fatalf("zero-axis Sum: %v", err) + for name, c := range map[string]struct { + f func(int, bool) (*Array, error) + want float64 + }{"Sum": {z.SumAxis, 0}, "Prod": {z.ProdAxis, 1}} { + r := mustArr(t, ok(c.f(0, false))) + eqInts(t, r.Shape(), []int{3}) + for i := 0; i < 3; i++ { + if r.At(i) != c.want { + t.Fatalf("zero-axis %s[%d] = %v, want %v", name, i, r.At(i), c.want) + } + } + rk := mustArr(t, ok(c.f(0, true))) + eqInts(t, rk.Shape(), []int{1, 3}) + } + for name, f := range map[string]func(int, bool) (*Array, error){ + "Mean": z.MeanAxis, "Max": z.MaxAxis, "Min": z.MinAxis, + "ArgMax": z.ArgMaxAxis, "ArgMin": z.ArgMinAxis, + } { + if _, err := f(0, false); !errors.Is(err, ErrShapeMismatch) { + t.Fatalf("zero-axis %s: %v", name, err) + } } - if _, err := z.MeanAxis(0, false); !errors.Is(err, ErrShapeMismatch) { - t.Fatalf("zero-axis Mean: %v", err) + if _, err := z.MeanAxis(9, false); !errors.Is(err, ErrAxis) { + t.Fatalf("MeanAxis bad axis: %v", err) } // Reducing a non-zero axis of an array that still has zero elements keeps an // empty result rather than erroring (axisLen > 0 but inner produces nothing). diff --git a/result_shape_test.go b/result_shape_test.go index 429ff96..4813e8e 100644 --- a/result_shape_test.go +++ b/result_shape_test.go @@ -70,7 +70,14 @@ func TestEmptyArraysWithHugeAxesReturnPromptly(t *testing.T) { ops := []func() (*Array, error){ func() (*Array, error) { return Concatenate([]*Array{a, a}, -1) }, func() (*Array, error) { return Concatenate([]*Array{a, a}, 0) }, - func() (*Array, error) { return a.SumAxis(0, false) }, + func() (*Array, error) { + if sh[0] == 0 { + // Summing the empty axis of (0, H) is H zeros, as in numpy: a + // real allocation of H elements, not a walk over nothing. + return nil, nil + } + return a.SumAxis(0, false) + }, func() (*Array, error) { return a.MaxAxis(0, true) }, func() (*Array, error) { return a.MeanAxis(0, false) }, func() (*Array, error) { return a.ArgMaxAxis(0, false) }, @@ -81,7 +88,7 @@ func TestEmptyArraysWithHugeAxesReturnPromptly(t *testing.T) { go func() { defer close(done) for _, op := range ops { - if r, err := op(); err == nil && r.Size() != 0 { + if r, err := op(); err == nil && r != nil && r.Size() != 0 { t.Errorf("%v: a result with %d elements from an empty array", sh, r.Size()) } } diff --git a/slicing.go b/slicing.go index 6138b76..f650d79 100644 --- a/slicing.go +++ b/slicing.go @@ -143,16 +143,21 @@ func normBound(b, length, lo, hi int) int { } // Slice returns a view of the array selected by per-axis indices, following -// NumPy basic-indexing semantics. It accepts exactly Ndim index arguments -// (build them with A, All, R, Rng, From, To). Integer indices (A) drop their +// NumPy basic-indexing semantics. It accepts up to Ndim index arguments (build +// them with A, All, R, Rng, From, To, Step); as in NumPy, the axes left +// unindexed are taken whole, so a.Slice(A(1)) of a matrix is its second row +// and a.Slice() is a view of all of a. Integer indices (A) drop their // axis; range indices keep the axis as a strided view that shares the receiver's // backing data — writes through the view are visible in the original and vice // versa. The result may have a non-zero offset and arbitrary strides. func (a *Array) Slice(idx ...Index) (*Array, error) { - if len(idx) != len(a.shape) { + if len(idx) > len(a.shape) { return nil, fmt.Errorf("%w: got %d indices for %d-dimensional array", ErrIndex, len(idx), len(a.shape)) } + for len(idx) < len(a.shape) { + idx = append(idx[:len(idx):len(idx)], All()) // never into the caller's array + } offset := a.offset shape := make([]int, 0, len(a.shape)) strides := make([]int, 0, len(a.shape)) diff --git a/slicing_test.go b/slicing_test.go index 8747644..368f4bf 100644 --- a/slicing_test.go +++ b/slicing_test.go @@ -119,8 +119,8 @@ func TestSliceIsViewWriteThrough(t *testing.T) { func TestSliceErrors(t *testing.T) { m := mustArr(t, ok(FromData([]float64{1, 2, 3, 4}, 2, 2))) - if _, err := m.Slice(A(0)); !errors.Is(err, ErrIndex) { - t.Fatalf("too few indices: %v", err) + if _, err := m.Slice(A(0), All(), All()); !errors.Is(err, ErrIndex) { + t.Fatalf("too many indices: %v", err) } if _, err := m.Slice(A(5), All()); !errors.Is(err, ErrIndex) { t.Fatalf("int index out of range: %v", err) @@ -132,3 +132,28 @@ func TestSliceErrors(t *testing.T) { t.Fatalf("zero step: %v", err) } } + +// As in NumPy, unindexed trailing axes are taken whole: m[1] is a row and m[()] +// is all of m. Filling them in must not write into the caller's index slice. +func TestSliceFewerIndices(t *testing.T) { + m := mustArr(t, ok(FromData([]float64{1, 2, 3, 4, 5, 6}, 2, 3))) + row := mustArr(t, ok(m.Slice(A(1)))) + eqInts(t, row.Shape(), []int{3}) + if row.At(0) != 4 || row.At(2) != 6 { + t.Fatalf("m[1] = %v", row) + } + all := mustArr(t, ok(m.Slice())) + eqInts(t, all.Shape(), []int{2, 3}) + if all.At(1, 2) != 6 { + t.Fatalf("m[()] = %v", all) + } + backing := []Index{Step(-1), A(0)} + rev := mustArr(t, ok(m.Slice(backing[:1]...))) + eqInts(t, rev.Shape(), []int{2, 3}) + if rev.At(0, 0) != 4 { + t.Fatalf("m[::-1][0,0] = %v", rev.At(0, 0)) + } + if backing[1] != A(0) { + t.Fatalf("Slice wrote into the caller's slice: %+v", backing[1]) + } +} diff --git a/testdata/fuzz/FuzzOps/404c93efded6ff79 b/testdata/fuzz/FuzzOps/404c93efded6ff79 new file mode 100644 index 0000000..19b59c9 --- /dev/null +++ b/testdata/fuzz/FuzzOps/404c93efded6ff79 @@ -0,0 +1,2 @@ +go test fuzz v1 +[]byte("7AAA0000000000000000000000000000000000000000000000000000000000010000000000081200700000090;")