python - Gridwise application of the bisection method -

i need find roots generalized state space. is, have discrete grid of dimensions grid=axbx(...)xx, of not know ex ante how many dimensions has (the solution should applicable grid.size) .

i want find roots (f(z) = 0) every state z inside grid using bisection method. remainder contains f(z), , know f'(z) < 0. need

• increase z if remainder > 0
• decrease z if remainder < 0

wlog, matrix historyof shape (grid.shape, t) contains history of earlier values of z every point in grid , need increase z (since remainder > 0). need select zalternative inside history[z, :] "smallest of those, larger z". in pseudo-code, is:

zalternative =  hist[z,:][hist[z,:] > z].min() 

i had asked earlier. solution given

b = sort(history[..., :-1], axis=-1) mask = b > history[..., -1:] index = argmax(mask, axis=-1) indices = tuple([arange(j) j in b.shape[:-1]]) indices = meshgrid(*indices, indexing='ij', sparse=true) indices.append(index) indices = tuple(indices) lowerz = history[indices]  b = sort(history[..., :-1], axis=-1) mask = b <= history[..., -1:] index = argmax(mask, axis=-1) indices = tuple([arange(j) j in b.shape[:-1]]) indices = meshgrid(*indices, indexing='ij', sparse=true) indices.append(index) indices = tuple(indices) higherz = history[indices]  newz = history[..., -1] criterion = 0.05 increase = remainder > 0 + criterion decrease = remainder < 0 - criterion newz[increase] = 0.5*(newz[increase] + higherz[increase]) newz[decrease] = 0.5*(newz[decrease] + lowerz[decrease]) 

however, code ceases work me. feel extremely bad admitting it, never understood magic happening indices, therefore unfortunately need help.

what code actually does, give me lowest respectively highest. is, if fix on 2 specific z values:

history[z1] = array([0.3, 0.2, 0.1]) history[z2] = array([0.1, 0.2, 0.3]) 

i higherz[z1] = 0.3 , lowerz[z2] = 0.1, is, extrema. correct value both cases have been 0.2. what's going wrong here?

if needed, in order generate testing data, can use along lines of

history = tile(array([0.1, 0.3, 0.2, 0.15, 0.13])[newaxis,newaxis,:], (10, 20, 1)) remainder = -1*ones((10, 20)) 

to test second case.

expected outcome

i adjusted history variable above, give test cases both upwards , downwards. expected outcome

lowerz = 0.1 * ones((10,20)) higherz = 0.15 * ones((10,20)) 

which is, every point z in history[z, :], next highest previous value (higherz) , next smallest previous value (lowerz). since points z have same history ([0.1, 0.3, 0.2, 0.15, 0.13]), have same values lowerz , higherz. of course, in general, histories each z different , hence 2 matrices contain potentially different values on every grid point.

i compared posted here the solution previous post , noticed differences.

for smaller z, said

mask = b > history[..., -1:] index = argmax(mask, axis=-1) 

they said:

mask = b >= a[..., -1:] index = np.argmax(mask, axis=-1) - 1 

for larger z, said

mask = b <= history[..., -1:] index = argmax(mask, axis=-1) 

they said:

mask = b > a[..., -1:] index = np.argmax(mask, axis=-1) 

using the solution previous post, get:

import numpy np history = np.tile(np.array([0.1, 0.3, 0.2, 0.15, 0.13])[np.newaxis,np.newaxis,:], (10, 20, 1)) remainder = -1*np.ones((10, 20))  = history  # b sorted ndarray excluding recent observation # sorted along observation axis b = np.sort(a[..., :-1], axis=-1)  # mask boolean array, comparing (sorted) # previous observations current observation - [..., -1:] mask = b > a[..., -1:]  # next 5 statements build indexing array. # true evaluates 1 , false evaluates zero. # argmax() return index of first true, # in case along last (observations) axis. # index array shape of z (2-d test data). # represents index of next greater # observation every 'element' of z. index = np.argmax(mask, axis=-1)  # next 2 statements construct arrays of indices # every element of z - first n-1 dimensions of history. indices = tuple([np.arange(j) j in b.shape[:-1]]) indices = np.meshgrid(*indices, indexing='ij', sparse=true)  # adding index end of indices (the last dimension of history) # produces 'group' of indices 'select' single observation # every 'element' of z indices.append(index) indices = tuple(indices) higherz = b[indices]  mask = b >= a[..., -1:] # since b excludes current observation, want # index before next highest observation lowerz, # hence minus one. index = np.argmax(mask, axis=-1) - 1 indices = tuple([np.arange(j) j in b.shape[:-1]]) indices = np.meshgrid(*indices, indexing='ij', sparse=true) indices.append(index) indices = tuple(indices) lowerz = b[indices] assert np.all(lowerz == .1) assert np.all(higherz == .15) 

which seems work