chronaxis.core.rate_numpy

源码基线：8904c2f744951857659abd40aaa3535d57faf53c

1"""用原dtype列、局部时间和中心化NumPy归约计算差分与等权拟合。23仅运行当前窗口，整数先在无符号幅度列计算精确差，再按IEEE表示检查提升。4拟合使用整列点积，最终只对两个标量归约结果还原倍率、秒与输出质量。5"""67import numpy as np8from .numeric import ScalarValue9from .numeric_numpy import DTYPES, _output, to_float_numpy10from .types import QualityFlag as Q, SignalDType as D111213def _relative_column(samples):14    """将已知数值打包为相对首样本的浮点列，给出纵轴缩放的二进制指数。1516    参数: samples为同型有效ScalarValue序列。返回: (相对列,指数,质量)。17    异常: 内存分配或内部结构错误传播。整数差转float不精确时返回精度质量。18    """19    dtype = samples[0].dtype20    raw = np.fromiter((v.value for v in samples),dtype=DTYPES[dtype],count=len(samples))21    if dtype is D.FLOAT64:22        # 先按2的幂缩放，异号相减不会产生无穷；单位倍率最后一次恢复。23        _, power = np.frexp(np.max(np.abs(raw)))24        scaled = np.ldexp(raw,-int(power))25        return scaled - scaled[0], int(power), Q.NONE26    # 对有符号整数翻转符号位后，uint64比较顺序恰好等于原带符号大小。27    ordered = raw.view(np.uint64)28    if dtype is D.INT64:29        ordered = ordered ^ np.uint64(1 << 63)30    base = ordered[0]31    magnitudes = np.maximum(ordered,base) - np.minimum(ordered,base)32    converted = to_float_numpy(_output(D.UINT64,magnitudes,np.zeros(len(samples),dtype=np.uint16)))33    quality = Q(int(np.bitwise_or.reduce(converted.quality,initial=0)))34    return np.where(ordered < base,-converted.values,converted.values), 0, quality353637def _restore(numerator, denominator, exponent, elapsed_ns):38    """还原局部纵轴数量级与所选纳秒跨度，不在中间乘法中溢出。3940    参数: 前两项为浮点归约的分子分母；exponent为纵轴2的幂；elapsed_ns为正数。41    返回: FLOAT64 ScalarValue。异常: 不可表示的最终数值以质量标识。42    """43    a,b = float(numerator).as_integer_ratio()44    c,d = float(denominator).as_integer_ratio()45    top = a * d * 1_000_000_000 * elapsed_ns.denominator46    bottom = b * c * elapsed_ns.numerator47    if exponent > 0:48        top <<= exponent49    else:50        bottom <<= -exponent51    try:52        value = top / bottom53    except OverflowError:54        return ScalarValue(None,D.FLOAT64,False,Q.NONFINITE)55    if value == 0 and top:56        return ScalarValue(None,D.FLOAT64,False,Q.PRECISION_LOSS)57    return ScalarValue(value,D.FLOAT64)585960def difference_numpy(left, right, elapsed_ns):61    """以原dtype两点列求差，恢复为所选周期或实际秒数上的导数。6263    参数: 同型有效的left/right及正elapsed_ns。返回: FLOAT64 ScalarValue。64    异常: 输入结构和分配错误传播。65    """66    values, exponent, quality = _relative_column((left,right))67    if quality:68        return ScalarValue(None,D.FLOAT64,False,quality)69    return _restore(values[-1],1.,exponent,elapsed_ns)707172def fit_numpy(positions, samples):73    """将实际或周期横轴先局部化，用整列中心积求等权回归斜率。7475    参数: positions为严格递增精确纳秒序列，samples为同行序有效样本。76    返回: FLOAT64 ScalarValue。异常: 结构和空间错误传播。77    """78    first, span = positions[0], positions[-1]-positions[0]79    if (all(type(t) is int for t in positions) and -(1 << 63) <= first80            and positions[-1] < (1 << 63) and span < (1 << 63)):81        # 普通实际时间在int64域先相减，靠近纪元上限的1ns仍然可区分。82        axis = np.asarray(positions,dtype=np.int64)83        x = (axis - axis[0]).astype(np.float64) / span84    else:85        # 介于两个驱动记录间的位置具有有理数表示，先在精确域移除起点；86        # 只打包相对坐标，后续中心化和每个样本的乘法均执行整列运算。87        x = np.fromiter((float((t-first)/span) for t in positions),dtype=np.float64,count=len(positions))88    y, exponent, quality = _relative_column(samples)89    if quality:90        return ScalarValue(None,D.FLOAT64,False,quality)91    x -= np.mean(x)92    y -= np.mean(y)93    return _restore(np.dot(x,y),np.dot(x,x),exponent,span)94