Splitting integrated probability density into two spatial regions

Question

I have some probability density function:

T = 10000
tmin = 0
tmax = 10**20
t = np.linspace(tmin, tmax, T) 
time = np.asarray(t)                 #this line may be redundant 
for j in range(T):
    timedep_PD[j]= probdensity_func(x,time[j],initial_state)

I want to integrate it over two distinct regions of x. I tried the following to split the timedep_PD array into two spatial regions and then proceeded to integrate:

step =  abs(xmin - xmax) / T
l1 = int(np.floor((abs(ab - xmin)* T ) / abs(xmin - xmax)))
l2 = int(np.floor((abs(bd - ab)* T ) / abs(xmin - xmax)))

#For spatial region 1
R1 = np.empty([l1])
R1 = x[:l1]
for i in range(T):
    Pd1[i] = Pd[i][:l1]

#For spatial region 2
Pd2 = np.empty([T,l2])
R2 = np.empty([l2])
R2 = x[l1:l1+l2]
for i in range(T):
    Pd2[i] = Pd[i][l1:l1+l2]

#Integrating over each spatial region 
for i in range(T):
    P[0][i]   = np.trapz(Pd1[i],R1) 
    P[1][i]   = np.trapz(Pd2[i],R2) 

Is there an easier/more clear way to go about splitting up a probability density function into two spatial regions and then integrating within each spatial region at each time-step?

Answer 1

The loops can be eliminated by using vectorized operations instead. It's not clear whether Pd is a 2D NumPy array; it it's something else (e.g., a list of lists), it should be converted to a 2D NumPy array with np.array(...). After that you can do this:

Pd1 = Pd[:, :l1]
Pd2 = Pd[:, l1:l1+l2]

No need to loop over the time index; the slicing happens for all times at once (having : in place of an index means "all valid indices").

Similarly, np.trapz can integrate all time slices at once:

P1 = np.trapz(Pd1, R1, axis=1)
P2 = np.trapz(Pd2, R2, axis=1)

Each P1 and P2 is now a time series of integrals. The axis parameter determines along which axis Pd1 gets integrated - it's the second axis, i.e., space.