def angle(P, Q):
    """
    P and Q are arrays of either dimension 1 or 2.

    If P is a matrix, then angles are computed wrt the rows of P. If Q
    is a matrix, then angles are computed wrt the columns of Q.

    """
    P = P.astype(np.float)
    Q = Q.astype(np.float)

    if P.ndim < 2:
        P = P.reshape((1,P.size))
    if Q.ndim < 2:
        Q = Q.reshape((Q.size,1))

    # Normalize P row-wise and Q column-wise
    P /= np.sqrt((P*P).sum(1)[:,np.newaxis])
    Q /= np.sqrt((Q*Q).sum(0)[np.newaxis,:])

    cos_PQ = np.dot(P,Q)
    
    # Rounding errors may result in values outside of the domain of
    # inverse cosine.
    cos_PQ[(cos_PQ > 1)] = 1.
    cos_PQ[(cos_PQ < -1)] = -1.

    out = np.arccos(cos_PQ)

    out = np.squeeze(out)
    if out.ndim == 0:
        return out[()]
    return out
dataset/ETHPy150Open inpho/vsm/vsm/spatial.py/angle/original