Slicing a NumPy array with another set of arrays having beginning and ending set of indices

Question

I am writing a code to simulate Continuous Time Random Walk phenomena with a function in python. My code works correctly so far, but I would like to exploit the indexing abilities of NumPy arrays and improve the speed. In the code, below I am generating an ensemble of trajectories, so I have to loop over each of them while generating it. Is it somehow possible to index the NumPy array x in such a way that I can get rid of the loop on Nens (the for loop in the below code snippet)

    for k in range(Nens):
            
        #start building the trajectory
        stop = 0
        i = 0
        while (stop < Nt):

            #reset the the stop time
            start = stop

            #increment the stop till the next waiting time
            stop += int(trand[i,k]) # A 2D numpy array

   
            #update the trajectory
            x[start:stop,k] = x[start-1,k] \ #x is also a 2D array
                            + (1-int(abs((x[start-1,k]+xrand[i,k])/(xmax))))* xrand[i,k] \
                            - int(abs((x[start-1,k]+xrand[i,k])/(xmax)))*np.sign(x[start-1,k]/xrand[i,k])* xrand[i,k] 
   
            i += 1

    print i
    return T, x

A plausible method that I can look to is as follows.

In this code, start and stop are scalar integers. However, I would like to index this is in way in which both start and stop are 1D Numpy integer arrays. But I have seen that if I can use only stop/start to slice the numpy array, but using slicing from a beginning to ending tuple of indices is not possible .

EDIT 1 (MWE): The following is the function that I have written, which produces random walk trajectory if given the appropriate parameters,

def ctrw_ens2d(sig,tau,sig2,tau2,alpha,xmax,Nens,Nt=1000,dt=1.0):
    #first define the array for time
    T = np.arange(0,Nt,1)*dt
    
    #generate at least Nt random time increments based on Poisson 
    # distribution (you will use only less than that)
    trand = np.random.exponential(tau, (2*Nt,Nens,1))
    xrand = np.random.normal(0.0,sig,(2*Nt,Nens,2))
    Xdist = np.random.lognormal(-1,0.9,(Nens))
    Xdist = np.clip(Xdist,2*sig,12*sig)
    
    trand2 = np.random.exponential(tau2, (2*Nt,Nens,1))
    xrand2 = np.random.normal(0.0,sig2,(2*Nt,Nens,2))
    
    
    #make a zero array of trajectory
    x1 = np.zeros((Nt,Nens))
    x2 = np.zeros((Nt,Nens))
    y1 = np.zeros((Nt,Nens))
    y2 = np.zeros((Nt,Nens))
    
    for k in range(Nens):
        #start building the trajectory
        stop = 0
        i = 0

        while (stop < Nt):

            #reset the the stop time
            start = stop

            #increment the stop till the next waiting time
            stop += int(trand[i,k,0])           

            #update the trajectory
            r1 = np.sqrt(x1[start-1,k]**2 + y1[start-1,k]**2)
            rr = np.linalg.norm(xrand[i,k])

            x1[start:stop,k] = x1[start-1,k] \
                               + (1-int(abs((r1+rr)/(Xdist[k]))))* xrand[i,k,0] \
                               - int(abs((r1+rr)/(Xdist[k])))* \
                               np.sign(x1[start-1,k]/xrand[i,k,0])* xrand[i,k,0] 
            
            y1[start:stop,k] = y1[start-1,k] \
                                 + (1-int(abs((r1+rr)/(Xdist[k]))))* xrand[i,k,1] \
                                - int(abs((r1+rr)/(Xdist[k])))* \
                                 np.sign(y1[start-1,k]/xrand[i,k,1])* xrand[i,k,1] 
            i += 1
        
        #randomly add jumps in between, at later stage
    
        stop = 1
        i = 0
        while (stop < Nt):

            #reset the the stop time
            start = stop

            #increment the stop till the next waiting time
            stop += int(trand2[i,k,0])

            #update the trajectory
            x2[start:stop,k] = x2[start-1,k] + xrand2[i,k,0]
            y2[start:stop,k] = y2[start-1,k] + xrand2[i,k,1]
            i += 1

        
    return T, (x1+x2), (y1+y2)

A simple run of the above function is given below,

    Tmin = 0.61 # in ps
    Tmax = 1000 # in ps
    NT = int(Tmax/Tmin)*10
    delt = (Tmax-0.0)/NT
    print "Delta T, No. of timesteps:",delt,NT
    
    Dint = 0.21 #give it Ang^2/ps
    sig = 0.3 #in Ang
    xmax = 5.*sig
    tau = sig**2/(2*Dint)/delt  # from ps, convert it into the required units according to delt
    print "Waiting time for confined motion (in Delta T units)",tau
    
    Dj = 0.03 # in Ang^2/ps
    tau2 = 10 # in ps
    sig2 = np.sqrt(2*Dj*tau2)
    print "Sigma 2:", sig2
    tau2 = tau2/delt
    alpha  = 1
    
    tim, xtall, ytall = ctrw_ens2d(sig,tau,sig2,tau2,alpha,xmax,100,Nt=NT,dt=delt)

The generated trajectories can be plotted as follows,

    rall = np.stack((xtall,ytall),axis=-1)
    print rall.shape
    print xtall.shape
    print rall[:,99,:].shape
    k = 19
    plt.plot(xtall[:,k],ytall[:,k])

it appears that the length of each slice can vary? Then you have use each separately. — hpaulj
– hpaulj, Commented Oct 8, 2021 at 8:31
Indeed, that is the main bottleneck in this problem. If not in NumPy, are there other data objects that can aid in this kind of a slicing? — user35952
– user35952, Commented Oct 8, 2021 at 8:48
Not every loop can be NumPy-vectorized. Since your motivation is speed/efficiency, maybe you can try to Numba-compile the loop. You didn't post a MWE so we can't try it out for you. — BatWannaBe
– BatWannaBe, Commented Oct 12, 2021 at 8:45
@BatWannaBe I have not used Numba so far, please have a look at the MWE that I have added and let me know about it's feasibility. — user35952
– user35952, Commented Oct 13, 2021 at 5:17
You can (and should) try it yourself: import numba as nb and decorate functions with @nb.njit. Be sure to use njit: it's an alias for jit(nopython=True), and compiling with "no Python" is how you get speed. But it's tricky: Numba does not inherently support all Python or NumPy functions; it has to re-implement things for compilation. @nb.njit is handy because it'll throw an error if there's something Numba can't compile. So what you would do is try to carve out the part of your function that is njit-compatible and test its speed (no guarantees of speedup). — BatWannaBe
– BatWannaBe, Commented Oct 13, 2021 at 6:39

Bob · Accepted Answer · 2021-10-13 21:16:16Z

Starting with a zero array, the loop

while stop < Nt:
  start = stop
  stop += randint();
  x[start:stop] = x[start-1] + rand();

will create a series of steps. A step can be achieved with the cumulative sum of the inpulse

while stop < Nt:
  start = stop
  stop += randint();
  x[start] = any();
np.cumsum(x, out=x)

This applies to both the first and second loop. The (x2, y2) are more easily vectorized because the increments do not depend on the previous values The (x2, y2) still require a while loop, but each iteration can be vectorized.

The final result is like this

def ctrw_ens2d_vectorized(sig,tau,sig2,tau2,alpha,xmax,Nens,Nt=1000,dt=1.0):
    # first define the array for time
    T = np.arange(0,Nt,1)*dt

    # generate at least Nt random time increments based on Poisson 
    # distribution (you will use only less than that)
    trand = np.random.exponential(tau, (2*Nt,Nens,1))
    xrand = np.random.normal(0.0,sig,(2*Nt,Nens,2))
    Xdist = np.random.lognormal(-1,0.9,(Nens))
    Xdist = np.clip(Xdist,2*sig,12*sig)

    trand2 = np.random.exponential(tau2, (2*Nt,Nens,1)).astype(np.int64)
    xrand2 = np.random.normal(0.0,sig2,(2*Nt,Nens,2))


    #make a zero array of trajectory
    x1 = np.zeros((Nt,Nens))
    x2 = np.zeros((Nt,Nens))
    y1 = np.zeros((Nt,Nens))
    y2 = np.zeros((Nt,Nens))
    
    #randomly add jumps in between, at later stage
    stop = 1 + np.cumsum(trand2[:,:,0], axis=0)
    # vectorize the indices
    I, J = np.indices(stop.shape)
    m = stop < Nt # Vectorized loop stopping condition
    I = I[m]; J = J[m]; # indices only for the executed iterations
    # update x
    x2[stop[I,J], J] = xrand2[I,J,0]
    y2[stop[I,J], J] = xrand2[I,J,1]
    np.cumsum(x2, axis=0, out=x2)
    np.cumsum(y2, axis=0, out=y2)
    
    
    # this part is more complicated and I vectorized on axis 1
    stop = np.zeros(Nens, dtype=np.int64)
    start = np.zeros(Nens, dtype=np.int64)
    k = np.arange(Nens)
    i = 0
    zx1 = np.zeros_like(x1[0])
    zy1 = np.zeros_like(y1[0])
    assert(np.all(trand > 0))
    m = k
    i = 0
    while np.any(stop < Nt):
        start[:] = stop;
        stop[m] += trand[i,m,0].astype(np.int64)
        m = k[stop < Nt];
        r1 = np.sqrt(zx1[m]**2 + zy1[m]**2)
        rr = np.linalg.norm(xrand[i,m,:],axis=-1) # axis requires numpy 1.8
        
        tx = (1-(abs((r1+rr)/(Xdist[m]))).astype(np.int64))* xrand[i,m,0] \
                - (abs((r1+rr)/(Xdist[m]))).astype(np.int64)* \
                 np.sign(zx1[m]/xrand[i,m,0])* xrand[i,m,0]
        ty = (1-(abs((r1+rr)/(Xdist[m]))).astype(np.int64))* xrand[i,m,1] \
                - (abs((r1+rr)/(Xdist[m]))).astype(np.int64)* \
                 np.sign(zy1[m]/xrand[i,m,1])* xrand[i,m,1] 
        zx1[m] += tx[:] * (start[m] < stop[m])
        zy1[m] += ty[:] * (start[m] < stop[m])
    
        x1[start[m],m] = tx[:]
        y1[start[m],m] = ty[:]   
        i += 1
    
    np.cumsum(x1, axis=0, out=x1)
    np.cumsum(y1, axis=0, out=y1)
    
    return T, (x1+x2), (y1+y2)

This runs ~8x faster than the original code here.

Collectives™ on Stack Overflow

Slicing a NumPy array with another set of arrays having beginning and ending set of indices

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related