010-通用函数ufunc

010-通用函数ufunc

学习目标

通过本章节的学习，你将掌握：

• 通用函数(ufunc)的概念和特性
• ufunc的方法和属性详解
• 自定义ufunc的创建和使用
• ufunc的性能优化技巧
• ufunc在向量化计算中的高级应用
• ufunc与广播机制的结合使用

1. ufunc基础概念

1.1 什么是ufunc

通用函数(Universal Function, ufunc)是NumPy中对数组进行逐元素操作的函数。它们支持数组广播、类型转换和其他标准数组操作。

import numpy as np
import time

# ufunc基础示例
print("=== ufunc基础示例 ===")

# 创建测试数组
arr = np.array([1, 2, 3, 4, 5])
print(f"输入数组: {arr}")

# 一元ufunc示例
print(f"\n一元ufunc:")
print(f"np.sqrt(arr): {np.sqrt(arr)}")
print(f"np.exp(arr): {np.exp(arr)}")
print(f"np.sin(arr): {np.sin(arr)}")
print(f"np.log(arr): {np.log(arr)}")

# 二元ufunc示例
arr2 = np.array([2, 3, 4, 5, 6])
print(f"\n二元ufunc:")
print(f"第二个数组: {arr2}")
print(f"np.add(arr, arr2): {np.add(arr, arr2)}")
print(f"np.multiply(arr, arr2): {np.multiply(arr, arr2)}")
print(f"np.power(arr, arr2): {np.power(arr, arr2)}")
print(f"np.maximum(arr, arr2): {np.maximum(arr, arr2)}")

# ufunc vs 普通函数的性能比较
print(f"\n=== 性能比较 ===")

# 创建大数组
large_arr = np.random.random(1000000)

# 使用Python循环
def python_sqrt(arr):
    return [x**0.5 for x in arr]

# 使用NumPy ufunc
def numpy_sqrt(arr):
    return np.sqrt(arr)

# 性能测试
test_arr = large_arr[:10000]  # 使用较小数组测试Python循环

start_time = time.time()
result1 = python_sqrt(test_arr)
python_time = time.time() - start_time

start_time = time.time()
result2 = numpy_sqrt(large_arr)
numpy_time = time.time() - start_time

print(f"Python循环 (10K元素): {python_time:.6f}秒")
print(f"NumPy ufunc (1M元素): {numpy_time:.6f}秒")
print(f"NumPy效率提升: {python_time/(numpy_time*100):.1f}倍 (考虑数组大小差异)")

1.2 ufunc的特性

# ufunc的特性
print("=== ufunc的特性 ===")

# 1. 广播支持
print("1. 广播支持:")
matrix = np.array([[1, 2, 3], [4, 5, 6]])
vector = np.array([10, 20, 30])
scalar = 100

print(f"矩阵:\n{matrix}")
print(f"向量: {vector}")
print(f"标量: {scalar}")

print(f"\n矩阵 + 向量 (广播):\n{np.add(matrix, vector)}")
print(f"矩阵 + 标量 (广播):\n{np.add(matrix, scalar)}")

# 2. 类型提升
print(f"\n2. 类型提升:")
int_arr = np.array([1, 2, 3], dtype=np.int32)
float_arr = np.array([1.5, 2.5, 3.5], dtype=np.float64)

print(f"整数数组类型: {int_arr.dtype}")
print(f"浮点数组类型: {float_arr.dtype}")

result = np.add(int_arr, float_arr)
print(f"相加结果类型: {result.dtype}")
print(f"相加结果: {result}")

# 3. 输出数组控制
print(f"\n3. 输出数组控制:")
a = np.array([1, 2, 3, 4])
b = np.array([2, 3, 4, 5])
output = np.zeros(4, dtype=np.float64)

print(f"输入a: {a}")
print(f"输入b: {b}")
print(f"输出数组(初始): {output}")

# 使用out参数
np.add(a, b, out=output)
print(f"输出数组(计算后): {output}")
print(f"输出数组类型: {output.dtype}")

# 4. 条件计算
print(f"\n4. 条件计算:")
x = np.array([-2, -1, 0, 1, 2])
y = np.array([1, 1, 1, 1, 1])
condition = x > 0

print(f"数组x: {x}")
print(f"数组y: {y}")
print(f"条件(x > 0): {condition}")

# 只在满足条件的地方计算
result = np.zeros_like(x)
np.add(x, y, out=result, where=condition)
print(f"条件计算结果: {result}")

2. ufunc的方法

2.1 reduce方法

# ufunc的reduce方法
print("=== ufunc的reduce方法 ===")

# 基本reduce操作
arr = np.array([1, 2, 3, 4, 5])
print(f"输入数组: {arr}")

# 不同ufunc的reduce操作
print(f"\n基本reduce操作:")
print(f"np.add.reduce(arr): {np.add.reduce(arr)}")
print(f"np.multiply.reduce(arr): {np.multiply.reduce(arr)}")
print(f"np.maximum.reduce(arr): {np.maximum.reduce(arr)}")
print(f"np.minimum.reduce(arr): {np.minimum.reduce(arr)}")

# 等价的内置函数
print(f"\n等价的内置函数:")
print(f"np.sum(arr): {np.sum(arr)}")
print(f"np.prod(arr): {np.prod(arr)}")
print(f"np.max(arr): {np.max(arr)}")
print(f"np.min(arr): {np.min(arr)}")

# 多维数组的reduce
print(f"\n=== 多维数组的reduce ===")
matrix = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(f"矩阵:\n{matrix}")

# 沿不同轴的reduce
print(f"\n沿不同轴的reduce:")
print(f"沿轴0 (列): {np.add.reduce(matrix, axis=0)}")
print(f"沿轴1 (行): {np.add.reduce(matrix, axis=1)}")
print(f"全部元素: {np.add.reduce(matrix, axis=None)}")

# 保持维度
print(f"\n保持维度:")
print(f"keepdims=True (轴0): {np.add.reduce(matrix, axis=0, keepdims=True)}")
print(f"keepdims=True (轴1): {np.add.reduce(matrix, axis=1, keepdims=True)}")

# 自定义reduce操作
print(f"\n=== 自定义reduce操作 ===")

# 计算几何平均数
def geometric_mean_reduce(arr):
    """使用reduce计算几何平均数"""
    product = np.multiply.reduce(arr)
    return np.power(product, 1.0/len(arr))

test_arr = np.array([2, 4, 8, 16])
print(f"测试数组: {test_arr}")
print(f"几何平均数: {geometric_mean_reduce(test_arr)}")
print(f"验证: {np.power(np.prod(test_arr), 1.0/len(test_arr))}")

# 逻辑运算的reduce
print(f"\n=== 逻辑运算的reduce ===")
bool_arr = np.array([True, False, True, True])
print(f"布尔数组: {bool_arr}")
print(f"logical_and.reduce: {np.logical_and.reduce(bool_arr)}")
print(f"logical_or.reduce: {np.logical_or.reduce(bool_arr)}")

# 等价操作
print(f"np.all(): {np.all(bool_arr)}")
print(f"np.any(): {np.any(bool_arr)}")

2.2 accumulate方法

# ufunc的accumulate方法
print("=== ufunc的accumulate方法 ===")

# 基本accumulate操作
arr = np.array([1, 2, 3, 4, 5])
print(f"输入数组: {arr}")

# 不同ufunc的accumulate操作
print(f"\n基本accumulate操作:")
print(f"np.add.accumulate(arr): {np.add.accumulate(arr)}")
print(f"np.multiply.accumulate(arr): {np.multiply.accumulate(arr)}")
print(f"np.maximum.accumulate(arr): {np.maximum.accumulate(arr)}")
print(f"np.minimum.accumulate(arr): {np.minimum.accumulate(arr)}")

# 等价的内置函数
print(f"\n等价的内置函数:")
print(f"np.cumsum(arr): {np.cumsum(arr)}")
print(f"np.cumprod(arr): {np.cumprod(arr)}")

# 多维数组的accumulate
print(f"\n=== 多维数组的accumulate ===")
matrix = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(f"矩阵:\n{matrix}")

# 沿不同轴的accumulate
print(f"\n沿轴0的累积和:\n{np.add.accumulate(matrix, axis=0)}")
print(f"\n沿轴1的累积和:\n{np.add.accumulate(matrix, axis=1)}")

# 实际应用：股价分析
print(f"\n=== 实际应用：股价分析 ===")

# 模拟股价变化
np.random.seed(42)
daily_returns = np.random.normal(0.001, 0.02, 20)  # 日收益率
initial_price = 100

print(f"日收益率: {daily_returns[:5]}...")
print(f"初始价格: {initial_price}")

# 计算累积收益
cumulative_returns = np.add.accumulate(daily_returns)
print(f"\n累积收益率: {cumulative_returns[:5]}...")

# 计算股价序列
prices = initial_price * (1 + cumulative_returns)
print(f"股价序列: {prices[:5]}...")

# 使用multiply.accumulate计算复利
compound_factors = 1 + daily_returns
compound_growth = np.multiply.accumulate(compound_factors)
prices_compound = initial_price * compound_growth

print(f"\n复利计算的股价: {prices_compound[:5]}...")
print(f"最终价格差异: {abs(prices[-1] - prices_compound[-1]):.6f}")

# 计算最大回撤
running_max = np.maximum.accumulate(prices_compound)
drawdowns = (prices_compound - running_max) / running_max
max_drawdown = np.minimum.accumulate(drawdowns)[-1]

print(f"\n风险指标:")
print(f"最高价格: {np.max(running_max):.2f}")
print(f"最大回撤: {max_drawdown*100:.2f}%")

2.3 reduceat方法

# ufunc的reduceat方法
print("=== ufunc的reduceat方法 ===")

# 基本reduceat操作
arr = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
indices = [0, 3, 6, 9]

print(f"输入数组: {arr}")
print(f"分割索引: {indices}")

# 使用reduceat进行分组求和
result = np.add.reduceat(arr, indices)
print(f"\n分组求和结果: {result}")
print(f"解释:")
print(f"  组1 [0:3): {arr[0:3]} -> 求和 = {np.sum(arr[0:3])}")
print(f"  组2 [3:6): {arr[3:6]} -> 求和 = {np.sum(arr[3:6])}")
print(f"  组3 [6:9): {arr[6:9]} -> 求和 = {np.sum(arr[6:9])}")
print(f"  组4 [9:10): {arr[9:10]} -> 求和 = {np.sum(arr[9:10])}")

# 不同ufunc的reduceat
print(f"\n不同ufunc的reduceat:")
print(f"multiply.reduceat: {np.multiply.reduceat(arr, indices)}")
print(f"maximum.reduceat: {np.maximum.reduceat(arr, indices)}")
print(f"minimum.reduceat: {np.minimum.reduceat(arr, indices)}")

# 多维数组的reduceat
print(f"\n=== 多维数组的reduceat ===")
matrix = np.arange(24).reshape(6, 4)
print(f"矩阵 (6x4):\n{matrix}")

# 沿轴0的reduceat
row_indices = [0, 2, 4]
print(f"\n行分组索引: {row_indices}")
row_result = np.add.reduceat(matrix, row_indices, axis=0)
print(f"行分组求和结果:\n{row_result}")

# 沿轴1的reduceat
col_indices = [0, 2]
print(f"\n列分组索引: {col_indices}")
col_result = np.add.reduceat(matrix, col_indices, axis=1)
print(f"列分组求和结果:\n{col_result}")

# 实际应用：时间序列数据聚合
print(f"\n=== 实际应用：时间序列聚合 ===")

# 模拟24小时的温度数据
np.random.seed(42)
hourly_temps = 20 + 10 * np.sin(np.linspace(0, 2*np.pi, 24)) + np.random.normal(0, 2, 24)
print(f"24小时温度数据: {hourly_temps[:6]}...")

# 按6小时分组计算平均温度
group_indices = [0, 6, 12, 18]
group_sums = np.add.reduceat(hourly_temps, group_indices)
group_counts = np.array([6, 6, 6, 6])  # 每组6小时
group_averages = group_sums / group_counts

print(f"\n6小时分组平均温度:")
time_periods = ['00-06', '06-12', '12-18', '18-24']
for period, avg_temp in zip(time_periods, group_averages):
    print(f"  {period}: {avg_temp:.2f}°C")

# 计算每组的最高和最低温度
group_max = np.maximum.reduceat(hourly_temps, group_indices)
group_min = np.minimum.reduceat(hourly_temps, group_indices)

print(f"\n温度范围:")
for period, min_t, max_t in zip(time_periods, group_min, group_max):
    print(f"  {period}: {min_t:.2f}°C ~ {max_t:.2f}°C")

2.4 outer方法

# ufunc的outer方法
print("=== ufunc的outer方法 ===")

# 基本outer操作
a = np.array([1, 2, 3])
b = np.array([10, 20, 30, 40])

print(f"数组a: {a}")
print(f"数组b: {b}")

# 不同ufunc的outer操作
print(f"\n加法outer:")
add_outer = np.add.outer(a, b)
print(f"np.add.outer(a, b):\n{add_outer}")

print(f"\n乘法outer:")
mul_outer = np.multiply.outer(a, b)
print(f"np.multiply.outer(a, b):\n{mul_outer}")

print(f"\n最大值outer:")
max_outer = np.maximum.outer(a, b)
print(f"np.maximum.outer(a, b):\n{max_outer}")

# 理解outer操作
print(f"\n=== 理解outer操作 ===")
print(f"outer操作相当于:")
print(f"对于每个a[i]和b[j]，计算ufunc(a[i], b[j])")
print(f"\n手动验证 add.outer:")
for i, ai in enumerate(a):
    for j, bj in enumerate(b):
        print(f"a[{i}] + b[{j}] = {ai} + {bj} = {ai + bj}")

# 高维数组的outer
print(f"\n=== 高维数组的outer ===")
matrix1 = np.array([[1, 2], [3, 4]])
vector = np.array([10, 20])

print(f"矩阵1:\n{matrix1}")
print(f"向量: {vector}")

# 注意：outer会展平输入数组
outer_result = np.add.outer(matrix1, vector)
print(f"\nouter结果形状: {outer_result.shape}")
print(f"outer结果:\n{outer_result}")

print(f"\n解释: matrix1被展平为 {matrix1.flatten()}")
print(f"然后与vector进行outer操作")

# 实际应用：距离矩阵计算
print(f"\n=== 实际应用：距离矩阵 ===")

# 一维坐标点
points_1d = np.array([1, 3, 7, 10])
print(f"一维坐标点: {points_1d}")

# 计算所有点对之间的距离
dist_matrix = np.abs(np.subtract.outer(points_1d, points_1d))
print(f"\n距离矩阵:\n{dist_matrix}")

# 二维坐标点的距离计算
print(f"\n=== 二维坐标距离 ===")
points_2d = np.array([[0, 0], [1, 1], [2, 0], [0, 2]])
print(f"二维坐标点:\n{points_2d}")

# 计算x坐标差和y坐标差
x_coords = points_2d[:, 0]
y_coords = points_2d[:, 1]

x_diff = np.subtract.outer(x_coords, x_coords)
y_diff = np.subtract.outer(y_coords, y_coords)

print(f"\nx坐标差矩阵:\n{x_diff}")
print(f"\ny坐标差矩阵:\n{y_diff}")

# 计算欧几里得距离
euclidean_dist = np.sqrt(x_diff**2 + y_diff**2)
print(f"\n欧几里得距离矩阵:\n{euclidean_dist}")

# 应用：相关性分析
print(f"\n=== 应用：相关性分析 ===")

# 模拟数据
np.random.seed(42)
data = np.random.randn(5, 100)  # 5个变量，100个观测
variable_names = ['A', 'B', 'C', 'D', 'E']

print(f"数据形状: {data.shape}")
print(f"变量名: {variable_names}")

# 计算均值
means = np.mean(data, axis=1)
print(f"\n各变量均值: {means}")

# 使用outer计算均值差异矩阵
mean_diff_matrix = np.subtract.outer(means, means)
print(f"\n均值差异矩阵:\n{mean_diff_matrix}")

# 找出差异最大的变量对
abs_diff = np.abs(mean_diff_matrix)
np.fill_diagonal(abs_diff, 0)  # 忽略对角线
max_diff_idx = np.unravel_index(np.argmax(abs_diff), abs_diff.shape)
print(f"\n差异最大的变量对: {variable_names[max_diff_idx[0]]} 和 {variable_names[max_diff_idx[1]]}")
print(f"差异值: {mean_diff_matrix[max_diff_idx]:.4f}")

3. 自定义ufunc

3.1 使用vectorize创建ufunc

# 使用vectorize创建ufunc
print("=== 使用vectorize创建ufunc ===")

# 定义标量函数
def sigmoid(x):
    """Sigmoid激活函数"""
    return 1 / (1 + np.exp(-x))

def relu(x):
    """ReLU激活函数"""
    return max(0, x)

def custom_distance(x1, y1, x2, y2):
    """计算两点间的欧几里得距离"""
    return np.sqrt((x2 - x1)**2 + (y2 - y1)**2)

# 向量化函数
print("创建向量化函数:")
vectorized_sigmoid = np.vectorize(sigmoid)
vectorized_relu = np.vectorize(relu)
vectorized_distance = np.vectorize(custom_distance)

print(f"sigmoid函数类型: {type(vectorized_sigmoid)}")
print(f"relu函数类型: {type(vectorized_relu)}")

# 测试向量化函数
test_input = np.array([-2, -1, 0, 1, 2])
print(f"\n测试输入: {test_input}")
print(f"sigmoid结果: {vectorized_sigmoid(test_input)}")
print(f"relu结果: {vectorized_relu(test_input)}")

# 测试多参数函数
print(f"\n=== 多参数函数测试 ===")
x1, y1 = np.array([0, 1, 2]), np.array([0, 1, 2])
x2, y2 = np.array([3, 4, 5]), np.array([4, 5, 6])

print(f"点1坐标: ({x1}, {y1})")
print(f"点2坐标: ({x2}, {y2})")

distances = vectorized_distance(x1, y1, x2, y2)
print(f"距离: {distances}")

# 指定输出类型
print(f"\n=== 指定输出类型 ===")

def int_division(a, b):
    """整数除法"""
    return a // b

# 不指定类型
vec_div1 = np.vectorize(int_division)
result1 = vec_div1(np.array([10, 20, 30]), np.array([3, 4, 5]))
print(f"默认类型结果: {result1}, 类型: {result1.dtype}")

# 指定输出类型
vec_div2 = np.vectorize(int_division, otypes=[np.int32])
result2 = vec_div2(np.array([10, 20, 30]), np.array([3, 4, 5]))
print(f"指定类型结果: {result2}, 类型: {result2.dtype}")

# 性能比较
print(f"\n=== 性能比较 ===")

large_array = np.random.randn(100000)

# 使用vectorize
start_time = time.time()
result_vec = vectorized_sigmoid(large_array)
vec_time = time.time() - start_time

# 使用NumPy内置函数
start_time = time.time()
result_numpy = 1 / (1 + np.exp(-large_array))
numpy_time = time.time() - start_time

print(f"vectorize时间: {vec_time:.6f}秒")
print(f"NumPy内置时间: {numpy_time:.6f}秒")
print(f"NumPy快: {vec_time/numpy_time:.1f}倍")
print(f"结果相同: {np.allclose(result_vec, result_numpy)}")

3.2 使用frompyfunc创建ufunc

# 使用frompyfunc创建ufunc
print("=== 使用frompyfunc创建ufunc ===")

# 定义Python函数
def add_with_threshold(a, b, threshold=10):
    """带阈值的加法：如果和超过阈值，返回阈值"""
    result = a + b
    return min(result, threshold)

def complex_operation(x, y):
    """复杂运算示例"""
    if x == 0:
        return 0
    return (x * y) / (x + y) if (x + y) != 0 else float('inf')

# 创建ufunc
print("创建ufunc:")
# frompyfunc(func, nin, nout)
# nin: 输入参数个数, nout: 输出参数个数
threshold_add_ufunc = np.frompyfunc(add_with_threshold, 3, 1)
complex_ufunc = np.frompyfunc(complex_operation, 2, 1)

print(f"threshold_add_ufunc类型: {type(threshold_add_ufunc)}")
print(f"complex_ufunc类型: {type(complex_ufunc)}")

# 测试ufunc
print(f"\n测试threshold_add_ufunc:")
a = np.array([1, 5, 8, 12])
b = np.array([2, 6, 7, 3])
threshold = np.array([10, 10, 10, 10])

result = threshold_add_ufunc(a, b, threshold)
print(f"输入a: {a}")
print(f"输入b: {b}")
print(f"阈值: {threshold}")
print(f"结果: {result}")
print(f"结果类型: {type(result[0])}")

# 注意：frompyfunc返回object数组
print(f"\n结果数组类型: {result.dtype}")

# 转换为数值类型
numeric_result = result.astype(float)
print(f"转换后类型: {numeric_result.dtype}")
print(f"转换后结果: {numeric_result}")

# 测试complex_ufunc
print(f"\n=== 测试complex_ufunc ===")
x = np.array([1, 2, 0, 4, 5])
y = np.array([2, 3, 5, -4, 5])

result_complex = complex_ufunc(x, y)
print(f"输入x: {x}")
print(f"输入y: {y}")
print(f"复杂运算结果: {result_complex}")

# 处理特殊值
print(f"\n处理特殊值:")
for i, (xi, yi, ri) in enumerate(zip(x, y, result_complex)):
    print(f"f({xi}, {yi}) = {ri}")

# ufunc方法测试
print(f"\n=== ufunc方法测试 ===")

# 简单的二元函数
def simple_max(a, b):
    return a if a > b else b

simple_max_ufunc = np.frompyfunc(simple_max, 2, 1)

test_arr = np.array([1, 5, 3, 8, 2])
print(f"测试数组: {test_arr}")

# reduce方法
max_result = simple_max_ufunc.reduce(test_arr)
print(f"reduce结果: {max_result}")
print(f"验证: {np.max(test_arr)}")

# accumulate方法
accum_result = simple_max_ufunc.accumulate(test_arr)
print(f"accumulate结果: {accum_result}")
print(f"验证: {np.maximum.accumulate(test_arr)}")

3.3 高性能自定义ufunc

# 高性能自定义ufunc
print("=== 高性能自定义ufunc ===")

# 使用numba加速（如果可用）
try:
    from numba import vectorize, float64
    
    @vectorize([float64(float64)])
    def fast_sigmoid(x):
        """使用numba加速的sigmoid函数"""
        return 1.0 / (1.0 + np.exp(-x))
    
    print("Numba可用，创建了加速版本")
    
    # 性能测试
    test_data = np.random.randn(1000000)
    
    # 测试numba版本
    start_time = time.time()
    result_numba = fast_sigmoid(test_data)
    numba_time = time.time() - start_time
    
    # 测试numpy版本
    start_time = time.time()
    result_numpy = 1.0 / (1.0 + np.exp(-test_data))
    numpy_time = time.time() - start_time
    
    print(f"\n性能比较 (100万元素):")
    print(f"Numba版本: {numba_time:.6f}秒")
    print(f"NumPy版本: {numpy_time:.6f}秒")
    print(f"加速比: {numpy_time/numba_time:.1f}倍")
    print(f"结果相同: {np.allclose(result_numba, result_numpy)}")
    
except ImportError:
    print("Numba不可用，跳过加速测试")

# 使用Cython的示例（概念性）
print(f"\n=== Cython优化示例（概念） ===")
print("""
# 在.pyx文件中定义Cython函数
import numpy as np
cimport numpy as cnp
from libc.math cimport exp

def cython_sigmoid(cnp.ndarray[double, ndim=1] x):
    cdef int i
    cdef int n = x.shape[0]
    cdef cnp.ndarray[double, ndim=1] result = np.empty(n)
    
    for i in range(n):
        result[i] = 1.0 / (1.0 + exp(-x[i]))
    
    return result

# 编译后可以获得接近C的性能
""")

# 内存高效的ufunc实现
print(f"\n=== 内存高效的实现 ===")

def memory_efficient_operation(input_array, chunk_size=10000):
    """分块处理大数组，节省内存"""
    n = len(input_array)
    result = np.empty_like(input_array)
    
    for i in range(0, n, chunk_size):
        end_idx = min(i + chunk_size, n)
        chunk = input_array[i:end_idx]
        
        # 在这里应用复杂的运算
        result[i:end_idx] = 1.0 / (1.0 + np.exp(-chunk))
    
    return result

# 测试内存效率
print("测试内存效率:")
large_data = np.random.randn(1000000)

start_time = time.time()
result_chunked = memory_efficient_operation(large_data)
chunked_time = time.time() - start_time

start_time = time.time()
result_direct = 1.0 / (1.0 + np.exp(-large_data))
direct_time = time.time() - start_time

print(f"分块处理时间: {chunked_time:.6f}秒")
print(f"直接处理时间: {direct_time:.6f}秒")
print(f"结果相同: {np.allclose(result_chunked, result_direct)}")
print(f"分块处理适用于内存受限的情况")

4. ufunc的高级应用

4.1 条件运算和掩码

# ufunc的条件运算和掩码
print("=== ufunc的条件运算和掩码 ===")

# 创建测试数据
data = np.array([-3, -1, 0, 2, 5, -2, 4, -4])
print(f"原始数据: {data}")

# 使用where参数进行条件运算
print(f"\n条件运算示例:")

# 只对正数进行平方运算
result1 = np.zeros_like(data, dtype=float)
condition1 = data > 0
np.square(data, out=result1, where=condition1)
print(f"只对正数平方: {result1}")
print(f"条件: {condition1}")

# 只对偶数进行开方运算
result2 = np.ones_like(data, dtype=float)  # 初始化为1
condition2 = (data % 2 == 0) & (data >= 0)
np.sqrt(np.abs(data), out=result2, where=condition2)
print(f"\n只对非负偶数开方: {result2}")
print(f"条件: {condition2}")

# 复杂条件运算
print(f"\n=== 复杂条件运算 ===")

# 创建2D数据
matrix = np.random.randn(4, 5)
print(f"随机矩阵:\n{matrix}")

# 条件：绝对值大于1的元素
condition = np.abs(matrix) > 1
print(f"\n条件(|x| > 1):\n{condition}")

# 对满足条件的元素应用sigmoid函数
result_matrix = np.zeros_like(matrix)
sigmoid_func = lambda x: 1 / (1 + np.exp(-x))
vectorized_sigmoid = np.vectorize(sigmoid_func)

# 使用where进行条件应用
np.copyto(result_matrix, vectorized_sigmoid(matrix), where=condition)
print(f"\n条件sigmoid结果:\n{result_matrix}")

# 使用np.where的替代方法
result_where = np.where(condition, vectorized_sigmoid(matrix), matrix)
print(f"\n使用np.where的结果:\n{result_where}")

# 多条件组合
print(f"\n=== 多条件组合 ===")

data_2d = np.random.randint(-10, 11, (3, 4))
print(f"整数矩阵:\n{data_2d}")

# 条件1：正数
cond_positive = data_2d > 0
# 条件2：偶数
cond_even = data_2d % 2 == 0
# 条件3：绝对值小于5
cond_small = np.abs(data_2d) < 5

# 组合条件：正偶数且绝对值小于5
combined_condition = cond_positive & cond_even & cond_small

print(f"\n组合条件结果:\n{combined_condition}")

# 对满足条件的元素乘以10
result_combined = data_2d.copy()
np.multiply(result_combined, 10, out=result_combined, where=combined_condition)
print(f"\n条件运算结果:\n{result_combined}")

4.2 类型转换和精度控制

# ufunc的类型转换和精度控制
print("=== ufunc的类型转换和精度控制 ===")

# 创建不同类型的数组
int_array = np.array([1, 2, 3, 4], dtype=np.int32)
float_array = np.array([1.5, 2.5, 3.5, 4.5], dtype=np.float64)

print(f"整数数组: {int_array}, 类型: {int_array.dtype}")
print(f"浮点数组: {float_array}, 类型: {float_array.dtype}")

# 自动类型提升
result_auto = np.add(int_array, float_array)
print(f"\n自动类型提升结果: {result_auto}, 类型: {result_auto.dtype}")

# 指定输出类型
print(f"\n=== 指定输出类型 ===")

# 强制输出为整数类型
result_int = np.add(int_array, float_array, dtype=np.int32)
print(f"强制整数输出: {result_int}, 类型: {result_int.dtype}")

# 强制输出为特定精度
result_float32 = np.add(int_array, float_array, dtype=np.float32)
print(f"float32输出: {result_float32}, 类型: {result_float32.dtype}")

# 使用casting参数控制类型转换
print(f"\n=== casting参数控制 ===")

try:
    # 'no' - 不允许任何类型转换
    result_no_cast = np.add(int_array, float_array, casting='no')
except TypeError as e:
    print(f"casting='no'失败: {e}")

# 'equiv' - 只允许等价类型转换
try:
    result_equiv = np.add(int_array, int_array, casting='equiv')
    print(f"casting='equiv'成功: {result_equiv}")
except TypeError as e:
    print(f"casting='equiv'失败: {e}")

# 'safe' - 只允许安全的类型转换
result_safe = np.add(int_array, float_array, casting='safe')
print(f"casting='safe': {result_safe}, 类型: {result_safe.dtype}")

# 'same_kind' - 允许同类型转换
result_same_kind = np.add(int_array, float_array, casting='same_kind')
print(f"casting='same_kind': {result_same_kind}, 类型: {result_same_kind.dtype}")

# 'unsafe' - 允许任何类型转换
result_unsafe = np.add(float_array, int_array, dtype=np.int8, casting='unsafe')
print(f"casting='unsafe': {result_unsafe}, 类型: {result_unsafe.dtype}")

# 精度损失示例
print(f"\n=== 精度损失示例 ===")

large_float = np.array([1e10, 1e15, 1e20], dtype=np.float64)
print(f"大浮点数: {large_float}")

# 转换为float32会损失精度
float32_result = large_float.astype(np.float32)
print(f"转换为float32: {float32_result}")
print(f"精度损失: {large_float - float32_result.astype(np.float64)}")

# 转换为整数会截断
int_result = large_float.astype(np.int64)
print(f"转换为int64: {int_result}")

# 溢出示例
print(f"\n=== 溢出示例 ===")

small_int = np.array([100, 200], dtype=np.int8)  # int8范围: -128到127
print(f"int8数组: {small_int}")

# 加法可能导致溢出
overflow_result = np.add(small_int, 50)
print(f"加50后: {overflow_result}")
print(f"是否溢出: {overflow_result[1] < 200 + 50}")

# 使用更大的数据类型避免溢出
safe_result = np.add(small_int.astype(np.int32), 50)
print(f"安全计算: {safe_result}")

4.3 性能优化技巧

# ufunc性能优化技巧
print("=== ufunc性能优化技巧 ===")

# 创建大数组进行性能测试
n = 1000000
array1 = np.random.random(n)
array2 = np.random.random(n)
array3 = np.random.random(n)

print(f"测试数组大小: {n:,} 元素")

# 技巧1：使用就地操作
print(f"\n=== 技巧1：就地操作 ===")

# 非就地操作
test_array = array1.copy()
start_time = time.time()
result1 = test_array + array2
result1 = result1 * array3
result1 = result1 - 1.0
non_inplace_time = time.time() - start_time

# 就地操作
test_array = array1.copy()
start_time = time.time()
test_array += array2
test_array *= array3
test_array -= 1.0
inplace_time = time.time() - start_time

print(f"非就地操作: {non_inplace_time:.6f}秒")
print(f"就地操作: {inplace_time:.6f}秒")
print(f"性能提升: {non_inplace_time/inplace_time:.1f}倍")

# 技巧2：使用out参数
print(f"\n=== 技巧2：使用out参数 ===")

# 不使用out参数
start_time = time.time()
for _ in range(10):
    temp1 = np.add(array1, array2)
    temp2 = np.multiply(temp1, array3)
    result = np.subtract(temp2, 1.0)
no_out_time = time.time() - start_time

# 使用out参数
temp_array = np.empty_like(array1)
result_array = np.empty_like(array1)
start_time = time.time()
for _ in range(10):
    np.add(array1, array2, out=temp_array)
    np.multiply(temp_array, array3, out=result_array)
    np.subtract(result_array, 1.0, out=result_array)
with_out_time = time.time() - start_time

print(f"不使用out: {no_out_time:.6f}秒")
print(f"使用out: {with_out_time:.6f}秒")
print(f"性能提升: {no_out_time/with_out_time:.1f}倍")

# 技巧3：选择合适的数据类型
print(f"\n=== 技巧3：数据类型优化 ===")

# float64 vs float32
array_f64 = np.random.random(n).astype(np.float64)
array_f32 = array_f64.astype(np.float32)

print(f"float64内存: {array_f64.nbytes / 1024**2:.2f} MB")
print(f"float32内存: {array_f32.nbytes / 1024**2:.2f} MB")

# 性能测试
start_time = time.time()
for _ in range(100):
    result_f64 = np.sqrt(array_f64)
f64_time = time.time() - start_time

start_time = time.time()
for _ in range(100):
    result_f32 = np.sqrt(array_f32)
f32_time = time.time() - start_time

print(f"\nfloat64计算: {f64_time:.6f}秒")
print(f"float32计算: {f32_time:.6f}秒")
print(f"float32快: {f64_time/f32_time:.1f}倍")

# 技巧4：避免不必要的数组创建
print(f"\n=== 技巧4：避免数组创建 ===")

# 低效：创建临时数组
start_time = time.time()
for _ in range(100):
    mask = array1 > 0.5
    filtered = array1[mask]
    result = np.sum(filtered)
inefficient_time = time.time() - start_time

# 高效：直接计算
start_time = time.time()
for _ in range(100):
    result = np.sum(array1, where=array1 > 0.5)
efficient_time = time.time() - start_time

print(f"创建临时数组: {inefficient_time:.6f}秒")
print(f"直接计算: {efficient_time:.6f}秒")
print(f"性能提升: {inefficient_time/efficient_time:.1f}倍")

# 技巧5：批量操作
print(f"\n=== 技巧5：批量操作 ===")

# 分别计算多个函数
start_time = time.time()
for _ in range(10):
    sin_result = np.sin(array1)
    cos_result = np.cos(array1)
    tan_result = np.tan(array1)
separate_time = time.time() - start_time

# 一次性计算（如果可能）
start_time = time.time()
for _ in range(10):
    # 虽然这个例子不能真正批量化，但展示了概念
    sin_result = np.sin(array1)
    cos_result = np.cos(array1)
    # 利用三角恒等式：tan = sin/cos
    tan_result = sin_result / cos_result
batch_time = time.time() - start_time

print(f"分别计算: {separate_time:.6f}秒")
print(f"利用关系: {batch_time:.6f}秒")
print(f"性能差异: {separate_time/batch_time:.1f}倍")

5. 总结

本章节深入介绍了NumPy的通用函数(ufunc)，包括：

核心概念

1. ufunc定义：对数组进行逐元素操作的向量化函数
2. ufunc特性：广播支持、类型提升、输出控制、条件计算
3. 性能优势：比Python循环快数十倍到数百倍
4. 内存效率：支持就地操作和输出数组控制

ufunc方法

1. reduce：沿指定轴聚合数组元素
2. accumulate：计算累积结果
3. reduceat：分组聚合操作
4. outer：计算两个数组的外积运算

自定义ufunc

1. vectorize：简单易用，适合原型开发
2. frompyfunc：更底层的控制，返回object数组
3. 性能考虑：自定义ufunc通常比内置函数慢
4. 优化策略：使用numba、Cython等工具加速

高级应用

1. 条件运算：使用where参数进行选择性计算
2. 类型控制：精确控制输入输出类型和转换规则
3. 性能优化：就地操作、out参数、数据类型选择
4. 内存管理：避免不必要的临时数组创建

最佳实践

1. 优先使用NumPy内置ufunc
2. 合理使用就地操作和out参数
3. 选择合适的数据类型平衡精度和性能
4. 利用条件运算避免不必要的计算
5. 注意类型转换和精度损失

性能优化要点

1. 使用就地操作减少内存分配
2. 利用out参数重用数组
3. 选择合适的数据类型
4. 避免创建临时数组
5. 批量处理相关操作

掌握ufunc的使用将大大提高你的NumPy编程效率和代码性能。

下一步学习

在下一章节 011-数组聚合操作 中，我们将深入学习：

• 各种聚合函数的使用和原理
• 多维数组的轴向聚合
• 分组聚合和条件聚合
• 聚合操作的性能优化
• 聚合在数据分析中的应用

这将帮助你掌握NumPy中数据汇总和统计分析的核心技能。