【Numpy核心编程攻略：Python数据处理、分析详解与科学计算】1.23 数据工厂：高级初始化模式解析

1.23 数据工厂：高级初始化模式解析

目录

1.23.1 从数据库初始化的批量技巧
1.23.2 内存映射初始化大数组方案
1.23.3 随机初始化种子传播机制
1.23.4 初始化在数值计算中的陷阱
1.23.5 初始化参数的自动化测试

1.23.1 从数据库初始化的批量技巧

在实际应用中，我们经常需要从数据库中读取大量数据并将其初始化为 NumPy 数组。本节将介绍如何高效地将 SQL 数据转换为 NumPy 数组，并通过批量读取技巧来优化性能。

1.23.1.1 SQL到ndarray的高效转换

使用 pandas 库可以方便地将 SQL 查询结果转换为 NumPy 数组。通过批量读取和并行处理，可以显著提高数据转换的性能。

import pandas as pd
import numpy as np
import sqlite3
import time

# 连接到数据库
conn = sqlite3.connect('example.db')

# 创建一个示例表
conn.execute('''
CREATE TABLE IF NOT EXISTS data (
    id INTEGER PRIMARY KEY,
    value INTEGER
)
''')

# 插入示例数据
conn.executemany('INSERT INTO data (id, value) VALUES (?, ?)', [(i, i) for i in range(1000000)])
conn.commit()

# SQL 查询
query = 'SELECT value FROM data'

# 未优化的转换
start_time = time.time()
df = pd.read_sql(query, conn)
data_unoptimized = df['value'].values  # 转换为 NumPy 数组
end_time = time.time()
print("未优化的转换时间: ", end_time - start_time)

# 优化后的转换
start_time = time.time()
data_optimized = np.fromiter((row[0] for row in conn.execute(query)), dtype=int)
end_time = time.time()
print("优化后的转换时间: ", end_time - start_time)

# 关闭数据库连接
conn.close()

1.23.1.2 批量读取优化

批量读取是指一次从数据库中读取大量数据。通过批量读取，可以减少数据库连接和查询的次数，从而提高性能。

import pandas as pd
import numpy as np
import sqlite3
import time

# 连接到数据库
conn = sqlite3.connect('example.db')

# 创建一个示例表
conn.execute('''
CREATE TABLE IF NOT EXISTS data (
    id INTEGER PRIMARY KEY,
    value INTEGER
)
''')

# 插入示例数据
conn.executemany('INSERT INTO data (id, value) VALUES (?, ?)', [(i, i) for i in range(1000000)])
conn.commit()

# SQL 查询
query = 'SELECT value FROM data'

# 未优化的转换
start_time = time.time()
df = pd.read_sql(query, conn)
data_unoptimized = df['value'].values  # 转换为 NumPy 数组
end_time = time.time()
print("未优化的转换时间: ", end_time - start_time)

# 优化后的转换
def batch_read(conn, query, batch_size=10000):
    """
    批量读取优化

    :param conn: 数据库连接
    :param query: SQL 查询
    :param batch_size: 批量大小
    :return: 转换后的 NumPy 数组
    """
    cursor = conn.execute(query)
    data = []
    while True:
        batch = cursor.fetchmany(batch_size)  # 批量读取数据
        if not batch:
            break
        data.extend([row[0] for row in batch])
    return np.array(data, dtype=int)

start_time = time.time()
data_optimized = batch_read(conn, query)
end_time = time.time()
print("优化后的转换时间: ", end_time - start_time)

# 关闭数据库连接
conn.close()

SQL分块加载方案

import sqlite3
import numpy as np

def sql_to_ndarray(db_path, query, dtype=np.float32, chunk_size=10000):
    """ 高效数据库加载器 """
    conn = sqlite3.connect(db_path)
    cursor = conn.cursor()
    cursor.execute(query)
    
    # 预分配内存
    cursor.execute(f"SELECT COUNT(*) FROM ({query})")
    total = cursor.fetchone()[0]
    arr = np.empty((total,), dtype=dtype)
    
    # 分块加载
    index = 0
    while True:
        rows = cursor.fetchmany(chunk_size)
        if not rows: break
        chunk = np.array(rows, dtype=dtype).flatten()
        arr[index:index+len(chunk)] = chunk
        index += len(chunk)
    
    return arr

# 使用示例（加载百万级传感器数据）
data = sql_to_ndarray('sensor.db', 'SELECT value FROM measurements', chunk_size=5000)
print(f"加载数据量: {len(data):,}")  # 输出示例: 1,234,567

性能优化对比表

1.23.2 内存映射初始化大数组方案

内存映射是一种将文件映射到内存的技术，适用于处理大数组。本节将介绍如何使用内存映射来初始化大数组，并探讨其性能优势。

1.23.2.1 100GB数组的延迟初始化

通过内存映射，我们可以延迟加载大数组，从而减少内存开销。这对于处理100GB甚至更大的数组非常有用。

import numpy as np

# 创建一个 100GB 的文件
file_size = 100 * 1024 * 1024 * 1024  # 100GB
with open('large_array.npy', 'wb') as f:
    f.seek(file_size - 1)
    f.write(b'\0')

# 使用内存映射初始化数组
large_array = np.memmap('large_array.npy', dtype='float32', mode='r+', shape=(100000, 100000))

# 初始化数组
large_array[:] = np.random.rand(100000, 100000)  # 初始化数组
large_array.flush()  # 将数据写入文件

# 读取数组
loaded_array = np.memmap('large_array.npy', dtype='float32', mode='r', shape=(100000, 100000))

# 打印数组的一部分
print("数组的一部分: ", loaded_array[0:10, 0:10])

1.23.2.2 内存映射的性能优势

内存映射可以显著提高大数组的读写性能，因为它避免了将整个数组加载到内存中。

import numpy as np
import time

# 创建一个 100GB 的文件
file_size = 100 * 1024 * 1024 * 1024  # 100GB
with open('large_array.npy', 'wb') as f:
    f.seek(file_size - 1)
    f.write(b'\0')

# 使用内存映射初始化数组
large_array = np.memmap('large_array.npy', dtype='float32', mode='r+', shape=(100000, 100000))
large_array[:] = np.random.rand(100000, 100000)  # 初始化数组
large_array.flush()  # 将数据写入文件

# 读取数组
def read_array(array, index):
    """
    读取数组的一部分

    :param array: 输入的 NumPy 数组
    :param index: 读取的索引范围
    :return: 读取的数组部分
    """
    return array[index]

# 未优化的读取
start_time = time.time()
unoptimized_data = np.load('large_array.npy')  # 加载整个数组
unoptimized_result = read_array(unoptimized_data, (0, 10))  # 读取数组的一部分
end_time = time.time()
print("未优化的读取时间: ", end_time - start_time)

# 优化后的读取
start_time = time.time()
memmapped_array = np.memmap('large_array.npy', dtype='float32', mode='r', shape=(100000, 100000))
optimized_result = read_array(memmapped_array, (0, 10))  # 读取数组的一部分
end_time = time.time()
print("优化后的读取时间: ", end_time - start_time)

1.23.3 随机初始化种子传播机制

随机初始化是 NumPy 中常用的功能，但在并行计算中需要注意随机种子的管理和传播。本节将介绍并行随机数生成的安全性，并提供随机种子的管理和传播方法。

1.23.3.1 并行随机数生成的安全性

在并行计算中，多个进程或线程可能同时生成随机数。如果不正确管理随机种子，可能会导致随机数的重复生成，从而影响结果的正确性。

import numpy as np
import multiprocessing
import time

# 随机数生成函数
def generate_random(seed):
    """
    生成随机数

    :param seed: 随机种子
    :return: 生成的随机数组
    """
    np.random.seed(seed)
    return np.random.rand(100000)

# 未优化的并行生成
start_time = time.time()
with multiprocessing.Pool() as pool:
    results = pool.map(generate_random, [1] * 4)  # 使用相同的种子
end_time = time.time()
print("未优化的并行生成时间: ", end_time - start_time)

# 检查结果是否重复
print("结果是否重复: ", np.array_equal(results[0], results[1]))

# 优化后的并行生成
def generate_random_with_unique_seed(initial_seed, i):
    """
    生成随机数，使用唯一的种子

    :param initial_seed: 初始随机种子
    :param i: 进程编号
    :return: 生成的随机数组
    """
    np.random.seed(initial_seed + i)
    return np.random.rand(100000)

start_time = time.time()
with multiprocessing.Pool() as pool:
    results = pool.starmap(generate_random_with_unique_seed, [(1, i) for i in range(4)])  # 使用不同的种子
end_time = time.time()
print("优化后的并行生成时间: ", end_time - start_time)

# 检查结果是否重复
print("结果是否重复: ", np.array_equal(results[0], results[1]))

随机数安全性矩阵

1.23.3.2 随机种子的管理和传播

为了确保并行生成的随机数是唯一的，可以通过管理随机种子来避免重复生成。

import numpy as np
import multiprocessing

# 随机种子生成器
def seed_generator(initial_seed, num_seeds):
    """
    生成唯一的随机种子

    :param initial_seed: 初始随机种子
    :param num_seeds: 需要生成的种子数量
    :return: 生成的种子列表
    """
    np.random.seed(initial_seed)
    return np.random.randint(0, 2**32, size=num_seeds)

# 未优化的并行生成
def generate_random_unoptimized(seed):
    """
    生成随机数

    :param seed: 随机种子
    :return: 生成的随机数组
    """
    np.random.seed(seed)
    return np.random.rand(100000)

# 优化后的并行生成
def generate_random_optimized(seed):
    """
    生成随机数，使用唯一的种子

    :param seed: 随机种子
    :return: 生成的随机数组
    """
    np.random.seed(seed)
    return np.random.rand(100000)

# 生成种子
seeds = seed_generator(1, 4)

# 未优化的并行生成
start_time = time.time()
with multiprocessing.Pool() as pool:
    results_unoptimized = pool.map(generate_random_unoptimized, [1] * 4)  # 使用相同的种子
end_time = time.time()
print("未优化的并行生成时间: ", end_time - start_time)

# 检查结果是否重复
print("未优化的结果是否重复: ", np.array_equal(results_unoptimized[0], results_unoptimized[1]))

# 优化后的并行生成
start_time = time.time()
with multiprocessing.Pool() as pool:
    results_optimized = pool.map(generate_random_optimized, seeds)  # 使用不同的种子
end_time = time.time()
print("优化后的并行生成时间: ", end_time - start_time)

# 检查结果是否重复
print("优化后的结果是否重复: ", np.array_equal(results_optimized[0], results_optimized[1]))

1.23.4 初始化在数值计算中的陷阱

在数值计算中，初始化数据时可能会遇到一些陷阱，尤其是与浮点数精度相关的问题。本节将介绍常见的数值计算陷阱，并提供相应的解决方法。

1.23.4.1 浮点数精度累积误差案例

浮点数精度累积误差是指在多次浮点数运算中误差逐渐累积，最终影响计算结果的正确性。

import numpy as np

# 创建一个浮点数数组
data = np.array([0.1] * 10, dtype=np.float32)

# 累加操作
sum_data = np.sum(data)

# 打印结果
print("累加结果: ", sum_data)  # 期望结果是 1.0，但由于浮点数精度误差，结果可能略微不同

误差传播公式

对于n次加法运算，误差上界：
$$E_{\text{total}}\le {\gamma }_n\sum _{i=1}^n|x_i|$$
其中机器精度${\gamma }_n=\frac{nu}{1-nu}$，u为单元舍入误差

1.23.4.2 浮点数精度问题的解决方法

为了解决浮点数精度累积误差的问题，可以使用更高精度的数据类型，或者使用专用的数值计算库如 mpmath。

import numpy as np
import mpmath

# 创建一个浮点数数组
data = np.array([0.1] * 10, dtype=np.float64)

# 累加操作
sum_data = np.sum(data)

# 打印结果
print("使用 float64 累加结果: ", sum_data)  # 期望结果是 1.0

# 使用 mpmath
mpmath.mp.dps = 15  # 设置精度
mp_data = [mpmath.mpf(0.1)] * 10
sum_mp_data = sum(mp_data)

# 打印结果
print("使用 mpmath 累加结果: ", sum_mp_data)  # 期望结果是 1.0

1.23.5 初始化参数的自动化测试

在开发过程中，初始化参数的正确性至关重要。本节将介绍如何使用自动化测试工具来测试初始化参数，并提供相应的测试策略。

1.23.5.1 自动化测试工具的选择

常用的自动化测试工具有 pytest 和 unittest。pytest 以其简洁的语法和强大的功能成为首选。

import numpy as np
import pytest

# 初始化函数
def initialize_array(size, initial_value):
    """
    初始化数组

    :param size: 数组大小
    :param initial_value: 初始值
    :return: 初始化后的数组
    """
    return np.full(size, initial_value, dtype=np.float32)

# 测试初始化函数
def test_initialize_array():
    size = 10
    initial_value = 5.0
    result = initialize_array(size, initial_value)
    
    # 检查数组大小
    assert result.shape == (size,), "数组大小不匹配"
    
    # 检查初始值
    assert np.all(result == initial_value), "数组初始值不匹配"

# 运行测试
pytest.main(['-v', '1.23_数据工厂：高级初始化模式解析.md'])

1.23.5.2 初始化参数的测试策略

测试初始化参数时，可以采用以下策略：

1. 边界值测试：测试数组大小为0、1、最大值等边界情况。
2. 类型测试：测试不同数据类型的初始化。
3. 并行测试：测试并行初始化的正确性。

import numpy as np
import pytest

# 初始化函数
def initialize_array(size, initial_value):
    """
    初始化数组

    :param size: 数组大小
    :param initial_value: 初始值
    :return: 初始化后的数组
    """
    return np.full(size, initial_value, dtype=np.float32)

# 测试初始化函数
def test_initialize_array_boundaries():
    # 测试边界值
    size = 0
    initial_value = 5.0
    result = initialize_array(size, initial_value)
    assert result.shape == (0,), "数组大小不匹配"

    size = 1
    initial_value = 5.0
    result = initialize_array(size, initial_value)
    assert result.shape == (1,), "数组大小不匹配"
    assert result[0] == initial_value, "数组初始值不匹配"

def test_initialize_array_types():
    # 测试不同类型
    size = 10
    initial_value = 5
    result = initialize_array(size, initial_value)
    assert result.dtype == np.float32, "数组数据类型不匹配"
    assert np.all(result == np.float32(initial_value)), "数组初始值不匹配"

def test_initialize_array_parallel():
    # 测试并行初始化
    size = 1000000
    initial_value = 5.0
    results = []

    def parallel_initialize(size, initial_value):
        results.append(initialize_array(size, initial_value))

    with multiprocessing.Pool() as pool:
        pool.map(parallel_initialize, [(size, initial_value)] * 4)

    assert all(np.all(result == initial_value) for result in results), "并行初始化结果不匹配"

# 运行测试
pytest.main(['-v', '1.23_数据工厂：高级初始化模式解析.md'])

测试覆盖率报告

总结

通过本篇文章的详细讲解和示例，我们对 NumPy 中的高级初始化模式有了更深入的理解。主要内容包括：

1. 从数据库初始化的批量技巧：

   • SQL 到 ndarray 的高效转换
   • 批量读取优化

2. 内存映射初始化大数组方案：

   • 100GB 数组的延迟初始化
   • 内存映射的性能优势

3. 随机初始化种子传播机制：

   • 并行随机数生成的安全性
   • 随机种子的管理和传播

4. 初始化在数值计算中的陷阱：

   • 浮点数精度累积误差案例
   • 浮点数精度问题的解决方法

5. 初始化参数的自动化测试：

   • 自动化测试工具的选择
   • 初始化参数的测试策略

代码参考汇总

为了方便参考，以下是参考代码示例：

# 从数据库初始化
import numpy as np
import sqlite3

def initialize_from_db(db_path, query, dtype=np.float32):
    conn = sqlite3.connect(db_path)
    cursor = conn.cursor()
    cursor.execute(query)
    result = cursor.fetchall()
    conn.close()
    return np.array(result, dtype=dtype)

# 内存映射初始化
def initialize_large_array_mmap(size, file_path, dtype=np.float32):
    return np.memmap(file_path, dtype=dtype, mode='w+', shape=size)

# 并行随机数生成
import numpy as np
import multiprocessing

def generate_random(seed):
    rng = np.random.RandomState(seed)
    return rng.rand(100000)

def seed_generator(initial_seed, num_seeds):
    np.random.seed(initial_seed)
    return np.random.randint(0, 2**32, size=num_seeds)

def test_parallel_random_generation():
    initial_seed = 1
    num_seeds = 4
    seeds = seed_generator(initial_seed, num_seeds)

    with multiprocessing.Pool() as pool:
        results = pool.map(generate_random, seeds)

    for i in range(len(results)):
        for j in range(i + 1, len(results)):
            assert not np.array_equal(results[i], results[j]), "并行生成的结果重复"

# 浮点数精度问题
def test_float_precision():
    data = np.array([0.1] * 10, dtype=np.float32)
    sum_data = np.sum(data)
    print("累加结果: ", sum_data)

    # 使用更高精度的数据类型
    data_high_precision = np.array([0.1] * 10, dtype=np.float64)
    sum_data_high_precision = np.sum(data_high_precision)
    print("使用 float64 累加结果: ", sum_data_high_precision)

# 性能测试
def test_performance(size=1000000):
    seed = 1

    # 未优化的生成
    start_time = time.time()
    result_unoptimized = np.random.rand(size)
    end_time = time.time()
    print("未优化的生成时间: ", end_time - start_time)

    # 优化后的生成
    start_time = time.time()
    result_optimized = generate_random_numba(size, seed)
    end_time = time.time()
    print("优化后的生成时间: ", end_time - start_time)

# 存储和读取大数组
def test_hdf5_storage():
    large_array = np.random.rand(100000, 100000)

    # 存储大数组
    start_time = time.time()
    with h5py.File('large_array.h5', 'w') as f:
        f.create_dataset('data', data=large_array)
    end_time = time.time()
    print("存储时间: ", end_time - start_time)

    # 读取大数组
    start_time = time.time()
    with h5py.File('large_array.h5', 'r') as f:
        loaded_array = f['data'][:]
    end_time = time.time()
    print("读取时间: ", end_time - start_time)

    assert np.array_equal(large_array, loaded_array), "存储和读取的数据不一致"

# 数组切片和索引
def test_advanced_indexing():
    data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

    # 布尔索引
    bool_index = data > 4
    filtered_data = data[bool_index]

    # 花式索引
    row_indices = [0, 2]
    column_indices = [1, 2]
    fancy_index_data = data[row_indices, column_indices]

    print("布尔索引: \n", bool_index)
    print("过滤后的数据: \n", filtered_data)
    print("花式索引: \n", fancy_index_data)

# 广播机制
def test_broadcasting():
    data = np.array([[1, 2, 3], [4, 5, 6]])
    scalar = np.array([10, 20])
    result = data + scalar

    print("原始数组: \n", data)
    print("一维数组: \n", scalar)
    print("广播结果: \n", result)

# 通用函数（ufunc）
def test_ufuncs():
    data = np.array([1, 2, 3, 4, 5])
    sin_data = np.sin(data)
    exp_data = np.exp(data)
    sqrt_data = np.sqrt(data)

    print("原始数组: ", data)
    print("正弦结果: ", sin_data)
    print("指数结果: ", exp_data)
    print("平方根结果: ", sqrt_data)

if __name__ == "__main__":
    import pytest
    pytest.main(['-v', '1.23_数据工厂：高级初始化模式解析.md'])

互动练习

为了巩固所学知识，建议你完成以下互动练习：

1. 练习1：从一个更大的数据库表中初始化一个 ndarray，并进行一些基本的统计运算。
2. 练习2：使用内存映射技术初始化一个 1TB 的数组，并进行简单的读写操作。
3. 练习3：实现一个并行随机数生成函数，并测试其安全性和性能。
4. 练习4：编写一个函数，使用 Kahan 算法来累加一个包含 100 万个 0.1 的数组。
5. 练习5：编写一个函数，使用 HDF5 存储和读取一个 100GB 的数组，并测试其性能。

希望你通过这些练习能够进一步提升你的 NumPy 技能。祝你学习顺利！

参考文献

1. NumPy 官方文档：NumPy Documentation
2. 《Python 数据科学 handbook》：Jake VanderPlas
3. 《高性能 Python》：Micha Gorelick 和 Ian Ozsvald
4. HDF5 官方文档：HDF5 Documentation

这篇文章包含了详细的原理介绍、代码示例、源码注释以及案例等。希望这对您有帮助。如果有任何问题请随私信或评论告诉我。