Evaluate array at specific subarray

Question

I warn in advance: I may be utterly confused at the moment. I tell a short story about what I actually try to achieve because that may clear things up. Say I have f(a,b,c,d,e), and I want to find arg max (d,e) f(a,b,c,d,e). Consider a (trivial example of a ) discretized grid F of f:

F = np.tile(np.arange(0,10,0.1)[newaxis,newaxis,:,newaxis,newaxis], [10, 10, 1, 10, 10])
maxE = F.max(axis=-1)
argmaxD = maxE.argmax(axis=-1)
maxD = F.max(axis=-2)
argmaxE = maxD.argmax(axis=-1)

This is the case how I typically solve the discretized version. But now assume instead, that I want to solve arg max d f(a,b,c,d,e=X): Instead of optimally chosen e for every other input, e is a fixed and given (of size AxBxCxD, which in this example would be 10x10x100x10). I have troubles solving this.

My naive approach was

X = np.tile(np.arange(0,10)[newaxis,newaxis,:,newaxis], [10,10,1,10])
maxX = F[X]
argmaxD = maxX.argmax(axis=-1)

However, the huge surge of memory that crashes my IDE implies that F[X] is apparently not what I was looking for.

Performance is key.

it seems what you want is something like: np.argmax(np.max(F, axis=-1), axis=-1) — Saullo G. P. Castro
– Saullo G. P. Castro, Commented Jul 13, 2014 at 17:18

score 2 · Accepted Answer · 2014-07-17 15:25:45Z

2

+50

I believe you can do it like this, but maybe there's a better way..

n = 10
F = np.tile(np.arange(0,n,0.1)[None,None,:,None,None], [n, n, 1, n, n])
X = np.tile(np.arange(0,n)[None,None,:,None], [n, n, 1, n])

a,b,c,d = np.ogrid[:n,:n,:n,:n]
argmaxD = F[a,b,c,d,X].argmax(axis=-1)

Above X doesn't occupy the whole space, as we discussed in the comments. If you would like to choose e for all a,b,c and d you could do e.g.:

X = np.tile(np.arange(0,n,0.1).astype(int)[None,None,:,None], [n, n, 1, n])
a,b,c,d = np.ogrid[:n,:n,:100,:n]
argmaxD = F[a,b,c,d,X].argmax(axis=-1)

Also, notice that instead of tile you could make use of broadcasting. But then F[a,b,c,d,X] has a singular dimension so you should provide something like axis=3:

X = np.arange(0,n,0.1).astype(int)[None,None,:,None]
a,b,c,d = np.ogrid[:n,:n,:100,:n]
argmaxD = F[a,b,c,d,X].argmax(axis=3)

edited Jul 17, 2014 at 15:25

answered Jul 8, 2014 at 8:24

user2379410

Sign up to request clarification or add additional context in comments.

5 Comments

FooBar Over a year ago

Without fully digesting what's happening here: my old argmaxD (from standard maximization) had shape (10L, 10L, 100L), while the argmaxD you provide has shape (10L, 10L, 10L). As it's supposed to be the optimal maximizer given (a,b,c), which have shape (10, 10, 100), this cannot be correct.

user2379410 Over a year ago

@FooBar, I had already noticed the different size in one of the dimensions, but that's the result of choosing X = np.tile(np.arange(0,n),.... I guess you should just input the 0.1 stepsize. Then you should also use np.ogrid[:n, :n, :100, :n]

FooBar Over a year ago

0.1 stepsize gives IndexError, since it needs to be integer. X needs to contain integer values for what we selected out of the grid for e. e contains 100 elements. That means, we need to select between (0, 99) for e, right? So then I tried X = np.tile(np.arange(0,100)...), but I get an IndexError... I know I'm making a mistake, I can't just spot it.

user2379410 Over a year ago

@FooBar, oh right they're supposed to be indices.. Well indices for e can only be 0-9 because that dimension is only length 10. I guess you need to choose X differently. Maybe floor(arange(0, 10, 0.1))?

FooBar Over a year ago

With X = np.tile(np.arange(0.1,10)[newaxis,newaxis,:,newaxis], [10,10,10,10]), I get shape (10,10,100,10). Then adjust ogrid[:n, :n, :100, :n], and the rest works. I'm still curious why this is so complicated. Especially compared to Matlab, where I started doing most of my matrix algebra.

MeisenBar · Accepted Answer · 2014-07-12 17:05:09Z

0

This would be my idea to solve this.

from itertools import product, starmap

f = lambda a,b,c,d,e : d / e

args_iterable = product([1],[2],[3],range(1,1000),range(1,1000))

max_val, max_args = max(starmap(lambda *args: (f(*args), args) , args_iterable))

print max_args

answered Jul 12, 2014 at 17:05

MeisenBar

547 bronze badges

Collectives™ on Stack Overflow

Evaluate array at specific subarray

2 Answers 2

5 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

5 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related