cumprod() with quaternions crashes
#225
Comments
|
Thanks for the bug report. This is pretty weird. I can consistently get the segfault with a variety of combinations of numpy and/or quaternion — but only with python 3.9 or later, whether I fix the numpy and quaternion versions across python versions, or update to the latest in each case. Also, it consistently works with 501 or fewer quaternions, and fails with 502 or more, on both macOS ARM, and Linux x86_64. Also worth noting that this behavior seems to happen with It may take some time for me to figure out this bug, and I imagine I'll need some insight from numpy devs. Meanwhile, a couple remarks. First, it's not entirely obvious to me that |
|
Thanks for the analysis and the nice library. I'd be surprised if this is an issue with A manual Yes, this is an integration. No, I'm not doing first order. It's an implementation of a high order geometric integrator based the Magnus expansion. That systematically beats the RK of |
To be clear, I'm not blaming numpy for being buggy, but there is no contract; there's just a minimally documented array finterface that's been going through huge changes these past few years.
Have you tried it? For small arrays, they're almost identical; for ~1,000 elements and up it's one order of magnitude. If that really matters, you might want to take a look at numba and quaternionic, which avoid many of the performance penalties inherent in numpy. (Of course, you'd have to roll your own there because
Ooh! I'd be interested to see it when you're done. |
Argh.
True. The manual impl vs
Will certainly publish once done. Basically, the core part is this. def magnus6(w):
"""exponential map for 6th order magnus expansion for integration of `w`
`q'(t) = q(t)*w(t)`
`w` (3, n): pure quaternion samples (given as vector so that `np.cross` works)
of w at the 3 nodes `_magnus6_nodes` in the n sample intervals `i`
i.e. supply `w[j, i] = h*w(h*(i + magnus6_nodes[j]))`
returns the optimal pure quaternion g[:n] such that
q[i + 1] = q[i]*exp(g[i])
Note the right multiplication for naturally ordered arrays where increasing index
is increasing time
Implementation of:
S. Blanes, F. Casas, and J. Ros. High order optimized geometric integrators
for linear differential equations. BIT, 42:262–284, 2002.
"""
a1 = w[1]
a2 = np.sqrt(5/3)*(w[0] - w[2])
a3 = 5/3*(w[0] - 2*w[1] + w[2])
s1 = np.cross(a1, a2)
s2 = 1/15*np.cross(a1, 2*a3 + s1)
g = a1 + 1/6*a3 + 1/60*np.cross(-10*a1 - a3 + s1, a2 - s2)
return g
# nodes
_magnus6_nodes = 1/2 + np.array([-1, 0, 1])*np.sqrt(3/5)/2
def vec2q(p):
"""Convert R3 vector to pure quaternion"""
shape = p.shape[:-1]
q = np.empty(shape + (4,), np.float64)
q[..., 0] = 0.
q[..., 1:] = p
return q.view(np.quaternion).reshape(shape)
# integrating
w0 = 1. # splitting
f0 = 1. # drive frequency
g0 = .8 # drive strength
# standard qubit drive (angular velocity)
def w(t):
return np.array([
g0*np.cos(f0*t),
g0*np.sin(f0*t),
w0/2*np.ones_like(t),
]).T
# t: time steps
# h: step size
h = 1e-2
wp = np.sqrt((w0 + f0)**2 + 4*g0**2)
t1 = 1000*2*np.pi/wp
t = np.arange(0, t1, h)
np.prod(np.exp(vec2q(magnus6(h*w((t[:, None] + h*(_magnus6_nodes - 1)).ravel()).reshape((t.shape[0], 3, 3)).swapaxes(0, 1)))))
# vs
ode = integrate.ode(lambda t, y: (quaternion.quaternion(*y)*quaternion.quaternion(*w(t))).components)
ode.set_initial_value([1, 0, 0, 0.], 0.)
ode.set_integrator("dop853", ...)On the toy problem of Blanes et al. it's about two orders of magnitude faster to get to machine precision than dop853. I don't need adaptive step size in my use case and I haven't done any optimization or tweaking. |
|
I've encountered this issue as well, running with Python 3.11.7 and numpy 1.26.3. Have you been able to make any progress locating the issue? |
Describe the bug
Calling
numpy.cumprod()on an array of quaternions segfaults.To Reproduce
Result:
Expected behavior
Runs without crashing
Environment (please complete the following information):
Additional context
gdb backtrace
The text was updated successfully, but these errors were encountered: