Different return value for multiple calls to function-CodePudding

I am encountering some strange behavior in my Python (v. 3.8.8) code: Multiple calls of the same function (with exactly the same input) returns different values. The code below reproduces the issue.

Expected behavior: Each call of the function with exactly the same input parameters should return the same output. In other words, the value of zmin should be constant for all iterations.

Observed behavior: Each call to the function gives a (slightly) different output, even though the input parameters are the same: zmin in the example is not constant. The difference is quite large, for instance the code below gives: zmin = array([1.94777151, 1.95983567, 1.97510284, 1.99350764, 2.01497125, 2.03940178, 2.0666946 , 2.09673286, 2.12938792, 2.16451995]).

Why does this happen?

The issue persists regardless of whether the helper function quarternion_mul is vectorized using @guvectorize.

import numpy as np
from numba import guvectorize
import matplotlib.pyplot as plt

# Update MC cycle parameters
def change_to(state,position=None,rotation=None):
    """Update the current state to the state described by the center of mass
    position and the applied rotation.
    """
    # If provided, apply rotation to all points:
    point = np.zeros(4,dtype=float)
    if rotation is not None:
        rotation_inverse = np.array([rotation[0], -rotation[1],
                                      -rotation[2], -rotation[3]])
        for i in range(1,state.shape[0]):
            # Translate center of mass to origin
            coords = state[i,:3] - state[0,:3]
            point[0] = 0.0
            point[1:] = coords
            # Apply rotation
            rp = quaternion_mul(rotation,point)
            rpr_inv = quaternion_mul(rp,rotation_inverse)
            state[i,:3] = rpr_inv[1:]   state[0,:3]
      
    # Move points to target location
    delta = np.zeros(len(state[0,:]),dtype=float)
    delta[:3] = position - state[0,:3]
    return state   delta

# Function to calculate Hamilton product
@guvectorize(['(float64[:], float64[:], float64[:])'], '(n),(n)->(n)', nopython=True, target="parallel")
def quaternion_mul(a,b,c):
    # Hamilton quaternion product
    c[0] = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3]
    c[1] = a[0]*b[1]   a[1]*b[0]   a[2]*b[3] - a[3]*b[2]
    c[2] = a[0]*b[2] - a[1]*b[3]   a[2]*b[0]   a[3]*b[1]
    c[3] = a[0]*b[3]   a[1]*b[2] - a[2]*b[1]   a[3]*b[0]

# def quaternion_mul(a,b):
#     # Hamilton quaternion product
#     c = np.zeros(4)
#     c[0] = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3]
#     c[1] = a[0]*b[1]   a[1]*b[0]   a[2]*b[3] - a[3]*b[2]
#     c[2] = a[0]*b[2] - a[1]*b[3]   a[2]*b[0]   a[3]*b[1]
#     c[3] = a[0]*b[3]   a[1]*b[2] - a[2]*b[1]   a[3]*b[0]
#     return c


start_state = np.array([[-0.053442, 2.05268,    1.92046,    0,  0,  0,  0],
[1.02709,   1.6784, 1.03511,    0.6365, 11310,  10156,  1],
[0.956003,  2.04269,    0.949338,   0.5853, 8225,   11248,  0],
[0.615004,  2.05092,    0.784356,   0.6365, 12229,  14318,  0],
[0.488786,  2.28627,    0.520153,   0.4509, 8414,   9459,   0],
[0.174876,  2.47819,    0.607007,   0.6563, 13127,  11105,  1],
[-0.0532915,    2.29712,    0.359422,   0.5477, 9610,   7634,   0],
[0.107118,  1.95669,    0.407201,   0.5927, 10666,  8236,   -1],
[0.0313548, 1.98191,    0.781065,   0.6187, 7632,   7583,   0]],dtype=float)

position = np.array([0,0,3.5])
rotation = np.array([0.9994646047913096,    0.0146483594354613, 
                     -0.0194742908484703,   0.0218330330268212])

iters = 10
zmin = np.zeros(iters)
for i in range(iters):
    end = change_to(start_state,
                    position=position,
                    rotation=rotation)
    zmin[i] = np.min(end[1:,2])

CodePudding user response：

You're modifying start_state.