Indexing routines: implement the 25 missing np.* indexing functions (diag, nditer, r_/c_, take_along_axis, ...)

[Nucs]

Overview

Indexing routines is NumSharp's largest partial gap: 10 of 35 implemented. This tracks the 25 missing functions, grouped into four work clusters. It is dominated by np.diag -- ~115K GitHub hits, one of the most-used NumPy functions overall and absent today (NumSharp has diagonal for extraction but not the diag constructor/dual).

Problem

The gap spans four very different kinds of work -- pure index arithmetic (diagonals/triangles), a slice-expression DSL (r_/c_/s_/ix_/ogrid), iterator objects, and along-axis selection. Splitting by cluster lets the cheap, high-value pieces (esp. diag) ship first.

Verified codebase state

• Implemented (np.*): take, put, compress, place, where, nonzero, diagonal, indices, ravel_multi_index, unravel_index (+ mgrid exists).
• Missing (25 below): verified absent from np.* (including static properties, so c_/r_/s_/ogrid are confirmed missing) and as NDArray methods (reflection + source grep).
• Nuance: the internal multi-operand iterator (public static class NDIter, Backends/Iterators/NDIter.cs) already exists -- so the nditer/ndindex/ndenumerate cluster is mostly a public wrapper, not new machinery.

Proposal

Add under np.*, cluster by cluster. Recommended order: A -> (D, C) -> B, leading with diag.

A. Diagonal & triangle index-math -- pure index arithmetic, no new kernels (do first)

• diag (114,816) -- extract a diagonal or construct a diagonal matrix (the dual of diagonal)
• fill_diagonal (26,032) -- write the main diagonal of an array in place
• triu_indices (8,944) -- upper-triangle indices of an (n,m) array
• triu_indices_from (4,600) -- upper-triangle indices for a given array
• diag_indices (4,340) -- indices to access the main diagonal
• tril_indices (3,800) -- lower-triangle indices of an (n,m) array
• diag_indices_from (3,488) -- main-diagonal indices for a given array
• tril_indices_from (3,140) -- lower-triangle indices for a given array
• mask_indices (2,156) -- indices from a generic triangle-mask function

D. Along-axis selection -- selection kernels (build on take/put/place)

• select (23,040) -- piecewise pick from choice arrays by a condition list
• take_along_axis (10,584) -- gather by matching 1-D index and data slices
• choose (9,416) -- build an array from an index array + a list of choice arrays
• putmask (3,512) -- masked in-place write (broadcasting values)
• put_along_axis (3,088) -- scatter by matching 1-D index and data slices

C. Iterator objects -- mostly surfacing the existing NDIter engine

• nditer (42,416) -- multi-operand iterator object -- internal NDIter engine already exists (Backends/Iterators/NDIter.cs)
• ndindex (20,160) -- N-dimensional index iterator
• ndenumerate (8,288) -- N-dimensional (index, value) iterator
• iterable (5,488) -- predicate: is the object iterable (trivial)
• nested_iters (1,664) -- nested multi-operand iterators
• flatiter (4,800) -- flat-iterator type behind .flat

B. Grid / slice-expression DSL -- needs C# indexer sugar

• r_ (44,016) -- row-wise concatenation / build from slice objects
• c_ (35,600) -- column-wise concatenation / build from slice objects
• ix_ (32,544) -- open-mesh index from multiple sequences
• s_ (19,104) -- slice-literal builder for index tuples
• ogrid (6,888) -- open (non-fleshed) meshgrid (mgrid already exists in NumSharp)

Evidence -- GitHub usage ranking

#	Function	GitHub hits	Cluster	What it does
1	np.diag	114,816	A	extract a diagonal or construct a diagonal matrix (the dual of diagonal)
2	np.r_[]	44,016	B	row-wise concatenation / build from slice objects
3	np.nditer	42,416	C	multi-operand iterator object -- internal NDIter engine already exists (Backends/Iterators/NDIter.cs)
4	np.c_[]	35,600	B	column-wise concatenation / build from slice objects
5	np.ix_	32,544	B	open-mesh index from multiple sequences
6	np.fill_diagonal	26,032	A	write the main diagonal of an array in place
7	np.select	23,040	D	piecewise pick from choice arrays by a condition list
8	np.ndindex	20,160	C	N-dimensional index iterator
9	np.s_[]	19,104	B	slice-literal builder for index tuples
10	np.take_along_axis	10,584	D	gather by matching 1-D index and data slices
11	np.choose	9,416	D	build an array from an index array + a list of choice arrays
12	np.triu_indices	8,944	A	upper-triangle indices of an (n,m) array
13	np.ndenumerate	8,288	C	N-dimensional (index, value) iterator
14	np.ogrid[]	6,888	B	open (non-fleshed) meshgrid (mgrid already exists in NumSharp)
15	np.iterable	5,488	C	predicate: is the object iterable (trivial)
16	np.flatiter	4,800	C	flat-iterator type behind .flat
17	np.triu_indices_from	4,600	A	upper-triangle indices for a given array
18	np.diag_indices	4,340	A	indices to access the main diagonal
19	np.tril_indices	3,800	A	lower-triangle indices of an (n,m) array
20	np.putmask	3,512	D	masked in-place write (broadcasting values)
21	np.diag_indices_from	3,488	A	main-diagonal indices for a given array
22	np.tril_indices_from	3,140	A	lower-triangle indices for a given array
23	np.put_along_axis	3,088	D	scatter by matching 1-D index and data slices
24	np.mask_indices	2,156	A	indices from a generic triangle-mask function
25	np.nested_iters	1,664	C	nested multi-operand iterators

Combined usage ~ 441,924 hits -- the largest of the three sweep issues, and diag (~114,816) is the single most-used missing function in the entire sweep.

_{Usage = GitHub REST code-search total_count in each function's canonical call form (np.<fn>(, or np.<fn>[ for the index-expression objects). It is an approximate estimate (noisy plus/minus ~10% on large sets) and counts only the dominant np. alias. Measured 2026-07-19; treat as relative demand, not exact.}

Scope / Non-goals

• np.matrix-returning behavior is out of scope (NumSharp has no matrix subclass; np.matlib is deprecated).
• nditer should expose the existing NDIter engine, not a re-implementation; full NumPy op_flags/buffering parity can be staged.
• The DSL objects (r_/c_/s_) require C# indexer ergonomics -- acceptable to expose a close-but-not-identical surface where C# syntax forces it (document divergences).

Related issues

• Part of the NumPy category-coverage sweep; sibling issues: Set routines: implement np.* set operations (isin, intersect1d, union1d, setdiff1d, setxor1d, unique_*) #622 (set) and Sorting & searching: implement partition/argpartition, lexsort, nanargmax/nanargmin, sort_complex #623 (sort). Umbrella / methodology: LinearAlgebra Missing Support #621 (linalg).

[Nucs]

Checked off the 16 functions verified implemented in src/NumSharp.Core (as of 297f883f):

• Cluster A — diagonal/triangle math: diag, fill_diagonal, triu_indices(+_from), tril_indices(+_from), diag_indices(+_from), mask_indices — Indexing/np.{diag,fill_diagonal,tril_indices,diag_indices}.cs.
• Cluster C — iterator objects: nditer, ndindex, ndenumerate — APIs/np.nditer.cs, Indexing/np.{ndindex,ndenumerate}.cs (plus the NumSharp-extension typed nditer<T>/nditer_chunks<T>).
• Cluster B — slice-expression DSL: r_, c_, ix_, s_ — Creation/np.{r_,c_}.cs, Indexing/np.{ix_,s_}.cs (index_exp also landed).

Still open (9):

• Cluster D (selection): select, take_along_axis, choose, putmask, put_along_axis — no selection kernels yet.
• Cluster C remainder: iterable, flatiter, nested_iters (.flat returns a raveled NDArray, not a flatiter type).
• Cluster B remainder: ogrid (mgrid exists).