Calculation - numpy python bug

I am using NumPy to do some calculation on finding the Y intercept, through an Aperture between a big box and a small box. I have over 100.000 particles in the big box, and around 1000 in the small one. And it's taking a lot of time to do so. All the self.YD, self.XD are very large arrays that i'm multiplying.

PS: The ind are indexes of the values that need to be multiplied. I had a nonzero condition before that line in my code.

Any ideas how I would do this calculation in a simpler way?

YD_zero = self.oldYD[ind] - ((self.oldYD[ind]-self.YD[ind]) * self.oldXD[ind])/(self.oldXD[ind]-self.XD[ind])

Thanks!

UPDATE

Would using multiply, divide, subtract and all that stuff of Numpy. make it faster? Or if maybe if i split the开发者_如何学运维 calculation. for example.

to do this first:

    YD_zero = self.oldYD[ind] - ((self.oldYD[ind]-self.YD[ind])*self.oldXD[ind])

and then the next line would be:

    YD_zero /= (self.oldXD[ind]-self.XD[ind])

Any suggestions?!

UPDATE 2

I have been trying to figure this out, in a while now, but not much progress. My concern is that the denominator :

    self.oldXL[ind]-self.XL[ind] == 0

and I am getting some weird results.

The other thing is the nonzero function. I have been testing it for a while now. Could anybody tell me that it is almost the same as find in Matlab

---

Perhaps I have got the wrong end of the stick but in Numpy you can perform vectorised calculations. Remove the enclosing while loop and just run this ...

YD_zero = self.oldYD - ((self.oldYD - self.YD) * self.oldXD) / (self.oldXD - self.XD)

It should be much faster.

Update: Iterative root finding using the Newton-Raphson method ...

unconverged_mask = np.abs(f(y_vals)) > CONVERGENCE_VALUE:
while np.any(unconverged_mask):
    y_vals[unconverged_mask] = y_vals[unconverged_mask] - f(y_vals[unconverged_mask]) / f_prime(y_vals[unconverged_mask])
    unconverged_mask = np.abs(f(y_vals)) > CONVERGENCE_VALUE:

This code is only illustrative but it shows how you can apply an iterative process using vectorised code to any function f which you can find the derivative of f_prime. The unconverged_mask means that the results of the current iteration will only be applied to those values that have not yet converged.

Note that in this case there is no need to iterate, Newton-Raphson will give you the correct answer in the first iteration since we are dealing with straight lines. What you have is an exact solution.

Second update

Ok, so you aren't using Newton-Raphson. To calculate YD_zero (the y intercept) in one go, you can use,

YD_zero = YD + (XD - X0) * dYdX

where dYdX is the gradient, which seems to be, in your case,

dYdX = (YD - oldYD) / (XD - oldXD)

I am assuming XD and YD are the current x,y values of the particle, oldXD and oldYD are the previous x,y values of the particle and X0 is the x value of the aperture.

Still not entirely clear why you have to iterate over all the particles, Numpy can do the calculation for all particles at once.

---

Since all the computations are done element-wise, it should be easy to re-write the expression in Cython. This will avoid all those very large temporary array that get created when you do oldYD-YD and such.

Another possibility is numexpr.

---

I would definitely go for numexpr. I'm not sure numexpr can handle indices, but I bet that the following (or something similar) would work:

import numexpr as ne

yold = self.oldYD[ind]
y = self.YD[ind]
xold = self.oldXD[ind]
x = self.XD[ind]
YD_zero = ne.evaluate("yold - ((yold - y) * xold)/(xold - x)")