# Integer-array indexing (e.g. numpy.take) as function composition: where can I find more resources?

For context, one of Numpy's features is that an an array of integers can be passed as an index to an array, and this selects values at all of the specified positions, in any order with possible duplicates.

>>> import numpy
>>> a = numpy.array([0.0, 1.1, 2.2, 3.3, 4.4, 5.5])
>>> b = numpy.array([2, 4, 3, 3, 0])
>>> a[b]
array([2.2, 4.4, 3.3, 3.3, 0. ])

>>> # alternative syntax
>>> numpy.take(a, b)
array([2.2, 4.4, 3.3, 3.3, 0. ])

I think this is a feature of most array libraries/languages, as well as SIMD instruction sets (specifically "gather").

If we consider the array as a function from a single integer index to its data type (i.e. $$[0, N) \to \mbox{dtype}$$ for an array of length $$N$$), then indexing with an array of integers is function composition. This, coupled with the fact that function composition is associative has been incredibly useful in my Numpy scripts and I'd like to know if there's a theoretical literature that might go into it in more depth.

To demonstrate what I mean with an example, consider any two functions from and to non-negative integers ($$\mathbb{Z}^{\ge 0} \to \mathbb{Z}^{\ge 0}$$):

def f(x):
return x**2 - 5*x + 10
def g(y):
return max(0, 2*y - 10) + 3

Sample these functions and their composition $$g(f(\cdot))$$ at the first few integers in their domain:

F   = numpy.array([f(i) for i in range(10)])     # F is f at 10 elements
G   = numpy.array([g(i) for i in range(100)])    # G is g at 100 elements (enough)
GoF = numpy.array([g(f(i)) for i in range(10)])  # GoF is g∘f at 10 elements

Now array-indexing G[F] returns the same result as the sampled composition GoF:

print("G\u2218F =", G[F])    # integer indexing
print("g\u2218f =", GoF)     # array of the composed functions
G∘F = [13  5  3  3  5 13 25 41 61 85]
g∘f = [13  5  3  3  5 13 25 41 61 85]

Since integer-array indexing is function composition, it's associative like function composition:

H = numpy.arange(1000)*1.1
print("H[G][F] =", H[G][F])  # index H by G, then the result by F
print("H[G[F]] =", H[G[F]])  # index H by the result of indexing G by F
H[G][F] = [14.3  5.5  3.3  3.3  5.5 14.3 27.5 45.1 67.1 93.5]
H[G[F]] = [14.3  5.5  3.3  3.3  5.5 14.3 27.5 45.1 67.1 93.5]

and this fact can be very useful (e.g. as an optimization if you have many H's or as a simplification if you're recursing over a tree of arrays). In some applications, it borders on magical.

Is there a paper or a sub-subbranch of computer science where this sort of thing is developed in more detail? What should I read to learn more?

• Welcome to CS.SE! I appreciate the thoughtful post and all of the background. I wanted to prepare you for the possibility that, unfortunately, your question might not fare so well here: in my experience the site tends to work better for more focused questions. If you can identify a more specific question about this, it's possible that might be more likely to get an answer. (For instance, you might ponder what you'd do if you found such a resource -- is there some specific task or problem you'd hope it'd teach you how to solve?) Hope someone can give you useful information! – D.W. Sep 16 '19 at 18:35
• Should I perhaps move it to cstheory.stackexchange.com? (The questions I see there seem to be deeper/more advanced, though...) – Jim Pivarski Sep 16 '19 at 22:05
• en.wikipedia.org/wiki/Map_(higher-order_function) – Yuval Filmus Sep 19 '19 at 13:37