Numpy Quaternion

import quaternion
import numpy as np

from quaternion import (
    as_rotation_matrix, from_rotation_matrix,
    as_euler_angles, from_euler_angles, 
    slerp, squad, integrate_angular_velocity
)

import quaternion
import numpy as np

# Create quaternions
q1 = quaternion.quaternion(1, 0, 0, 0)  # Identity quaternion
q2 = quaternion.quaternion(0, 1, 0, 0)  # Pure quaternion (i)

# Use predefined constants
identity = quaternion.one
x_rotation = quaternion.x
y_rotation = quaternion.y
z_rotation = quaternion.z

# Create quaternion arrays
q_array = np.array([q1, q2, quaternion.y, quaternion.z])

# Basic arithmetic
result = q1 * q2  # Quaternion multiplication
conjugate = q1.conjugate()  # Quaternion conjugate
magnitude = abs(q1)  # Quaternion norm

# Convert to rotation matrix
rotation_matrix = quaternion.as_rotation_matrix(q1)

# Create from Euler angles (radians)
alpha, beta, gamma = 0.1, 0.2, 0.3
q_from_euler = quaternion.from_euler_angles(alpha, beta, gamma)

# Interpolate between quaternions
t_in = np.linspace(0, 1, 5)
t_out = np.linspace(0, 1, 10)
R_in = np.array([quaternion.one, quaternion.x, quaternion.y])
R_out = quaternion.slerp(R_in, t_in, t_out)

# Core quaternion type and constants
quaternion(w, x, y, z)  # Constructor
quaternion.zero         # quaternion(0, 0, 0, 0)
quaternion.one          # quaternion(1, 0, 0, 0)  
quaternion.x            # quaternion(0, 1, 0, 0)
quaternion.y            # quaternion(0, 0, 1, 0)
quaternion.z            # quaternion(0, 0, 0, 1)

# Array conversion functions
def as_float_array(a): ...          # -> ndarray with shape (..., 4)
def as_quat_array(a): ...           # -> quaternion array
def as_vector_part(q): ...          # -> ndarray with shape (..., 3)
def from_vector_part(v, vector_axis=-1): ...  # -> quaternion array
def as_spinor_array(a): ...         # -> complex array (..., 2)

def as_rotation_matrix(q): ...      # -> ndarray (..., 3, 3)
def from_rotation_matrix(rot, nonorthogonal=True): ...  # -> quaternion array
def as_rotation_vector(q): ...      # -> ndarray (..., 3)
def from_rotation_vector(rot): ...  # -> quaternion array
def as_euler_angles(q): ...         # -> ndarray (..., 3)
def from_euler_angles(alpha_beta_gamma, beta=None, gamma=None): ...  # -> quaternion array
def as_spherical_coords(q): ...     # -> ndarray (..., 2)  
def from_spherical_coords(theta_phi, phi=None): ...  # -> quaternion array
def rotate_vectors(R, v, axis=-1): ...  # -> ndarray

def slerp(R1, R2, t1, t2, t_out): ...  # -> quaternion array
def squad(R_in, t_in, t_out, unflip_input_rotors=False): ...  # -> quaternion array
def integrate_angular_velocity(Omega, t0, t1, R0=None, tolerance=1e-12): ...  # -> quaternion array
def unflip_rotors(q, axis=-1, inplace=False): ...  # -> quaternion array
def angular_velocity(R, t): ...     # -> ndarray (..., 3)
def minimal_rotation(R, t, iterations=2): ...  # -> quaternion array

# Low-level evaluation functions
def slerp_evaluate(R1, R2, t):
    """
    Low-level spherical linear interpolation evaluation.
    
    Args:
        R1 (quaternion): Starting quaternion (must be normalized)
        R2 (quaternion): Ending quaternion (must be normalized)  
        t (float or array): Interpolation parameter(s) in [0, 1]
        
    Returns:
        quaternion or quaternion array: Interpolated quaternion(s)
        
    Notes:
        Direct evaluation function used internally by slerp().
        Assumes inputs are properly normalized unit quaternions.
    """

def squad_evaluate(R0, R1, R2, R3, t):
    """
    Low-level spherical quadrangle interpolation evaluation.
    
    Args:
        R0 (quaternion): Control point before start (normalized)
        R1 (quaternion): Start quaternion (normalized)
        R2 (quaternion): End quaternion (normalized)
        R3 (quaternion): Control point after end (normalized)
        t (float or array): Interpolation parameter(s) in [0, 1]
        
    Returns:
        quaternion or quaternion array: Interpolated quaternion(s)
        
    Notes:
        Direct evaluation function used internally by squad().
        Provides C1 continuous interpolation using spherical Bézier curves.
    """

def derivative(f, t):
    """
    Primary derivative function using best available method.
    
    Args:
        f (array_like): Function values to differentiate
        t (array_like): Time values
        
    Returns:
        ndarray: Derivative values with same shape as f
        
    Notes:
        Uses spline_derivative if SciPy available, otherwise finite difference.
        Automatically selects most accurate method for the given data.
    """

def definite_integral(f, t):
    """
    Primary definite integral function using best available method.
    
    Args:
        f (array_like): Function values to integrate
        t (array_like): Time values
        
    Returns:
        ndarray: Definite integral values
        
    Notes:
        Uses spline integration if SciPy available, otherwise trapezoidal rule.
        Automatically selects most accurate method.
    """

def indefinite_integral(f, t):
    """
    Primary indefinite integral function using best available method.
    
    Args:
        f (array_like): Function values to integrate
        t (array_like): Time values
        
    Returns:
        ndarray: Indefinite integral values with same shape as f
        
    Notes:
        Uses spline integration if SciPy available, otherwise trapezoidal rule.
        Automatically selects most accurate method.
    """

# Finite difference implementations
def fd_derivative(f, t):
    """
    Fourth-order finite difference derivative for non-uniform time steps.
    
    Args:
        f (array_like): Function values (1D, 2D, or 3D arrays)
        t (array_like): Time values (must be 1D)
        
    Returns:
        ndarray: Derivative values with same shape as f
        
    Notes:
        Uses fourth-order finite difference formula from Bowen & Smith.
        Requires at least 5 points for full accuracy.
    """

def fd_definite_integral(f, t):
    """
    Finite difference definite integral using trapezoidal rule.
    
    Args:
        f (array_like): Function values to integrate
        t (array_like): Time values
        
    Returns:
        ndarray: Definite integral values
    """

def fd_indefinite_integral(f, t):
    """
    Finite difference indefinite integral (cumulative integration).
    
    Args:
        f (array_like): Function values to integrate
        t (array_like): Time values
        
    Returns:
        ndarray: Indefinite integral values with same shape as f
    """

# Spline-based implementations (require SciPy)
def spline_derivative(f, t):
    """
    Spline-based derivative using cubic spline interpolation.
    
    Args:
        f (array_like): Function values
        t (array_like): Time values
        
    Returns:
        ndarray: Derivative values computed from spline
        
    Notes:
        More accurate than finite differences for smooth data.
        Requires SciPy. Falls back to finite differences if unavailable.
    """

def spline_definite_integral(f, t):
    """
    Spline-based definite integral using cubic spline interpolation.
    
    Args:
        f (array_like): Function values to integrate
        t (array_like): Time values
        
    Returns:
        ndarray: Definite integral computed from spline
        
    Notes:
        More accurate than trapezoidal rule for smooth data.
    """

def spline_indefinite_integral(f, t):
    """
    Spline-based indefinite integral using cubic spline interpolation.
    
    Args:
        f (array_like): Function values to integrate
        t (array_like): Time values
        
    Returns:
        ndarray: Indefinite integral values computed from spline
    """

def spline(f, t, t_out=None):
    """
    Create cubic spline interpolation of quaternion time series.
    
    Args:
        f (array_like): Function values to interpolate
        t (array_like): Time values corresponding to f
        t_out (array_like, optional): Output time values. If None, returns spline object.
        
    Returns:
        ndarray or spline object: Interpolated values or spline interpolator
        
    Notes:
        Available only when SciPy is installed.
        Returns InterpolatedUnivariateSpline object if t_out is None.
    """

def allclose(a, b, rtol=4*np.finfo(float).eps, atol=0.0, equal_nan=False, verbose=False):
    """
    Test if all quaternions in arrays are equal within tolerance.
    
    Args:
        a (array_like): First quaternion array
        b (array_like): Second quaternion array
        rtol (float): Relative tolerance parameter (default: 4*machine epsilon)  
        atol (float): Absolute tolerance parameter (default: 0.0)
        equal_nan (bool): Whether to consider NaN values as equal (default: False)
        verbose (bool): Print non-close values if False is returned (default: False)
        
    Returns:
        bool: True if all elements are close within tolerance, False otherwise
    """

def rotor_intrinsic_distance(a, b):
    """
    Intrinsic (geodesic) distance between unit quaternions.
    
    Args:
        a (quaternion array): First quaternion array (must be normalized)
        b (quaternion array): Second quaternion array (must be normalized)
        
    Returns:
        ndarray: Intrinsic distances with same broadcasting shape as inputs
        
    Notes:
        Computes geodesic distance on quaternion manifold (4D unit sphere).
        Distance = 2 * arccos(|⟨a,b⟩|) where ⟨⟩ is quaternion inner product.
    """

def rotor_chordal_distance(a, b):
    """
    Chordal (Euclidean) distance between unit quaternions.
    
    Args:
        a (quaternion array): First quaternion array (must be normalized)
        b (quaternion array): Second quaternion array (must be normalized)
        
    Returns:
        ndarray: Chordal distances with same broadcasting shape as inputs
        
    Notes:
        Standard Euclidean distance in 4D space: |a - b|.
        Simpler to compute than intrinsic distance.
    """

def rotation_intrinsic_distance(a, b):
    """
    Intrinsic distance between rotations represented by quaternions.
    
    Args:
        a (quaternion array): First quaternion array (must be normalized)
        b (quaternion array): Second quaternion array (must be normalized)
        
    Returns:
        ndarray: Rotation distances with same broadcasting shape as inputs
        
    Notes:
        Accounts for quaternion double-cover: ±q represent same rotation.
        Distance = min(rotor_intrinsic_distance(a,b), rotor_intrinsic_distance(a,-b)).
    """

def rotation_chordal_distance(a, b):
    """
    Chordal distance between rotations represented by quaternions.
    
    Args:
        a (quaternion array): First quaternion array (must be normalized)
        b (quaternion array): Second quaternion array (must be normalized)
        
    Returns:
        ndarray: Rotation chordal distances with same broadcasting shape as inputs
        
    Notes:
        Accounts for quaternion double-cover: ±q represent same rotation.
        Distance = min(|a - b|, |a + b|).
    """

def mean_rotor_in_chordal_metric(R, t=None):
    """
    Compute mean quaternion using chordal metric (least-squares sense).
    
    Args:
        R (quaternion array): Input quaternions (should be normalized)
        t (array_like, optional): Time weights for integral. If None, uses simple sum.
        
    Returns:
        quaternion: Normalized mean quaternion
        
    Raises:
        ValueError: If t is provided but len(t) < 4 or len(R) < 4
        
    Notes:
        Uses spline integration when t is provided (requires ≥4 points).
        Without t, computes simple arithmetic mean and normalizes.
    """

def optimal_alignment_in_chordal_metric(Ra, Rb, t=None):
    """
    Find rotation Rd such that Rd*Rb is closest to Ra in chordal metric.
    
    Args:
        Ra (quaternion array): Target quaternions
        Rb (quaternion array): Base quaternions to align
        t (array_like, optional): Time weights for integral
        
    Returns:
        quaternion: Optimal alignment rotation
        
    Notes:
        Minimizes ∫|Rd*Rb - Ra|² dt by computing mean of Ra/Rb.
        Equivalent to mean_rotor_in_chordal_metric(Ra / Rb, t).
    """

def optimal_alignment_in_Euclidean_metric(avec, bvec, t=None):
    """
    Solve Wahba's problem: find R such that R*b⃗*R̄ ≈ a⃗ in Euclidean metric.
    
    Args:
        avec (array_like): Target 3D vectors with shape (..., 3)
        bvec (array_like): Base 3D vectors with shape (..., 3)
        t (array_like, optional): Time weights for integral
        
    Returns:
        quaternion: Optimal rotation using Davenport's method
        
    Notes:
        Uses Davenport's eigenvector method for optimal vector alignment.
        Constructs matrix from vector cross-correlations and finds dominant eigenvector.
    """

# Core quaternion type (implemented in C extension)
class quaternion:
    """
    Quaternion number with w, x, y, z components representing q = w + xi + yj + zk.
    
    Attributes:
        w (float): Scalar component (real part)
        x (float): i component (first imaginary part)
        y (float): j component (second imaginary part)
        z (float): k component (third imaginary part)
        
    Notes:
        Components are stored as double-precision floating point.
        Supports all standard arithmetic operations and mathematical functions.
        Can be used in NumPy arrays with dtype=quaternion.
    """
    def __init__(self, w=0.0, x=0.0, y=0.0, z=0.0):
        """Initialize quaternion with given components."""
        
    def conjugate(self):
        """Return quaternion conjugate: q* = w - xi - yj - zk."""
        
    def norm(self):
        """Return quaternion norm (magnitude): |q| = sqrt(w² + x² + y² + z²)."""
        
    def normalized(self):
        """Return unit quaternion: q/|q|. Raises error if norm is zero."""
        
    def __add__(self, other):
        """Add quaternions or scalars component-wise."""
        
    def __mul__(self, other):
        """Multiply quaternions using Hamilton product or scale by scalar."""
        
    def __truediv__(self, other):
        """Divide by quaternion (multiply by inverse) or scalar."""
        
    def __pow__(self, other):
        """Raise quaternion to scalar power using q^p = exp(p * log(q))."""
        
    def exp(self):
        """Exponential: exp(q) = exp(w) * (cos(|v|) + sin(|v|)*v/|v|) where v=(x,y,z)."""
        
    def log(self):
        """Natural logarithm: log(q) = ln(|q|) + arccos(w/|q|) * v/|v| where v=(x,y,z)."""