管理大数据应用程序的云存储成本
降低使用基于云的存储成本的技巧
·
关注 发表在 Towards Data Science · 11 分钟阅读 · 2023 年 6 月 26 日
–
图片由 JOSHUA COLEMAN 提供，来源于 Unsplash
随着对不断增长的数据量的依赖，现代公司比以往任何时候都更加依赖高容量和高度可扩展的数据存储解决方案。对于许多公司来说，这种解决方案通常表现为云存储服务，如Amazon S3、Google Cloud Storage和Azure Blob Storage，这些服务都提供了丰富的 API 和功能（例如，多层存储），支持各种数据存储设计。当然，云存储服务也有相关的成本。这些成本通常包括多个组成部分，包括你使用的存储空间的总体大小，以及如将数据传入、传出或在云存储中传输等活动。例如，Amazon S3 的价格（截至本文撰写时）包括六个成本组成部分，每个部分都需要考虑。可以看出，管理云存储成本可能变得复杂，因此开发了指定的计算器（例如，这里）来协助处理这一问题。
在一篇近期文章中，我们扩展了设计数据和数据使用的重要性，以降低与数据存储相关的成本。我们在那里关注的是使用数据压缩作为减少数据总体大小的一种方法。在这篇文章中，我们关注的是一个有时被忽视的云存储成本组成部分 —— 对你的云存储桶和数据对象进行的 API 请求成本。我们将通过示例展示为什么这个组成部分常常被低估，以及如果管理不当，它如何可能成为你大数据应用成本的重要组成部分。然后我们将讨论几种简单的方法来控制这一成本。
免责声明
尽管我们的演示将使用 Amazon S3，但本文的内容同样适用于任何其他云存储服务。请不要将我们选择 Amazon S3 或我们提到的任何其他工具、服务或库解读为对其使用的支持。最适合你的选项将取决于你自己项目的具体细节。此外，请记住，关于如何存储和使用数据的任何设计选择都会有其利弊，这些利弊应根据你自己项目的细节权衡。
本文将包括在 Amazon EC2 c5.4xlarge 实例（具有 16 个 vCPUs 和 “高达 10 Gbps” 的网络带宽）上运行的一系列实验。我们将分享这些实验的输出，作为你可能看到的比较结果的示例。请注意，输出结果可能会因实验运行的环境而有很大不同。请不要将这里呈现的结果作为你自己设计决策的依据。我们强烈建议你在决定自己项目的最佳方案之前，运行这些以及其他额外的实验。
一个简单的思考实验
假设你有一个数据转换应用程序，处理来自 S3 的 1 MB 数据样本，并生成 1 MB 数据输出，然后将其上传到 S3。假设你的任务是通过在合适的 Amazon EC2 实例 上运行你的应用程序来转换 10 亿个数据样本（与 S3 存储桶位于同一地区，以避免数据传输成本）。现在假设 Amazon S3 收费 每 1000 次 GET 操作收费 $0.0004，每 1000 次 PUT 操作收费 $0.005（截至本文撰写时）。乍一看，这些成本似乎如此低，以至于与数据转换相关的其他成本相比，几乎可以忽略不计。然而，简单的计算显示，单是我们的 Amazon S3 API 调用就会产生**$5,400**的费用!! 这可能是你项目中最主要的成本因素，甚至比计算实例的成本还要高。我们将在本文最后回到这个思考实验。
将数据批量归为大文件
降低 API 调用成本的显而易见方法是将样本归为大文件，并在样本批次上进行转换。将我们的批量大小记作N，这种策略可能将我们的成本降低N倍（假设没有使用多部分文件传输——见下文）。这种技术不仅在 PUT 和 GET 调用上节省了费用，还能在所有与对象文件数量相关的 Amazon S3 成本组成部分上节省开支，而不是数据的整体大小（例如，生命周期转换请求)。
将样本组合在一起有一些缺点。例如，当你单独存储样本时，你可以随时访问其中的任何一个。当样本组合在一起时，这就变得更加具有挑战性。（有关将样本批处理成大文件的优缺点，请参见这篇文章）。如果你选择将样本组合在一起，那么最大的疑问是如何选择大小 N。较大的 N 可能会降低存储成本，但可能会引入延迟，增加计算时间，并由此增加计算成本。找到最佳数量可能需要一些实验，考虑到这些以及其他额外的因素。
但我们不要自欺欺人。进行这种更改并不容易。你的数据可能有很多消费者（包括人类和人工智能），每个人都有自己特定的需求和限制。将样本存储在单独的文件中可以更容易地让每个人满意。找到一个能够满足所有人需求的批处理策略将会很困难。
可能的折衷方案：批量 PUT，单独 GET
你可以考虑的一种折衷方案是上传包含分组样本的大文件，同时允许访问单个样本。实现这一点的一种方法是维护一个索引文件，其中包含每个样本的位置（其所在的文件、起始偏移量和结束偏移量），并向每个消费者提供一个轻量的 API 层，以便他们可以自由下载单个样本。该 API 将使用索引文件和一个允许从对象文件中提取特定范围的 S3 API（例如，Boto3 的 get_object 函数）来实现。虽然这种解决方案不会节省 GET 调用的费用（因为我们仍然要提取相同数量的单个样本），但由于我们上传的是较少的较大文件，因此更昂贵的 PUT 调用将会减少。请注意，这种解决方案对我们与 S3 交互所使用的库有一些限制，因为它依赖于一个允许从大型文件对象中提取部分块的 API。在之前的帖子中（例如，这里），我们讨论了与 S3 接口的不同方式，其中许多方式不支持此功能。
下面的代码块演示了如何实现一个简单的 PyTorch 数据集（使用 PyTorch 版本 1.13），该数据集使用 Boto3 get_object API 从包含分组样本的大文件中提取单个 1 MB 样本。我们将这种方式迭代数据的速度与迭代存储在单独文件中的样本的速度进行比较。
import os, boto3, time, numpy as np
import torch
from torch.utils.data import Dataset
from statistics import mean, variance
KB = 1024
MB = KB * KB
GB = KB ** 3
sample_size = MB
num_samples = 100000
# modify to vary the size of the files
samples_per_file = 2000 # for 2GB files
num_files = num_samples//samples_per_file
bucket = '<s3 bucket>'
single_sample_path = '<path in s3>'
large_file_path = '<path in s3>'
class SingleSampleDataset(Dataset):
def __init__(self):
super().__init__()
self.bucket = bucket
self.path = single_sample_path
self.client = boto3.client("s3")
def __len__(self):
return num_samples
def get_bytes(self, key):
response = self.client.get_object(
Bucket=self.bucket,
Key=key
)
return response['Body'].read()
def __getitem__(self, index: int):
key = f'{self.path}/{index}.image'
image = np.frombuffer(self.get_bytes(key),np.uint8)
return {"image": image}
class LargeFileDataset(Dataset):
def __init__(self):
super().__init__()
self.bucket = bucket
self.path = large_file_path
self.client = boto3.client("s3")
def __len__(self):
return num_samples
def get_bytes(self, file_index, sample_index):
response = self.client.get_object(
Bucket=self.bucket,
Key=f'{self.path}/{file_index}.bin',
Range=f'bytes={sample_index*MB}-{(sample_index+1)*MB-1}'
)
return response['Body'].read()
def __getitem__(self, index: int):
file_index = index // num_files
sample_index = index % samples_per_file
image = np.frombuffer(self.get_bytes(file_index, sample_index),
np.uint8)
return {"image": image}
# toggle between single sample files and large files
use_grouped_samples = True
if use_grouped_samples:
dataset = LargeFileDataset()
else:
dataset = SingleSampleDataset()
# set the number of parallel workers according to the number of vCPUs
dl = torch.utils.data.DataLoader(dataset, shuffle=True,
batch_size=4, num_workers=16)
stats_lst = []
t0 = time.perf_counter()
for batch_idx, batch in enumerate(dl, start=1):
if batch_idx % 100 == 0:
t = time.perf_counter() - t0
stats_lst.append(t)
t0 = time.perf_counter()
mean_calc = mean(stats_lst)
var_calc = variance(stats_lst)
print(f'mean {mean_calc} variance {var_calc}')
下表总结了不同样本分组大小 N 的数据遍历速度。
不同分组策略对数据遍历时间的影响（作者）
注意，尽管这些结果强烈暗示将样本分组到大文件中对单独提取它们的性能影响相对较小，但我们发现比较结果会根据样本大小、文件大小、文件偏移量的值、对同一文件的并发读取次数等有所变化。虽然我们无法了解 Amazon S3 服务的内部工作原理，但考虑到内存大小、内存对齐和限制等因素对性能的影响也不足为奇。找到适合你数据的最佳配置可能需要一些实验。
一个可能干扰我们在这里描述的节省成本的分组策略的重要因素是多部分下载和上传的使用，我们将在下一节中讨论。
使用可以控制多部分数据传输的工具
许多云存储服务提供商支持对象文件的多部分上传和下载选项。在多部分数据传输中，大于某个阈值的文件被划分为多个部分并同时传输。如果你想加速大文件的数据传输，这是一项关键功能。AWS 建议对大于 100 MB 的文件使用多部分上传。在以下简单示例中，我们比较了将 2 GB 文件的多部分 阈值 和 分块大小 设置为不同值时的下载时间：
import boto3, time
KB = 1024
MB = KB * KB
GB = KB ** 3
s3 = boto3.client('s3')
bucket = '<bucket name>'
key = '<key of 2 GB file>'
local_path = '/tmp/2GB.bin'
num_trials = 10
for size in [8*MB, 100*MB, 500*MB, 2*GB]:
print(f'multi-part size: {size}')
stats = []
for i in range(num_trials):
config = boto3.s3.transfer.TransferConfig(multipart_threshold=size,
multipart_chunksize=size)
t0 = time.time()
s3.download_file(bucket, key, local_path, Config=config)
stats.append(time.time()-t0)
print(f'multi-part size {size} mean {mean(stats)}')
此实验的结果总结在下面的表格中：
多部分分块大小对下载时间的影响（作者）
请注意，相对比较将大大依赖于测试环境，特别是实例与 S3 存储桶之间的通信速度和带宽。我们的实验是在与存储桶位于同一区域的实例上运行的。然而，随着距离的增加，使用多部分下载的影响也会增加。
关于我们讨论的话题，重要的是要注意多部分数据传输的成本影响。具体来说，当你使用多部分数据传输时，你会为每个文件部分的 API 操作收费。因此，使用多部分上传/下载将限制将数据样本批量处理成大文件的成本节省潜力。
许多 API默认使用多部分下载。如果你的主要关注点是减少与 S3 的交互延迟，这种做法是很好的。但是，如果你的关心是限制成本，那么这种默认行为对你来说并不有利。例如，Boto3 是一个流行的 Python API，用于从 S3 上传和下载文件。如果未指定，boto3 的 S3 API 如 upload_file 和 download_file 将使用默认的 TransferConfig，该配置将对任何大于 8 MB 的文件应用 8 MB 的块大小的多部分上传/下载。如果你负责控制组织中的云成本，你可能会惊讶地发现这些 API 被广泛使用其默认设置。在许多情况下，你可能会发现这些设置是不必要的，增加多部分阈值和块大小值，或完全禁用多部分数据传输，对应用程序的性能影响很小。
示例 — 多部分文件传输大小对速度和成本的影响
在下面的代码块中，我们创建一个简单的多进程转换函数，并测量多部分块大小对其性能和成本的影响：
import os, boto3, time, math
from multiprocessing import Pool
from statistics import mean, variance
KB = 1024
MB = KB * KB
sample_size = MB
num_files = 64
samples_per_file = 500
file_size = sample_size*samples_per_file
num_processes = 16
bucket = '<s3 bucket>'
large_file_path = '<path in s3>'
local_path = '/tmp'
num_trials = 5
cost_per_get = 4e-7
cost_per_put = 5e-6
for multipart_chunksize in [1*MB, 8*MB, 100*MB, 200*MB, 500*MB]:
def empty_transform(file_index):
s3 = boto3.client('s3')
config = boto3.s3.transfer.TransferConfig(
multipart_threshold=multipart_chunksize,
multipart_chunksize=multipart_chunksize
)
s3.download_file(bucket,
f'{large_file_path}/{file_index}.bin',
f'{local_path}/{file_index}.bin',
Config=config)
s3.upload_file(f'{local_path}/{file_index}.bin',
bucket,
f'{large_file_path}/{file_index}.out.bin',
Config=config)
stats = []
for i in range(num_trials):
with Pool(processes=num_processes) as pool:
t0 = time.perf_counter()
pool.map(empty_transform, range(num_files))
transform_time = time.perf_counter() - t0
stats.append(transform_time)
num_chunks = math.ceil(file_size/multipart_chunksize)
num_operations = num_files*num_chunks
transform_cost = num_operations * (cost_per_get + cost_per_put)
if num_chunks > 1:
# if multi-part is used add cost of
# CreateMultipartUpload and CompleteMultipartUpload API calls
transform_cost += 2 * num_files * cost_per_put
print(f'chunk size {multipart_chunksize}')
print(f'transform time {mean(stats)} variance {variance(stats)}
print(f'cost of API calls {transform_cost}')
在这个例子中，我们将文件大小固定为 500 MB，并对下载和上传应用了相同的多部分设置。一个更完整的分析将会变化数据文件的大小和多部分设置。
在下表中，我们总结了实验的结果。
多部分块大小对数据转换速度和成本的影响（作者）
结果表明，至多达到 500 MB 的多部分块大小（我们文件的大小），数据转换时间的影响是最小的。另一方面，与使用 Boto3 的默认块大小（8MB）相比，云存储 API 成本的潜在节省是显著的，最多可达 98.4%。这个例子不仅展示了将样本分组的成本效益，而且还暗示了通过适当配置多部分数据传输设置获得额外节省的机会。
结论
让我们将上一个示例的结果应用到我们在本文开头介绍的思维实验中。我们展示了如果数据样本存储在单独的文件中，对 1 亿个数据样本应用简单转换的费用将是$5,400。如果我们将样本分组为 200 万个文件，每个文件有 500 个样本，并且在不使用多部分数据传输的情况下（如上述示例的最后一次试验），API 调用的成本将减少到$10.8！！同时，假设测试环境相同，我们预期的整体运行时间的影响（基于我们的实验）将相对较低。我认为这是一个相当不错的交易。你觉得呢？
摘要
在开发基于云的大数据应用时，了解我们活动的所有成本细节至关重要。在这篇文章中，我们专注于 Amazon S3 定价策略中的“请求与数据检索”组件。我们演示了这个组件如何成为大数据应用总体成本的重要组成部分。我们讨论了可能影响这一成本的两个因素：数据样本的分组方式和多部分数据传输的使用方式。
自然，仅优化一个成本组件可能会以某种方式增加其他组件，从而提高整体成本。数据存储的适当设计需要考虑所有潜在成本因素并且将大大依赖于你的具体数据需求和使用模式。
一如既往，请随时提出评论和修正。
机器学习系统的技术债务管理
探索用于可持续减轻快速交付成本的实践（设计模式、版本控制和监控系统）——包括实现代码
·
关注 发表在 Towards Data Science ·10 分钟阅读·2023 年 9 月 26 日
–
随着机器学习（ML）社区的进步，开发 ML 项目的资源非常丰富。例如，我们可以依赖通用的 Python 包 scikit-learn，该包建立在 NumPy、SciPy 和 matplotlib 之上，用于数据预处理和基本预测任务。或者，我们可以利用来自 Hugging Face 的预训练模型的开源集合来分析各种类型的数据集。这些工具使得当前的数据科学家能够迅速而轻松地处理标准 ML 任务，同时取得适度良好的模型性能。
然而，ML 工具的丰富性往往使业务利益相关者甚至从业者低估了构建企业级 ML 系统所需的努力。尤其是在面对紧迫的项目截止日期时，团队可能会加快将系统部署到生产环境的速度，而没有充分考虑技术因素。因此，ML 系统往往不能以技术上可持续和可维护的方式满足业务需求。
随着系统的发展和部署，技术债务会不断累积——隐含的成本拖延不解决，修复的成本会越来越高。
在机器学习（ML）系统中存在多种技术债务来源。以下是一些例子。
#1 硬性代码设计以应对不可预见的需求
为了验证 ML 是否能解决当前的企业挑战，许多 ML 项目开始时会进行 概念验证（PoC）。我们最初创建了一个 Jupyter Notebook 或 Google Colab 环境来探索数据，然后开发了几个临时函数，给利益相关者制造了项目即将完成的假象。直接从 PoC 构建的系统最终可能大部分由 胶水代码 构成——这是一种连接特定不兼容组件的支持代码，但本身并不具备数据分析的功能。它们可能像意大利面一样，难以维护且容易出错。
Jakob Owens 拍摄于 Unsplash
业务利益相关者经常提出不同的要求或希望扩大项目，例如尝试新的数据源或新的算法。因此，我们经常发现自己频繁地重新访问涵盖当前预处理管道和模型开发过程的代码库。灵活性差的代码设计可能导致对新变化的反应困难，甚至需要重写大部分代码以进行小的调整。
#2 ML 系统配置中的混乱
软件工程编程通过定义计算机执行的规则来自动化任务，确保相同输入的精确输出。软件工程师还关注每个角落情况的正确性。另一方面，ML 系统通过收集和输入特征数据到模型中来自动化任务，以实现期望的目标结果。这个实验过程接受不确定性和可变性。随着 ML 系统的成熟，它们通常包含多个版本的配置选项，如具有不同特征组合的数据集和特定算法的学习设置。
照片来自 Ricardo Viana 在 Unsplash
ML 系统中的输入特征本质上是相互关联的。考虑一种情况，如果特征 A 在生产环境中不再可用，你需要重新评估剩余特征的权重。然而，2 个月后，特征 A 又变得可用。如果你没有系统地保存或甚至错误地修改了原始配置，恢复性能下降将需要额外的计算资源和时间精力。
#3 适应不断变化的外部世界的能力有限
ML 系统通常依赖于外部世界，各种隐藏因素不断演变，但未得到适当考虑和监控。
导致模型性能潜在下降的外部因素（图像来源：作者）
-
来自上游生产者的不稳定数据输出：我们的 ML 系统的输入信号可能来自一个随着时间更新的其他机器学习模型。此外，系统可能依赖于来自物联网（IoT）设备的信号、网页抓取数据或音频到文本转换器的输出数据。如果上游生产者的这些工具的维护没有得到适当声明或实施了有缺陷的修补，它可能会意外地降低 ML 系统的性能。
-
输入数据的漂移：以零售行业的需求预测为例。输入数据可以周期性地（例如购买行为的季节性循环）、渐进性地（如供应商商品的通货膨胀成本）和突然地（新竞争者的进入）表现出新的分布。
在接下来的部分中，我们将深入探讨构建机器学习系统的一些最佳实践，并提供示例 Python 代码以演示其实现。
想象一下你希望为城市建立一个强大的交通系统，准备应对交通高峰，因此你收集了过去两年的传感器交通数据。你的目标是预测未来六个月的交通模式（即车辆数量）。
-
训练数据集：ID、日期时间以及车辆数量
-
测试数据集：ID 和日期时间
图片来源：Joey Kyber 在 Unsplash
#1 对于机器学习代码库使用设计模式
为了使代码设计在未来需求中更灵活和可重用，我们可以利用设计模式。这些模式作为解决各种情况中常见问题的模板，使我们能够解耦代码库的不同部分。因此，这有助于提高对代码库的理解，并建立一种共同的语言，以快速沟通解决方案。
机器学习项目中的两个主要组成部分是数据和算法，这些部分可以从设计模式中受益。
- 工厂模式
这种创建模式为在运行时生成对象提供了一层抽象。在机器学习系统中，我们可以通过创建一个数据加载类（例如CSVDataLoader）来实现这种模式，该类处理训练/测试数据的加载、保存和返回，并保持一致的数据类型。然后我们可以声明DataProcessor接口，而不指定具体实现。
# CSV Data loader class
class CSVDataLoader:
def __init__(self, file_path):
self.file_path = file_path
def get_data(self):
return pd.read_csv(self.file_path)
def save_data(self, df):
df.to_csv(f'data/transformed_{self.file_path}', index=False)
# Interface
class DataProcessor:
def __init__(self, train_data_loader, test_data_loader):
self.train_data_loader = train_data_loader
self.test_data_loader = test_data_loader
def run(self):
# Load training and test data using data loaders
train_df = self.train_data_loader.get_data()
test_df = self.test_data_loader.get_data()
# Create a data processor instance
process = DataProcessor(
train_data_loader=CSVDataLoader(file_path='train.csv'),
test_data_loader=CSVDataLoader(file_path='test.csv')
)
# Run the data processing pipeline
process.run()
这种方法允许你扩展代码，而不必重新实现DataProcessor。例如，如果你想从 JSON 文件加载数据集，你可以简单地创建一个新的类JSONDataLoader作为数据加载器的实例进行声明。
- 策略模式
由于没有适用于所有机器学习问题的通用算法，我们经常在原型设计或项目增强过程中切换和试验不同的算法。我们可以应用策略模式，通过创建一个新的类DataTransformer进行特征工程，以及另一个类LGBMModel来封装使用 LightGBM 模型进行拟合和预测的策略。
class DataTransformer:
def transform_data(self, train_df, test_df):
for idx, df in enumerate([train_df, test_df]):
df[‘DateTime’] = pd.to_datetime(df['DateTime'])
# Build ‘Time’ column
df['Time'] = [date.hour * 3600 + date.minute * 60 + date.second for date in df['DateTime']]
# Convert DateTime to Unix timestamp
unixtime = [time.mktime(date.timetuple()) for date in df['DateTime']]
df['DateTime'] = unixtime
# Perform one-hot encoding on the DataFrame
df = pd.get_dummies(df)
if idx == 0:
# Split training DataFrame into features (X_train) and target (y_train)
X_train = df.drop(['Vehicles'], axis=1)
y_train = df[['Vehicles']]
elif idx == 1:
# Store test DataFrame
X_test = df
return X_train, y_train, X_test
class LGBMModel:
def __init__(self, num_leaves, n_estimators):
self.model = lgb.LGBMRegressor(
num_leaves=num_leaves,
n_estimators=n_estimators
)
def fit(self, X, y):
self.model.fit(X, y)
self.model.booster_.save_model('model/lgbm_model.txt')
return self
def predict(self, X):
predictions = self.model.predict(X)
return predictions
接口DataProcessor的实现和声明如下。这是一个端到端的过程，包括分别使用train_data_loader和test_data_loader加载训练和测试数据，使用data_transformer转换数据，并使用model将模型拟合到转换后的数据。因此，我们可以预测测试数据集中每条记录中的车辆数量。
# Interface
class DataProcessor:
def __init__(self, train_data_loader, test_data_loader, data_transformer, model):
self.train_data_loader = train_data_loader
self.test_data_loader = test_data_loader
self.data_transformer = data_transformer
self.model = model
def run(self):
# Load train and test data using data loaders
train_df = self.train_data_loader.get_data()
test_df = self.test_data_loader.get_data()
# Transform the data using the data transformer
X_train, y_train, X_test = self.data_transformer.transform_data(train_df, test_df)
# Fit the model and make prediction
self.model.fit(X_train, y_train)
test_df['Vehicles'] = self.model.predict(X_test)
# Save the transformed training data and test data
self.train_data_loader.save_data(pd.concat([X_train, y_train], axis=1))
self.test_data_loader.save_data(test_df)
# Create a data processor instance
process = DataProcessor(
train_data_loader=CSVDataLoader(file_path='train.csv'),
test_data_loader=CSVDataLoader(file_path='test.csv'),
data_transformer=DataTransformer(),
model=LGBMModel(num_leaves=16, n_estimators=80)
)
# Run the data processing pipeline
process.run()
你可以轻松地添加新的代码块来实现其他数据转换想法或算法。类似于工厂模式，这些更改不需要你修改接口DataProcessor。这种设计使得即使有很长的算法列表，也更容易维护代码。ML 系统的行为可以根据选择的策略动态变化。
当然，上述代码实现仅是开发的初步模板。例如，我们可以进一步通过涵盖数据验证、实施超参数调优机制以及评估模型来增强代码。
#2 ML 系统的版本控制
在复杂的模型开发和管理过程中，我们需要适当的版本控制。这使我们能够维护自己或团队成员所做的修改历史，并跟踪本地环境中相对于 ML 系统组件的数据、训练模型和超参数的版本。因此，我们可以解决一些常见问题，包括：
-
是什么更改导致了模型的失败？
-
哪些修改导致了模型性能的提升？
-
最近发布的是哪个版本的模型？
在这里，我们展示了如何利用DVC中的版本控制功能，它在Git库中效果最佳，用于跟踪我们的原始交通数据、转换后的交通数据和 LGBM 模型。
# Initialise a Git and DVC project in the current working directory
git init
dvc init
# Capture the current state of transformed data in folder 'data' and latest LGBM model in folder 'model'
dvc add data model
# Commit the current state of 1st version
git add data.dvc model.dvc .gitignore
git commit -m “First LGBM model, with Time feature”
git tag -a “v1.0” -m “model v1.0, Time feature”
让我们考虑一个场景，我们在第二个版本的数据处理过程中进行了以下更改：
- 在
DataTransformer类中添加 Weekday 特性
df[‘Weekday’] = [datetime.weekday(date) for date in df.DateTime]
- 在
DataProcessor接口中设置 LGBM 模型的新配置参数
model=LGBMModel(num_leaves=20, n_estimators=90)
使用以下命令，我们可以在 DVC 中跟踪数据和模型的第二个版本，并用 Git 提交指向它们的 .dvc 文件。
git add data.dvc model.dvc
git commit -m “Second LGBM model, with Time and Weekday features”
git tag -a “v2.0” -m “model v2.0, Time and Weakday features”
尽管工作区当前定位于数据和模型的第二个版本，但我们可以在必要时轻松地切换回并恢复到第一个快照。
git checkout v1.0
dvc checkout
上述命令涵盖了基本操作。我们可以进一步利用该工具进行项目组织和协作。使用案例的例子包括理解数据最初是如何构建的，并比较实验中的模型指标。
#3 持续测试和监控 ML 系统
一旦增强的模型能够生成预测，必须在将其投入生产之前执行 合理性检查。这可以通过将一组随机的在线数据集拟合到最新模型中进行离线检查，并从不同角度审视结果来实现。
-
确保正确的访问权限：模型结果可以存储在目标路径中（例如，将它们写入 Hive 表）。
-
消除 语义错误：可视化模型拟合的变换特征的分布，以识别任何异常行为。
-
评估模型性能：使用最新模型重新评分，并使用适当的性能指标（例如，对于不平衡类问题，F1-score 是比准确度更好的衡量指标）比较结果。
即使最新版本的机器学习系统发布后，持续监控仍然是必要的，以应对不断变化的外部环境。
-
监控数据漂移和模型漂移：通过模型性能指标、统计测试和自适应窗口技术来检测漂移条件。
-
跟踪上游生产者：了解上游流程中的变化，并定期监控它们，以确保它们符合服务水平目标。
总结
我们已经探索了几种有效的实践，这些实践可以用来解决在开发和部署机器学习系统时出现的技术债务。
-
使用设计模式，创建一个模块化且灵活的数据处理管道，以适应不可预见的需求。
-
利用版本控制，跟踪和管理机器学习工件，如数据和模型，确保工作流程更少混乱。
-
测试和监控机器学习系统，以便及时顺利地处理动态外部世界中的变化。
离开之前
如果你喜欢这篇阅读，我邀请你关注我的 Medium 页面 和 Linkedin 页面。这样，你可以及时获得有关数据科学侧项目、机器学习运营（MLOps）演示以及项目管理方法的精彩内容。
我们的模型在不断变化的世界中如何受到影响？一个集中于漂移示例的分析，并实现基于 Python 的…
[towardsdatascience.com ## 利用超越 A/B 测试的方法优化策略
对经典 A/B 测试的深入解释：Epsilon-greedy、Thompson Sampling、Contextual…
[towardsdatascience.com
使用 Rclone 管理你的云数据存储
如何优化多个对象存储系统之间的数据传输
·
关注 发表于 Towards Data Science ·7 分钟阅读·2023 年 11 月 22 日
–
图片由 Tom Podmore 提供，来源于 Unsplash
随着公司越来越依赖基于云的存储解决方案，拥有合适的工具和技术以有效管理其大数据变得至关重要。在之前的文章中（例如，这里和这里），我们探讨了几种从云存储中检索数据的方法，并展示了它们在不同任务中的有效性。我们发现，最优的工具可能会根据具体任务（例如文件格式、数据文件的大小、数据访问模式）以及我们希望优化的指标（例如延迟、速度或成本）而有所不同。在这篇文章中，我们探讨了另一个流行的云存储管理工具——有时被称为“云存储的瑞士军刀”——rclone命令行工具。rclone 支持超过70 种存储服务提供商，具备类似于供应商特定存储管理应用程序（如 AWS CLI（用于 Amazon S3）和gsutil（用于 Google Storage））的功能。但它是否足够出色以构成一个可行的替代方案？在什么情况下 rclone 会是首选工具？在接下来的章节中，我们将展示 rclone 的使用，评估其性能，并突出其在特定用例中的价值——在不同对象存储系统之间转移数据。
免责声明
本文绝不是为了取代官方的rclone 文档。也不意图为使用 rclone 或我们提到的其他工具提供认可。您在进行云数据管理时的最佳选择将大大依赖于项目的详细信息，并应根据详细的、特定用例的测试来做出。请确保在阅读时重新评估我们所做的陈述，并对照您手头的最新工具。
从云存储中检索数据与 Rclone
以下命令行使用rclone sync来同步云对象存储路径与本地目录的内容。这个例子演示了Amazon S3存储服务的使用，但同样可以使用其他云存储服务。
rclone sync -P \
--transfers 4 \
--multi-thread-streams 4 \
S3store:my-bucket/my_files ./my_files
rclone 命令有数十个 标志 用于编程其行为。* -P * 标志输出数据传输的进度，包括传输速率和总体时间。在上述命令中，我们包含了两个（众多）可能影响 rclone 运行时性能的控制：transfers 标志确定要并发下载的最大文件数，而 multi-thread-streams 确定用于传输单个文件的最大线程数。这里我们将两者都保留在默认值（4）。
rclone 的功能依赖于对 rclone 配置文件 的适当定义。下面我们演示了上面命令行中使用的远程 S3store 对象存储位置的定义。
[S3store]
type = s3
provider = AWS
access_key_id = <id>
secret_access_key = <key>
region = us-east-1
现在我们已经看到 rclone 的实际操作，接下来的问题是它是否相较于其他云存储管理工具（例如流行的 AWS CLI）提供了任何额外的价值。在接下来的两个部分中，我们将评估 rclone 相比于其一些替代品在我们之前帖子中详细探讨的两个场景中的性能：1) 下载一个 2 GB 文件和 2) 下载数百个 1 MB 文件。
用例 1：下载大文件
以下命令行使用 AWS CLI 从 Amazon S3 下载一个 2 GB 文件。这只是我们在 之前的帖子 中评估的众多方法之一。我们使用 linux 的 time 命令来测量性能。
time aws s3 cp s3://my-bucket/2GB.bin .
报告的下载时间约为 26 秒（即 ~79 MB/s）。请注意，该值是在我们自己本地 PC 上计算的，可能会因运行时环境的不同而有很大差异。等效的 rclone copy 命令如下：
rclone sync -P S3store:my-bucket/2GB.bin .
在我们的设置中，我们发现 rclone 的下载时间比标准的 AWS CLI 慢两倍多。通过适当调整 rclone 控制标志，这一性能有可能得到显著改善。
用例 2：下载大量小文件
在这个用例中，我们评估了下载 800 个相对较小的 1 MB 文件的运行时性能。在 之前的博客文章 中，我们讨论了在将数据样本流式传输到深度学习训练工作负载的背景下的这个用例，并展示了 s5cmd beast 模式的优越性能。在 beast 模式下，我们创建一个包含对象文件操作列表的文件，s5cmd 使用 多个并行工作线程（默认值为 256）来执行这些操作。下面展示了 s5cmd beast 模式选项：
time s5cmd --run cmds.txt
cmds.txt 文件包含 800 行，格式如下：
cp s3://my-bucket/small_files/<i>.jpg <local_path>/<i>.jpg
s5cmd 命令平均耗时 9.3 秒（十次试验的平均值）。
Rclone 支持类似于 s5cmd 的 beast mode 的功能，通过 files-from 命令行选项实现。下面我们使用 transfers 值设置为 256 在我们的 800 个文件上运行 rclone copy，以匹配 s5cmd 的默认 并发设置。
rclone -P --transfers 256 --files-from files.txt S3store:my-bucket /my-local
files.txt 文件包含 800 行，格式如下：
small_files/<i>.jpg
我们的 800 个文件的 rclone copy 平均耗时 8.5 秒，比 s5cmd 略少（十次试验的平均值）。
我们承认目前展示的结果可能不足以说服您选择 rclone 而非现有工具。在下一节中，我们将描述一个用例，突出 rclone 的潜在优势。
对象存储系统之间的数据传输
现如今，开发团队维护多个对象存储并不罕见。这可能是为了防范存储故障的风险，或是决定使用多个云服务提供商的数据处理服务。例如，您可能会在 AWS 中使用存储在 Amazon S3 中的数据进行 AI 模型训练，同时在 Microsoft Azure 中运行数据分析，分析的数据存储在 Azure Storage 中。此外，您可能还希望在 FlashBlade、Cloudian 或 VAST 等本地存储基础设施中保持数据的备份。这些情况需要在安全、可靠和及时的方式下，在多个对象存储之间传输和同步数据的能力。
一些云服务提供商提供专门用于此类目的的服务。然而，这些服务并不总是满足您项目的具体需求，或者可能无法提供您所期望的控制级别。例如，Google Storage Transfer 在指定存储文件夹内的 所有数据 迁移方面表现出色，但（截至本文撰写时）不支持从其中传输特定子集的文件。
我们可以考虑的另一个选项是将现有的数据管理方案应用于这个目的。问题在于，像 AWS CLI 和 s5cmd 这样的工具（截至目前）不支持为源存储系统和目标存储系统指定不同的访问设置和安全凭证。因此，在存储位置之间迁移数据需要将数据转移到一个中间（临时）位置。下面的命令中，我们结合使用了 s5cmd 和 AWS CLI，通过系统内存和 Linux 管道将文件从 Amazon S3 复制到 Google Storage：
s5cmd cat s3://my-bucket/file \
| aws s3 cp --endpoint-url https://storage.googleapis.com
--profile gcp - s3://gs-bucket/file
尽管这是一种合法的、尽管笨拙的传输单个文件的方式，但在实际操作中，我们可能需要能够传输数百万个文件。为了支持这一点，我们需要添加一个额外的层来启动和管理多个并行的工作进程/处理器。事情可能会迅速变得棘手。
使用 Rclone 进行数据传输
与 AWS CLI 和 s5cmd 等工具不同，rclone 使我们能够为源和目标指定不同的访问设置。在以下的 rclone 配置文件中，我们添加了 Google Cloud Storage (GCS) 访问设置。（在这里，我们将 GCS 视为 S3 提供商。请查看这里了解其他配置选项。）
[S3store]
type = s3
provider = AWS
access_key_id = <id>
secret_access_key = <key>
[GSstore]
type = s3
provider = GCS
access_key_id = <id>
secret_access_key = <key>
endpoint = https://storage.googleapis.com
在存储系统之间传输单个文件的格式与复制到本地目录相同：
rclone copy -P S3store:my-bucket/file GSstore:gs-bucket/file
然而，rclone 的真正力量来自于将上述的files-from选项与这个功能结合。我们无需为数据迁移的并行化设计自定义解决方案，只需用一个命令就能传输一长串文件：
rclone copy -P --transfers 256 --files-from files.txt \
S3store:my-bucket/file GSstore:gs-bucket/file
实际上，我们可以通过将对象文件列表解析为较小的列表（例如，每个列表包含 10,000 个文件）并在单独的计算资源上运行每个列表，从而进一步加速数据迁移。虽然这种解决方案的具体影响因项目而异，但它可以显著提高开发的速度和效率。
总结
在这篇文章中，我们探讨了使用 rclone 进行基于云的存储管理，并展示了它在维护和同步多个存储系统数据方面的应用。数据传输确实有许多替代解决方案，但 rclone 基于的方法的便利性和优雅性无可置疑。
这只是我们在最大化基于云的存储解决方案效率方面撰写的众多文章之一。请务必查看一些我们的其他文章以了解这一重要话题。
在编写 Apache Beam 管道时使用示例进行 Map、Filter 和 CombinePerKey 转换
让我们用一些真实数据进行练习
·发表于 Towards Data Science ·8 分钟阅读·2023 年 7 月 12 日
–
Apache Beam 作为统一的编程模型，在高效和可移植的大数据处理管道中越来越受欢迎。它可以处理批量数据和流数据。这也是名字的由来。Beam 是 Batch 和 Stream 两个词的组合：
B（来自Batch）+ eam（来自 stream）= Beam
便携性也是一个很棒的特性。你只需专注于运行管道，它可以在任何地方运行，例如 Spark、Flink、Apex 或 Cloud Dataflow。你不需要更改逻辑或语法。
在这篇文章中，我们将专注于学习如何通过示例编写一些 ETL 管道。我们将尝试使用一个好的数据集进行一些转换操作，希望你会发现这些转换操作在工作中也非常有用。
请随意下载这个公共数据集并跟随练习：
这个练习使用了 Google Colab notebook。因此，安装非常简单。只需使用这一行代码：
!pip install --quiet apache_beam
安装完成后，我为这个练习创建了一个名为‘data’的目录：
mkdir -p data
让我们深入探讨今天的话题——转换操作。首先，我们将处理一个最简单的管道，即读取 CSV 文件并将其写入文本文件。
这不像 Padas 的 read_csv() 方法那么简单。它需要一个 coder() 操作。首先，在这里定义了一个 CustomCoder() 类，该类首先将对象编码为字节字符串，然后将字节解码为其对应的对象，并最终指定这个 coder 是否保证对值进行确定性编码。请查看 文档这里。
如果这是你的第一个管道，请注意管道的语法。在 CustomCoder() 类之后是最简单的管道。我们首先将空管道初始化为‘p1’。然后我们编写了‘sales’管道，其中首先从我们之前创建的数据文件夹中读取 CSV 文件。在 Apache beam 中，管道中的每个转换操作都以 | 符号开始。读取 CSV 文件中的数据后，我们只是将其写入文本文件。最后，通过 run() 方法我们运行了管道。这是 Apache beam 中标准和常用的管道语法。
import apache_beam as beam
from apache_beam.coders.coders import Coder
class CustomCoder(Coder):
"""A custom coder used for reading and writing strings as UTF-8."""
def encode(self, value):
return value.encode("utf-8", "replace")
def decode(self, value):
return value.decode("utf-8", "ignore")
def is_deterministic(self):
return True
p1 = beam.Pipeline()
sales = (p1
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.io.WriteToText('data/output'))
p1.run()
如果你现在检查你的‘data’文件夹，你会看到一个‘output-00000-of-00001’文件。从这个文件中打印前 5 行以检查数据：
!head -n 5 data/output-00000-of-00001
输出：
10107,30,95.7,2,2871,2/24/2003 0:00,Shipped,1,2,2003,Motorcycles,95,S10_1678,Land of Toys Inc.,2125557818,897 Long Airport Avenue,,NYC,NY,10022,USA,NA,Yu,Kwai,Small
10121,34,81.35,5,2765.9,5/7/2003 0:00,Shipped,2,5,2003,Motorcycles,95,S10_1678,Reims Collectables,26.47.1555,59 rue de l'Abbaye,,Reims,,51100,France,EMEA,Henriot,Paul,Small
10134,41,94.74,2,3884.34,7/1/2003 0:00,Shipped,3,7,2003,Motorcycles,95,S10_1678,Lyon Souveniers,+33 1 46 62 7555,27 rue du Colonel Pierre Avia,,Paris,,75508,France,EMEA,Da Cunha,Daniel,Medium
10145,45,83.26,6,3746.7,8/25/2003 0:00,Shipped,3,8,2003,Motorcycles,95,S10_1678,Toys4GrownUps.com,6265557265,78934 Hillside Dr.,,Pasadena,CA,90003,USA,NA,Young,Julie,Medium
10159,49,100,14,5205.27,10/10/2003 0:00,Shipped,4,10,2003,Motorcycles,95,S10_1678,Corporate Gift Ideas Co.,6505551386,7734 Strong St.,,San Francisco,CA,,USA,NA,Brown,Julie,Medium
Map
让我们来看一下如何在上述管道中使用 Map 转换。这是最常见的转换操作。你在 Map 中指定的转换将应用于 PCollection 中的每一个元素。
例如，我想添加一个 split 方法以从 PCollection 中的每个元素创建列表。在这里，我们将使用 lambda 进行 Map 转换。如果你不熟悉 lambda，请查看这里的 lambda 代码。lambda 后我们提到‘row’，任何其他变量名也可以。我们对‘row’应用的任何函数或方法将应用于 PCollection 中的每个元素。
p2 = beam.Pipeline()
sales = (p2
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.Map(lambda row: row.split(','))
|beam.io.WriteToText('data/output2'))
p2.run()
看，它是完全相同的语法。只是我在读取和写入操作之间多了一行代码。再次打印输出的前 5 行进行检查：
!head -n 5 data/output2-00000-of-00001
输出：
['10107', '30', '95.7', '2', '2871', '2/24/2003 0:00', 'Shipped', '1', '2', '2003', 'Motorcycles', '95', 'S10_1678', 'Land of Toys Inc.', '2125557818', '897 Long Airport Avenue', '', 'NYC', 'NY', '10022', 'USA', 'NA', 'Yu', 'Kwai', 'Small']
['10121', '34', '81.35', '5', '2765.9', '5/7/2003 0:00', 'Shipped', '2', '5', '2003', 'Motorcycles', '95', 'S10_1678', 'Reims Collectables', '26.47.1555', "59 rue de l'Abbaye", '', 'Reims', '', '51100', 'France', 'EMEA', 'Henriot', 'Paul', 'Small']
['10134', '41', '94.74', '2', '3884.34', '7/1/2003 0:00', 'Shipped', '3', '7', '2003', 'Motorcycles', '95', 'S10_1678', 'Lyon Souveniers', '+33 1 46 62 7555', '27 rue du Colonel Pierre Avia', '', 'Paris', '', '75508', 'France', 'EMEA', 'Da Cunha', 'Daniel', 'Medium']
['10145', '45', '83.26', '6', '3746.7', '8/25/2003 0:00', 'Shipped', '3', '8', '2003', 'Motorcycles', '95', 'S10_1678', 'Toys4GrownUps.com', '6265557265', '78934 Hillside Dr.', '', 'Pasadena', 'CA', '90003', 'USA', 'NA', 'Young', 'Julie', 'Medium']
['10159', '49', '100', '14', '5205.27', '10/10/2003 0:00', 'Shipped', '4', '10', '2003', 'Motorcycles', '95', 'S10_1678', 'Corporate Gift Ideas Co.', '6505551386', '7734 Strong St.', '', 'San Francisco', 'CA', '', 'USA', 'NA', 'Brown', 'Julie', 'Medium']
看，每个元素都变成了一个列表。
Filter
接下来，我将把 Filter 转换也添加到上述代码块中。这里的 lambda 也可以用于过滤。我们将过滤掉所有数据，只保留 Produc line 中的‘经典汽车’数据。数据集的第 11 列是产品线。正如你所知，Python 是零索引的。所以，列号的计数也从零开始。
p3 = beam.Pipeline()
sales = (p3
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.Map(lambda row: row.split(','))
|beam.Filter(lambda row:row[10] == "Classic Cars")
|beam.io.WriteToText('data/output3'))
p3.run()
如之前所示，打印前 5 行进行检查：
!head -n 5 data/output3-00000-of-00001
输出：
['10103', '26', '100', '11', '5404.62', '1/29/2003 0:00', 'Shipped', '1', '1', '2003', 'Classic Cars', '214', 'S10_1949', 'Baane Mini Imports', '07-98 9555', 'Erling Skakkes gate 78', '', 'Stavern', '', '4110', 'Norway', 'EMEA', 'Bergulfsen', 'Jonas', 'Medium']
['10112', '29', '100', '1', '7209.11', '3/24/2003 0:00', 'Shipped', '1', '3', '2003', 'Classic Cars', '214', 'S10_1949', '"Volvo Model Replicas', ' Co"', '0921-12 3555', 'Berguvsvgen 8', '', 'Lule', '', 'S-958 22', 'Sweden', 'EMEA', 'Berglund', 'Christina', 'Large']
['10126', '38', '100', '11', '7329.06', '5/28/2003 0:00', 'Shipped', '2', '5', '2003', 'Classic Cars', '214', 'S10_1949', '"Corrida Auto Replicas', ' Ltd"', '(91) 555 22 82', '"C/ Araquil', ' 67"', '', 'Madrid', '', '28023', 'Spain', 'EMEA', 'Sommer', 'Martn', 'Large']
['10140', '37', '100', '11', '7374.1', '7/24/2003 0:00', 'Shipped', '3', '7', '2003', 'Classic Cars', '214', 'S10_1949', 'Technics Stores Inc.', '6505556809', '9408 Furth Circle', '', 'Burlingame', 'CA', '94217', 'USA', 'NA', 'Hirano', 'Juri', 'Large']
['10150', '45', '100', '8', '10993.5', '9/19/2003 0:00', 'Shipped', '3', '9', '2003', 'Classic Cars', '214', 'S10_1949', '"Dragon Souveniers', ' Ltd."', '+65 221 7555', '"Bronz Sok.', ' Bronz Apt. 3/6 Tesvikiye"', '', 'Singapore', '', '79903', 'Singapore', 'Japan', 'Natividad', 'Eric', 'Large']
查看上面输出中每个列表的第 11 个元素。它是‘经典汽车’。
回答一些问题
每种类型的汽车订购了多少数量？
为了找出这一点，我们首先将创建元组，其中第一个元素或键将来自数据集的第 11 个元素，第二个元素即值将是数据集的第二个元素，即‘订购数量’。在下一步中，我们将使用 CombinePerKey() 方法。正如名字所示，它将为每个键结合具有聚合函数的值。
当你看到代码时会更清楚。这里是代码。
p3a = beam.Pipeline()
sales = (p3a
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.Map(lambda row: row.split(','))
#|beam.Filter(lambda row:row[10] == "Classic Cars")
|beam.Map(lambda row: (row[10], int(row[1])))
|beam.io.WriteToText('data/output3a'))
p3a.run()
!head -n 10 data/output3a-00000-of-00001
如你所见，我们在这里使用了两次 Map 函数。第一次是分割并像之前一样生成列表，然后从每行数据中提取第 10 列的产品线和第二列的数量。
这里是输出：
('Motorcycles', 30)
('Motorcycles', 34)
('Motorcycles', 41)
('Motorcycles', 45)
('Motorcycles', 49)
('Motorcycles', 36)
('Motorcycles', 29)
('Motorcycles', 48)
('Motorcycles', 22)
('Motorcycles', 41)
只是打印了输出的前 10 行。如你所见，这里每行数据的订单数量都列出了。回答上述问题的下一步也是最后一步是将每项的所有值结合起来。Apache Beam 中有 CombinePerKey 方法可以实现。顾名思义，它会为每个键使用聚合函数来结合值。在这种情况下，我们需要的是“总和”。
p4 = beam.Pipeline()
sales = (p4
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.Map(lambda row: row.split(','))
#|beam.Filter(lambda row:row[10] == "Classic Cars")
|beam.Map(lambda row: (row[10], int(row[1])))
|beam.CombinePerKey(sum)
|beam.io.WriteToText('data/output4'))
p4.run()
!head -n 10 data/output4-00000-of-00001
输出：
('Motorcycles', 11663)
('Classic Cars', 33992)
('Trucks and Buses', 10777)
('Vintage Cars', 21069)
('Planes', 10727)
('Ships', 8127)
('Trains', 2712)
所以，我们得到了每个产品的总订单数量。
哪些州的订单数量超过了 2000 个？
这是一个有趣的问题，我们需要每一个之前做过的变换加上另一个过滤变换。我们需要计算每个州的总订单数量，就像在前面的例子中计算每个产品的总订单数量一样。然后，将数量超过 2000 的订单进行过滤。
在所有之前的例子中，lambda 函数被用于 Map 和 Filter 变换。这里我们将看到如何定义一个函数并在 Map 或 Filter 函数中使用它。这里定义了一个函数 quantity_filter()，它返回值数量大于 2000 的项。
def quantity_filter(row):
name, count = row
if count > 2000:
return row
p7 = beam.Pipeline()
sales = (p7
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.Map(lambda row: row.split(','))
|beam.Map(lambda row: (row[17], int(row[1])))
|beam.CombinePerKey(sum)
|beam.Map(quantity_filter)
|beam.io.WriteToText('data/output7'))
p7.run()
!head -n 10 data/output7-00000-of-00001
输出：
('NYC', 5294)
None
None
None
('San Francisco', 2139)
None
('', 33574)
None
None
None
这是输出，其中如果数量不超过 2000，则返回‘None’。我不喜欢看到所有这些‘None’值。我将添加另一个过滤变换来过滤掉这些‘None’值。
p8 = beam.Pipeline()
sales = (p8
|beam.io.ReadFromText('data/sales_data_sample.csv', coder=CustomCoder(), skip_header_lines=1)
|beam.Map(lambda row: row.split(','))
|beam.Map(lambda row: (row[17], int(row[1])))
|beam.CombinePerKey(sum)
|beam.Map(quantity_filter)
|beam.Filter(lambda row: row != None)
|beam.io.WriteToText('data/output8'))
p8.run()
!head -n 10 data/output8-00000-of-00001
输出：
('NYC', 5294)
('San Francisco', 2139)
('', 33574)
('New Bedford', 2043)
('San Rafael', 6366)
所以，我们总共有 5 个返回的值，其中订单数量大于 2000。
结论
在本教程中，我想展示如何在 Apache Beam 中使用 Map、Filter 和 CombinePerKey 变换来编写 ETL 管道。希望它们足够清晰，能够在你的项目中使用。我将在下一篇文章中解释如何使用 ParDo。
随意关注我 Twitter 并点赞我的 Facebook 页面。
相关阅读
关于 Python 中多项式回归的详细教程、概述、实现和过拟合 | 作者：Rashida Nasrin Sucky | 2023 年 6 月 | Towards AI
迷你 VGG 网络图像识别的完整实现 | 作者：Rashida Nasrin Sucky | Towards Data Science
如何在 TensorFlow 中定义自定义层、激活函数和损失函数 | 作者：Rashida Nasrin Sucky | Towards Data Science
在 TensorFlow 中开发多输出模型的逐步教程 | 作者：拉希达·纳斯林·苏基 | 数据科学前沿
OpenCV Python 中的简单边缘检测方法 | 作者：拉希达·纳斯林·苏基 | 数据科学前沿
轨迹预测中的地图匹配
你要去哪里？你是否应该朝那个方向前进？
·
关注 发布于 Towards Data Science ·16 min read·2023 年 7 月 20 日
–
照片由Mateusz Wacławek拍摄，来源于Unsplash
本文提出了一种方法，利用从嘈杂的 GPS 传感器采样的历史旅行数据库来预测数字道路网络上的车辆轨迹。除了预测未来方向外，该方法还为任意位置序列分配一个概率。
这一理念的核心是使用数字地图，我们将所有采样位置投影到地图上，通过将它们聚合到单独的轨迹中并与地图匹配。这一匹配过程将连续的 GPS 空间离散化为预定位置和序列。将这些位置编码为唯一的地理空间标记后，我们可以更容易地预测序列，评估当前观察的概率，并估计未来的方向。这是本文的要旨。
问题
我在这里尝试解决哪些问题？如果你需要分析车辆路径数据，可能需要回答文章小标题中的问题。
你要去哪里？你应该那样走吗？
如何评估观察到的路径是否遵循经常行驶的方向的概率？这是一个重要的问题，因为通过回答它，你可以编程一个自动化系统来根据观察到的频率对行程进行分类。新轨迹的低分将引起关注并促使立即标记。
如何预测车辆接下来会进行哪些操作？它会继续直行，还是在下一个交叉口右转？你期望在接下来的十分钟或十英里内看到车辆在哪里？对这些问题的快速回答将帮助在线跟踪软件解决方案向配送规划者、在线路线优化器，甚至机会充电系统提供答案和见解。
解决方案
我在这里提出的解决方案使用了一个历史轨迹数据库，每个轨迹由特定车辆运动生成的时间序列位置组成。每个位置记录必须包含时间、位置信息、车辆标识符的参考以及轨迹标识符。一辆车有很多轨迹，每个轨迹有很多位置记录。我们输入数据的样本如图 1所示。
图 1 — 上表显示了来自扩展车辆能量数据集的一个轨迹的小样本。虽然每一行包含比展示的更多属性，但我们只需要隐含顺序和位置。注意，由于采样策略，有很多重复的位置。我们稍后将处理这个问题。（图片来源：作者）
我从扩展车辆能量数据集（EVED）[1] 文章中提取了上述数据。你可以通过参考我之前文章中的代码来构建相应的数据库。
本文解释了如何使用通过 Quadkeys 索引的已知速度向量来估计旅行时间。
towardsdatascience.com
我们的第一项工作是将这些轨迹匹配到支持的数字地图上。这个步骤的目的不仅是为了消除 GPS 数据采样误差，更重要的是，将获取的旅行数据强制转换为已知和固定的现有道路网络，其中每个节点和边都已知。每条记录的轨迹因此被转换为与现有数字地图节点相符的一系列数值标记。在这里，我们将使用开源数据和软件，地图数据来源于OpenStreetMap（由Geofabrik编译）、Valhalla地图匹配包和H3作为地理空间分词器。
边缘与节点匹配
地图匹配比乍看起来要复杂。为了说明这个概念的含义，我们来看一下下面的图 2。
图 2 — 上图展示了从 EVED 中获得的噪声轨迹（蓝色）。如你所见，它与最近的道路不匹配，需要与地图进行匹配。一旦我们将蓝线的顶点投影到地图上，就可以得到原始点到推断地图边缘的一系列投影，你可以看到结果轨迹为红色。然而，这条路径在某些地方仍然未能覆盖底层地图，特别是在图像的中心，红线在道路之间跳跃。我们的目标是重建地图上的旅行路径，如绿色线所示，遵循底层地图节点。（图片来源：作者使用 Folium 和 OpenStreetMap 图像）
上述图 2显示，我们可以从原始 GPS 序列中推导出两条轨迹。我们通过将原始 GPS 位置投影到最近（且最可能的）道路网络段来获得第一条轨迹。如你所见，结果折线有时只会沿着道路，因为地图使用图形节点来定义其基本形状。通过将原始位置投影到地图边缘，我们获得了属于地图的新点，但在通过直线连接到下一个点时，可能会偏离地图的几何形状。
通过将 GPS 轨迹投影到地图上的节点，我们得到了一条完美覆盖地图的路径，如图 2中的绿色线所示。虽然这条路径更好地表示了最初驱动的轨迹，但它不一定与原始位置一一对应。幸运的是，这对我们来说没有问题，因为我们将始终将任何轨迹与地图节点进行匹配，因此我们将始终获得一致的数据，只有一个例外。Valhalla 地图匹配代码始终将初始和最终轨迹点进行边缘投影，因此我们会系统性地丢弃它们，因为它们与地图节点不对应。
H3 分词
不幸的是，Valhalla 不报告唯一的道路网络节点标识符，因此我们必须将节点坐标转换为唯一的整数令牌，以便后续的序列频率计算。这就是 H3 介入的地方，它允许我们将节点坐标编码为一个独特的 64 位整数。我们选择 Valhalla 生成的多段线，去掉初始和最终点（这些点不是节点，而是边的投影），并将所有剩余的坐标映射到 level 15 H3 indices。
对偶图
使用上述过程，我们将每个历史轨迹转换为一系列 H3 令牌。下一步是将每个轨迹转换为令牌三元组序列。序列中的三个值表示预测图的两个连续边，我们希望知道这些的频率，因为它们将是预测和概率评估的核心数据。图 3 下面直观地描绘了这一过程。
图 3 — 左侧的地理空间令牌列表扩展为另一列表的三元组，表示隐含图的对偶视图。每个令牌是地理空间图上的一个节点，其序列表示边。转换后的列表将每个边视为对偶图中的一个节点，而中间的令牌是新的边，如右列所示。（图片来源：作者）
上述转换计算了道路图的对偶图，颠倒了原始节点和边的角色。
现在我们可以开始回答提出的问题。
你确定要走那条路吗？
要回答这个问题，我们需要知道车辆轨迹直到某个时刻。我们使用上述相同的过程进行映射匹配和令牌化轨迹，然后使用已知的历史频率计算每个轨迹三元组的频率。最终结果是所有个体频率的乘积。如果输入轨迹中有未知的三元组，它的频率将是零，作为最终路径概率。
三元组概率是特定序列 (A, B, C) 的计数与所有 *(A, B, ) 三元组计数的比率，如下图 图 4 所示。
图 4 — 三元组概率是其频率与所有具有相同两个初始令牌的三元组频率的比率。（图片来源：作者）
旅行概率只是各个旅行三元组的乘积，如下图图 5所示。
图 5 — 旅行概率是所有匹配三元组的简单乘积。（图片来源：作者）
你要去哪里？
我们使用相同的原则来回答这个问题，但仅从最后一个已知的三元组开始。我们可以使用这个三元组作为输入，通过列举所有以输入的最后两个令牌作为前两个令牌的三元组，预测最可能的 k 个后继者。下方的图 6展示了三元组序列生成和评估的过程。
图 6 — 在这个虚拟的案例中，最可能的下一个三元组是观察到的频率最高的三元组（B, C, D）。(图片来源：作者)
我们可以提取前 k 个后继三元组并重复该过程，以预测最可能的旅行。
实现
我们准备讨论实现细节，从地图匹配和一些相关概念开始。接下来，我们将学习如何从 Python 使用 Valhalla 工具集，提取匹配的路径并生成令牌序列。一旦我们将结果存储在数据库中，数据预处理步骤就完成了。
最后，我展示了一个使用 Streamlit 的简单用户界面，该界面计算任何手绘轨迹的概率，并将其投射到未来。
地图匹配
地图匹配 将从移动物体路径中采样的 GPS 坐标转换为现有的道路图。道路图是一个离散的模型，表示物理道路网络，由 节点 和连接的 边 组成。每个节点对应于沿道路已知的地理位置，编码为纬度、经度和高度元组。每个 有向边 连接沿着基础道路的相邻节点，并包含许多属性，如航向、最高速度、道路类型等。下方的图 7用一个简单的例子说明了这个概念。
图 7 — 上图展示了一个小型数字道路网络，突出了一个交叉口。每个红点代表现有道路上的已知地理位置。蓝色线条表示节点之间的连接边。请注意，这些边通常是有向的，并且可能是多个的。(图片来源：作者)
成功时，地图匹配过程会生成有关采样轨迹的相关和有价值的信息。一方面，该过程将采样的 GPS 点投射到最可能的道路图 边 上。地图匹配过程通过将观察到的点准确地放置到推断的道路图 边 上来“纠正”观察到的点。另一方面，该方法还通过提供最可能的路径，通过道路图重建图 节点 的序列，以对应于采样的 GPS 位置。请注意，如前所述，这些输出是不同的。第一个输出包含沿着最可能路径的 边 的坐标，而第二个输出由重建的图 节点 序列组成。下方的图 8展示了这个过程。
图 8 — 上图展示了地图匹配过程，其中绿色点代表观察到的 GPS 坐标，橙色钻石代表沿已知边的投影位置。请注意，对于上述简化的示例，我们只能安全地重建节点 2 和 3 之间的路径。这种困境并不像看起来那么严重，因为在实际地图中，轨迹匹配的边缘远远超过一个，因此信息丢失最小。（图片来源：作者）
地图匹配过程的一个副产品是使用共享道路网络表示法对输入位置进行标准化，特别是在考虑第二种输出类型时：最可能的节点序列。在将采样的 GPS 轨迹转换为一系列节点时，我们通过将推断路径减少为一系列节点标识符，使其可比拟。我们可以将这些节点序列视为已知语言的短语，其中每个推断的节点标识符是一个词，其排列传达了行为信息。
这是第五篇文章，我在其中探讨了扩展车辆能源数据集¹ (EVED) [1]。该数据集是对先前工作的扩展和评审，并提供了原始 GPS 采样位置的地图匹配版本（上方图 8中的橙色钻石）。
不幸的是，EVED 仅包含投影的 GPS 位置，而缺少重建的道路网络节点序列。在我之前的两篇文章中，我解决了从转换后的 GPS 位置重建道路段序列而无需地图匹配的问题。我发现结果有些令人失望，因为我预期的缺陷重建率低于观察到的 16%。您可以从以下文章中跟踪这一讨论。
三角形在地理空间查询中具有强大的特性
towardsdatascience.com ## 更多关于道路网络匹配
道路网络匹配的趣事
towardsdatascience.com
现在，我正在查看源地图匹配工具，以了解它在纠正缺陷重建方面的能力。因此，让我们对Valhalla进行测试。以下是我用来在 Docker 容器中运行 Valhalla 的步骤、参考文献和代码。
Valhalla 设置
在这里，我紧跟桑迪普·潘迪 [2]在其博客上提供的说明。
首先，确保你的机器上已安装 Docker。要安装 Docker 引擎，请参阅 在线说明。如果你使用的是 Mac，另一个很好的选择是 Colima。
安装完成后，你必须从 GitHub 拉取一个 Valhalla 镜像，方法是按照下面的 Figure 9 所示在命令行中发出以下命令。
Figure 9 — 从命令行拉取 Valhalla 的 Docker 镜像。（图片来源：作者）
在执行上述命令时，你可能需要输入你的 GitHub 凭据。此外，确保你已经克隆了本文的 GitHub 仓库，因为一些文件和文件夹结构会引用它。
完成后，你应该打开一个新的终端窗口，并发出以下命令以启动 Valhalla API 服务器（MacOS、Linux、WSL）：
Figure 10 — 上述命令在 Docker 容器中运行拉取的 Valhalla 镜像。首次执行时，该命令还会在启动之前下载和准备最新的 Geofabrik Michigan 数据文件。（图片来源：作者）
上述命令行明确指定了要从 Geofabrik 服务下载哪个 OSM 文件，即最新的 Michigan 文件。这一指定意味着第一次执行时，服务器将下载并处理该文件并生成优化后的数据库。在后续调用中，服务器将省略这些步骤。需要时，删除目标目录下的所有内容以刷新下载的数据，并重新启动 Docker。
我们现在可以使用专用客户端调用 Valhalla API。
输入 PyValhalla
这个衍生项目简单地提供了对精彩的 Valhalla 项目 的打包 Python 绑定。
使用 PyValhalla Python 包非常简单。我们从使用以下命令行进行的简洁安装过程开始。
Figure 11 — 你可以使用 PIP 安装 PyValhalla。（图片来源：作者）
在你的 Python 代码中，你必须导入所需的引用，从处理后的 GeoFabrik 文件中实例化配置，最后创建一个 Actor 对象，这是你访问 Valhalla API 的入口。
Figure 12 — 上述代码展示了如何在 Python 应用程序或笔记本上轻松设置 PyValhalla。（图片来源：作者）
在我们调用 Meili 地图匹配服务之前，我们必须使用 Figure 13 中列出的函数获取轨迹 GPS 位置。
Figure 13 — 上述函数加载车辆轨迹的唯一位置，返回一个包含纬度、经度和时间戳的 Pandas DataFrame。（图片来源：作者）
我们现在可以设置参数字典以传递给 PyValhalla 调用以追踪路线。有关这些参数的更多细节，请参见 Valhalla 文档。下面的函数调用了 Valhalla（Meili）中的地图匹配功能，并包含在 数据准备脚本 中。它展示了如何从包含观测 GPS 位置的 Pandas 数据框中推断路线，这些位置编码为纬度、经度和时间元组。
图 14 — 上述函数接受一个 PyValhalla Actor 对象和一个包含源路径的 Pandas DataFrame，并返回一个地图匹配的字符串编码折线。这个字符串随后会解码为与数字地图网络节点对应的地理空间位置列表，极端位置除外，这些位置被投影到边缘上。（图片来源：作者）
上述函数返回的匹配路径是字符串编码的折线。如下面的数据准备代码所示，我们可以使用 PyValhalla 库调用轻松解码返回的字符串。请注意，这个函数返回的是一条折线，其第一个和最后一个位置被投影到边缘，而不是图节点。您将看到这些极端位置在本文后面的代码中被去除。
现在让我们来看看数据准备阶段，我们将 EVED 数据库中的所有轨迹转换为一组地图边缘序列，从中我们可以导出模式频率。
数据准备
数据准备的目的是将噪声较大的 GPS 获取轨迹转换为对应于已知地图位置的地理空间标记序列。主要代码遍历现有的行程，一次处理一个。
我在本文中使用 SQLite 数据库来存储所有数据处理结果。我们从填充匹配的轨迹路径开始。您可以参考下面的 图 15 中的代码描述。
图 15 — 上述代码包含了预处理数据的循环。这个循环遍历已知的轨迹，计算它们的地图匹配路径（如果有的话），将节点分词，并将其扩展为三元组。代码将所有中间结果和最终结果存储在数据库中。（图片来源：作者）
对于每条轨迹，我们实例化一个Actor类型的对象（第 9 行）。这是一个未明确说明的要求，因为每次调用地图匹配服务都需要一个新实例。接下来，我们加载（第 13 行）车辆 GPS 接收器获取的轨迹点，这些点带有原始 VED 文章中提到的添加噪声。在第 14 行，我们调用 Valhalla 进行地图匹配，检索编码后的匹配路径并将其保存到数据库。接着，我们将编码后的字符串解码成一组地理空间坐标列表，去除两端极值（第 17 行），然后将其转换为计算在 15 级别 H3 网格上的 H3 索引列表（第 19 行）。在第 23 行，我们将转换后的 H3 索引和原始坐标保存到数据库，以便后续的反向映射。最后，在第 25至27行，我们基于 H3 索引列表生成一系列 3 元组，并保存它们以便后续推断计算。
让我们逐步分析每一个步骤，并详细解释它们。
轨迹加载
我们已经看到如何从数据库中加载每条轨迹（参见图 13）。一条轨迹是一个按时间顺序排列的采样 GPS 位置序列，编码为纬度和经度对。请注意，我们没有使用 EVED 数据提供的匹配版本的这些位置。在这里，我们使用最初 VED 数据库中存在的带有噪声和原始坐标。
地图匹配
调用地图匹配服务的代码已经在上文的图 14中介绍过。其核心问题在于配置设置；除此之外，这是一个非常直接的调用。将结果编码后的字符串保存到数据库中也很简单。
图 16 — 上述代码将编码后的折线字符串保存到新的数据库中。（图片来源：作者）
在主循环的第 15 图中的第 17 行，我们将几何字符串解码为纬度和经度元组列表。请注意，这是我们剥离初始和最终位置的地方，因为它们没有投影到节点上。接下来，在第 19 行，我们将此列表转换为相应的 H3 标记列表。我们使用最大详细级别来尝试避免重叠，并确保 H3 标记与地图图形节点之间的一对一关系。在接下来的两行中，我们将这些标记插入数据库中。首先，我们保存整个标记列表，并将其与轨迹关联。
图 17 — 上述函数将轨迹的 H3 标记列表插入数据库中。（图片来源：作者）
接下来，我们插入节点坐标到 H3 标记的映射，以便从给定的标记列表绘制折线。这一功能在推断未来行程方向时将会很有帮助。
图 18 — 我们插入 H3 标记和节点坐标之间的映射，以便从给定的推断标记重构轨迹。（图片来源：作者）
我们现在可以生成并保存相应的令牌三元组。下面的函数使用新生成的 H3 令牌列表并将其扩展为另一个三元组列表，如图 3中详细说明的那样。扩展代码在图 19中展示。
图 19 — 上述代码将 H3 令牌列表转换为相应三元组的列表。（图片来源：作者）
在三重扩展之后，我们可以将最终产品保存到数据库中，如下方的图 20所示。通过巧妙地查询这个表，我们将推断当前的三重概率和未来最可能的轨迹。
图 20 — 上述函数将 H3 三元组保存到数据库中。这是数据准备阶段的最后一步。我们现在可以开始探索我们收集的信息。（图片来源：作者）
我们已经完成了一轮数据准备循环。一旦外循环完成，我们将拥有一个新的数据库，所有轨迹都转换为令牌序列，我们可以随意探索。
你可以在 GitHub 仓库中找到完整的数据准备代码。
概率和预测
我们现在转向估计现有行程概率和预测未来方向的问题。让我们先定义一下“现有行程概率”的含义。
行程概率
我们从通过地图匹配投影到道路网络节点的任意路径开始。因此，我们有一个来自地图的节点序列，并希望评估该序列的概率，使用已知的行程数据库作为频率参考。我们使用图 5中的公式。简而言之，我们计算所有单独的三重组概率的乘积。
为了说明这一功能，我实现了一个简单的Streamlit 应用程序，允许用户在覆盖的安娜堡区域绘制任意行程并立即计算其概率。
一旦用户在地图上绘制表示行程或假设 GPS 样本的点，代码会将它们进行地图匹配以检索底层的 H3 令牌。从那时起，只需计算单个三重组的频率并将其相乘即可计算总概率。图 21中的函数计算任意行程的概率。
图 21 — 上述函数从三重频率数据库中计算任意路径的概率。（图片来源：作者）
该代码得到另一个函数的支持，该函数检索任何现有的 H3 令牌对的后继。下列函数在图 22中列出，查询频率数据库并返回一个 Python Counter 对象，包含输入令牌对所有后继的计数。当查询未找到后继时，函数返回None 常量。注意该函数如何使用缓存以提高数据库访问性能（此处未列出代码）。
图 22 — 上述函数查询频率数据库以获取任何一对 H3 令牌的已知后继，并返回一个包含所有后继计数的 Counter 对象。(图片来源：作者)
我设计了这两个函数，使得在任何给定节点没有已知后继时，计算的概率为零。
让我们看看如何预测轨迹的最可能未来路径。
预测方向
我们只需要从给定的运行轨迹中获取最后两个令牌，以预测其最可能的未来方向。这个想法涉及扩展该令牌对的所有后继，并选择最频繁的那些。下面的代码展示了作为方向预测服务入口点的函数。
图 23 — 上述函数从 Folium 中填充一个 FeatureGroup 对象，包含现有用户提供的轨迹的预测路径。(图片来源：作者)
上述函数首先通过检索用户绘制的轨迹作为一组与地图匹配的 H3 令牌，并提取最后一对。我们将这一对令牌称为种子，并将在代码中进一步扩展它。在第 9 行，我们调用种子扩展函数，该函数返回一个与输入扩展标准对应的折线列表：每次迭代的最大分支数和总迭代次数。
让我们通过查看下列代码了解种子扩展函数的工作原理，如图 24所示。
图 24 — 种子扩展函数使用 PredictedPath 类来管理每次迭代。有关该类的更多详细信息，请见下文。(图片来源：作者)
通过调用生成最佳后继路径的路径扩展函数，种子扩展函数迭代地扩展路径，从初始路径开始。路径扩展通过选择一条路径并生成最可能的扩展，如图 25所示。
图 25 — 上述路径扩展函数迭代当前路径的最频繁后继。它为每个最频繁的后继创建一条新路径，使用一个专门的函数（见下文）。(图片来源：作者)
该代码通过将后继节点附加到源路径上来生成新路径，如下图 26所示。
图 26 — 要生成“子”路径，我们只需将后继节点附加到现有路径上，如下所示。注意，代码在附加新节点之前创建了原始路径的副本。(图片来源：作者)
该代码使用一个专门的类来实现预测路径，如图 27所示。
图 27 — 上述类实现了一个具有概率排序支持的预测路径，基于种子令牌对进行创建，并生成地图折线。（图片来源：作者）
应用
现在可以在下面的图 28中查看结果 Streamlit 应用程序。
图 28 — Streamlit 应用程序展示了两个描述的功能。输入轨迹为蓝色，可以使用地图左侧的工具菜单绘制。一旦绘制完成，代码将计算其概率并在底部显示。三条红色轨迹是源轨迹可能演变的三个最可能的五十边预测。点击每条轨迹可以弹出计算出的概率。（图片来源：作者）
结论
在这篇文章中，我介绍了一种预测车辆在数字地图道路网络中未来轨迹的方法。利用历史轨迹数据库，该方法为任何行程分配一个概率，并预测近期最可能的方向。因此，这种方法可以检测到不太可能的或甚至是前所未见的新轨迹。
我们从感兴趣区域的大量车辆轨迹数据库开始。每条路径都是地理坐标（纬度和经度）的时间顺序序列，以及其他相关属性，如速度。我们通常从车载 GPS 接收器收集这些轨迹，并将其集中编入数据库。
GPS 样本由于信号测量过程中不可避免的误差而产生噪声。自然和人工障碍物，如城市峡谷，可能显著降低信号的接收精度并增加地理位置误差。幸运的是，实用的解决方案通过概率匹配将 GPS 样本与数字地图对齐来解决这个问题。这就是地图匹配的全部内容。
通过将嘈杂的 GPS 样本与已知数字地图进行匹配，我们不仅通过将每个实例投影到地图上最可能的道路段来纠正精度问题，还获得了一系列车辆最可能经过的现有地图定义的位置。这个最后的结果对我们的轨迹预测至关重要，因为它本质上将一组嘈杂的 GPS 坐标转换为数字地图中干净且已知的点集。这些数字标记是固定的，不会改变，通过将 GPS 样本序列投影到这些标记中，我们得到一串已知的令牌，稍后可以用于预测。
我们使用已知令牌序列频率来计算所有概率，这些序列代表任意轨迹及其未来演变。结果是几个 Python 脚本，一个用于数据准备，另一个用于使用 Streamlit 平台进行数据输入和可视化。
备注
参考文献
[1] 张三、Fatih、Abdulqadir、Schwarz、和马晓（2022）。扩展车辆能源数据集（eVED）：一个用于深度学习车辆行程能源消耗的增强型大规模数据集。arXiv。 doi.org/10.48550/arXiv.2203.08630
[2] Valhalla 的高效快速地图匹配 — Sandeep Pandey（ikespand.github.io）
[3] 使用 Valhalla 的 Meili 正确完成地图匹配 | by Serge Zotov | Towards Data Science
João Paulo Figueira 是 tb.lx by Daimler Truck 在葡萄牙里斯本的数据科学家。
使用 R 绘制南美洲地图：深入探讨地理可视化
导航数据集、地缘政治细节和编码挑战，描绘大陆的全貌
·
关注 发表在 Towards Data Science · 8 分钟阅读 · 2023 年 8 月 30 日
–
图片由 Alexander Schimmeck 提供，来源于 Unsplash
所以你是那种从小就喜欢地图和地理的数据科学家和业余 Medium 作者。你正在寻找一个适合你后续图表工作，特别是地图的好主题，当你意识到你所在的国家巴西的官方统计局发布了最新的普查数据时。为什么不呢？为什么不对比一下巴西与其南美洲邻国的情况呢？这可能是使用 R 和所有优秀包的简单任务。让我们来做吧。
在做出这个决定的瞬间，意识到这个简单的任务实际上是一个英雄的旅程，包含了发现最合适的数据集与形状文件、信息缺乏、形状文件互操作性、纬度和经度数学、地理概念中的文化差异，甚至地缘政治问题，如如何将法国海外领土的地图和数据正确地放在南美洲中。
接下来的段落解释了在世界地图的限定区域绘制人口信息的一些可能路径。下面描述的逐步过程可能对所有那些对国际比较有兴趣的地理可视化方法者有用，即使他们的目的只是比较非洲国家的水资源获取情况或北美的肥胖率。
帕查玛玛
让我们从整体图像开始：R 版的世界地图。请见下图和代码。
世界地图：图片由作者提供
library(readxl)
library(geobr)
library(tidyverse)
library(sf)
library(spData)
library(ggrepel)
library(colorspace)
data("world")
#mapa mundi
world %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
scale_fill_continuous_sequential(palette= "Heat 2" )+
theme_void() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões de habitantes", 10)
)
我使用{spData}包作为具有全球领土形状文件几何信息的数据框的参考。aes 函数使用人口信息填充形状。众所周知，中国和印度是世界上人口最多的国家，每个国家都有超过 10 亿人。热度颜色显示了与其他国家的对比。大多数顺序颜色较弱。我们几乎无法理解图片中的颜色渐变。如果你想要更好的颜色分布，对数是最佳选择。见下文。
对数刻度的世界地图。图片由作者提供
world %>%
ggplot() +
geom_sf(aes(fill=pop)) +
scale_fill_continuous_sequential(palette= "Heat 2", trans= "log2" )+
theme_void() +
theme(
panel.background = element_rect(fill="#0077be"),
legend.position = "none"
)
在代码中，你可以看到 scale_fill_continuous_sequential 函数中的对数变换。
在世界数据框结构中，有一个“Continent”列。因此，使用该列筛选数据以获取南美洲地图是显而易见的。请查看代码，紧接着是地图。
南美洲地图：第一版。图片由作者提供
world %>%
filter(continent == "South America") %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
scale_fill_continuous_sequential(palette= "Heat 2" )+
theme_void() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões de habitantes", 10)
)
如你所见，dplyr 过滤函数运作良好；这正是我们想要看到的地图。但这真的正确吗？
气候变化是一个巨大的问题，但海平面尚未上升到淹没曾经出现在南美洲北部的明显区域的程度。这到底发生了什么？让我们现在借助坐标绘制另一张地图，并命名多边形。
南美洲地图：第二版。作者图片
southamerica<-
world %>%
filter(continent=="South America")
southamerica$lon<- sf::st_coordinates(sf::st_centroid(southamerica$geom))[,1]
southamerica$lat<- sf::st_coordinates(sf::st_centroid(southamerica$geom))[,2]
southamerica %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
scale_fill_continuous_sequential(palette= "Heat 2" )+
theme_light() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões de habitantes", 10)
)+
geom_text_repel(aes(x=lon, y=lat, label= str_wrap(name_long,20)),
color = "black",
fontface = "bold",
size = 2.8)
theme_light 代替 theme_void 足以显示坐标。多边形命名花费了更多工作。我们必须计算每个多边形的质心，然后将这些信息用作 geom_text_repel 函数中的 x 和 y 坐标。
使用这个新地图版本和一些先前的知识，我们发现缺失的领土是法属圭亚那，它位于北纬 0º 和 10º 之间，西经 53º 和 55º 之间。我们的下一个任务是了解如何获取法属圭亚那的信息：多边形、人口以及一些坐标来填补我们的地图。
La Mer
我必须将法国从世界其他地方隔离开来，以理解 {spData} 包如何处理这个国家地图的数据。见下文结果。
法国地图。作者图片
france<-
world %>%
filter(iso_a2 == "FR")
france %>%
ggplot() +
geom_sf(aes(fill=pop)) +
scale_fill_continuous_sequential(palette= "Heat 2", trans= "log2" )+
theme_light() +
theme(
#panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População", 30)
)
法国有许多所谓的海外领土。{spData} 包的方法是仅表示主要领土，加上科西嘉岛（地中海中的一个岛屿）和法属圭亚那，它位于精确的坐标范围内，这个范围正好填补了我们南美洲最后一张地图中的空白。
我的下一个尝试是将包含法国几何数据的数据框添加到我的南美洲过滤器中，但我知道我还需要更多。见下文
南美洲 + 法国。作者图片
southamerica %>%
bind_rows(france) %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
scale_fill_continuous_sequential(palette= "Heat 2" )+
theme_light() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões de habitantes", 10)
)+
geom_text_repel(aes(x=lon, y=lat, label= str_wrap(name_long,20)),
color = "black",
fontface = "bold",
size = 2.8)
正如你在代码中看到的，我使用 bind_row 将南美洲领土与法国的 shapefile 结合起来。这样我们现在就有了良好定位的法属圭亚那。另一方面，地图上没有人口信息，而法国在殖民历史的反面作为南美洲的一部分。
换句话说，我想要的是这张地图。
法属圭亚那在南美洲地图上。作者图片
data_guiana<-
insee::get_idbank_list('TCRED-ESTIMATIONS-POPULATION') %>%
filter(str_detect(REF_AREA_label_fr,"Guyane")) %>%
filter(AGE == "00-") %>% #all ages
filter(SEXE == 0) %>% #men and women
pull(idbank) %>%
insee::get_insee_idbank() %>%
filter(TIME_PERIOD == "2023") %>%
select(TITLE_EN,OBS_VALUE) %>%
mutate(iso_a2 = "FR")
data_guiana <- janitor::clean_names(data_guiana)
southamerica %>%
bind_rows(france) %>%
left_join(data_guiana) %>%
mutate(pop=ifelse(iso_a2=="FR",obs_value,pop))%>%
mutate(lon= ifelse(iso_a2=="FR", france[[11]][[1]][[1]][[1]][1,1], lon),
lat= ifelse(iso_a2=="FR",france[[11]][[1]][[1]][[1]][1,2], lat)) %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
scale_fill_continuous_sequential(palette= "Heat 2" )+
geom_text_repel(aes(x=lon, y=lat, label= str_wrap(name_long,20)),
color = "black",
fontface = "bold",
size = 2.8)+
coord_sf(xlim = c(-82,-35), ylim=c(-60,15))+
theme_light() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões de habitantes", 10)
)
正如你所看到的，我使用了由法国官方统计办公室制作的 R 包 来获取圭亚那的人口。此外，我将地图限制在适当的坐标范围内，以查看南美洲。
Emolduram e aquarelam o meu Brasil
现在地图英雄终于解决了南美洲的问题，并与法国演奏了 pipes of peace，是时候回到巴西的数据和地图了。记住，我想将一些巴西的人口普查细节与巴拿马以南其他国家和地区进行比较。
数据普查可以在 R 包 或 API 地址上找到。我选择了使用 API 的更具挑战性的选项。另一次使用其他选项可能是个好主意。查看下面的代码和地图，我展示了巴西各州的人口与其他南美洲领土的对比。
南美洲 + 巴西各州。图片由作者提供
central_america<-
world %>%
filter(subregion == "Central America")
brasil<- geobr::read_country()
estados<- geobr::read_state()
#dados de população
ibge2022<-
get_municipalies_data()
estados<-
estados %>%
inner_join(
ibge2022 %>%
rename(abbrev_state = uf) %>%
summarise(.by=abbrev_state,
pop = sum(populacao_residente)
)
)
southamerica %>%
filter(iso_a2!="BR") %>%
bind_rows(france) %>%
left_join(data_guiana) %>%
mutate(pop=ifelse(iso_a2=="FR",obs_value,pop))%>%
mutate(lon= ifelse(iso_a2=="FR", france[[11]][[1]][[1]][[1]][1,1], lon),
lat= ifelse(iso_a2=="FR",france[[11]][[1]][[1]][[1]][1,2], lat)) %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
geom_sf(data=estados, aes(fill=pop/10⁶)) +
geom_sf(data=brasil,fill=NA, color="#00A859", lwd=1.2)+
geom_sf(data= central_america,fill= "#808080")+
scale_fill_continuous_sequential(palette= "Heat 2" )+
geom_text_repel(aes(x=lon, y=lat,
label= str_wrap(name_long,20)),
color = "black",
fontface = "bold",
size = 2.8)+
coord_sf(xlim = c(-82,-35), ylim=c(-60,15))+
theme_void() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões", 10)
)
我使用上述 API 编写了函数 get_municipalites_data。代码可以在我的 gist 上找到。还请注意提供用于绘制巴西及其子区域边界的两个函数：read_country 和 read_states。这些函数在 {geobr} 包中。
我使用了世界数据框中的另一个筛选器。在这种情况下，目的是显示中美洲次大陆的起始部分，并用灰色阴影绘制其地图。在这里，我们面临了一种文化差异，因为我们在巴西了解到美洲有三个次大陆：北美、中美洲和南美洲。对于数据集的作者来说，中美洲是北美的一个子区域。
现在是时候结束我的工作了。我想在地图上显示八个最人口稠密的领土的名称。即使在最后的冲刺阶段，也有一些代码技巧。
最人口密集的领土。图片由作者提供
estados$lon<- sf::st_coordinates(sf::st_centroid(estados$geom))[,1]
estados$lat<- sf::st_coordinates(sf::st_centroid(estados$geom))[,2]
most_populated<-
southamerica %>%
filter(iso_a2 !="BR") %>%
rename(name= name_long) %>%
as_tibble() %>%
select(name, pop, lat, lon) %>%
bind_rows(
estados %>%
rename(name= name_state) %>%
as_tibble() %>%
select(name, pop, lat, lon)
) %>%
slice_max(order_by = pop, n=8)
southamerica %>%
filter(iso_a2!="BR") %>%
bind_rows(france) %>%
left_join(data_guiana) %>%
mutate(pop=ifelse(iso_a2=="FR",obs_value,pop))%>%
mutate(lon= ifelse(iso_a2=="FR", france[[11]][[1]][[1]][[1]][1,1], lon),
lat= ifelse(iso_a2=="FR",france[[11]][[1]][[1]][[1]][1,2], lat)) %>%
ggplot() +
geom_sf(aes(fill=pop/10⁶)) +
geom_sf(data=estados, aes(fill=pop/10⁶)) +
geom_sf(data=brasil,fill=NA, color="#00A859", lwd=1.2)+
geom_sf(data= central_america,fill= "#808080")+
scale_fill_continuous_sequential(palette= "Heat 2" )+
geom_text_repel(data= most_populated,
aes(x=lon, y=lat,
label= str_c(str_wrap(name,10),": ",round(pop/10⁶,1))),
color = "black",
fontface = "bold",
size = 2.9)+
coord_sf(xlim = c(-82,-35), ylim=c(-60,15))+
theme_void() +
theme(
panel.background = element_rect(fill="#0077be")
) +
labs(
fill= str_wrap("População em milhões", 10)
)
三个巴西州位于南美洲最人口稠密的八个领土中。实际上，圣保罗是地图上第二个最有人口的地区，仅次于哥伦比亚。
现在，专注于代码，你可以看到我创建了一个新的数据框，通过结合两个不同的 sf 对象来建立这个排名。我选择了一部分列，并将类型从 sf 更改为 tibble，以便进行行绑定。
就这样。英雄完成了一个可能的路径，并留下了下一次旅程的足迹。现在轮到你了。记住所有可能通过地图表示得到显著改善的项目。根据上述操作步骤，收集有关人口、社会经济问题等的所有数据，只需选择一个变量来填充多边形。
代码和数据
完整代码可以在 gist 上找到。
所有巴西数据集被认为是公共领域的数据，因为这些数据是由联邦政府机构生产的，作为主动透明度在互联网上提供，并且受到巴西信息公开法的约束。
IBGE: 巴西人口普查数据
IPEA: 巴西 shapefiles
法国的数据可以从开放数据门户获取，并标注为开放许可，允许将信息用于商业目的。
映射全球自然再造林项目的潜力
原文：
towardsdatascience.com/mapping-the-global-potential-of-natural-reforestation-projects-c425d2998fa5
使用地面观测、遥感和机器学习
·发布于Towards Data Science ·阅读时长 11 分钟·2023 年 8 月 21 日
–
作者：Stephen Klosterman 和 Earthshot 科学团队。内容最初在2022 年 12 月美国地球物理联盟秋季会议上展示。
介绍
生态恢复项目通常需要投资以启动活动。为了为森林增长和保护项目创造碳融资机会，必须能够预测木质生物量中的碳积累，或在防止森林砍伐的情况下避免的排放。此外，还需要尝试理解广泛的其他生态系统属性的可能变化，例如植物和动物物种组成及水质。为了创建碳积累预测，常见的方法是对特定地点的项目给予个别关注和研究努力，这些项目可能分布在全球各地。因此，拥有一张局部准确且全球适用的生长率或其他感兴趣参数值的地图，将便于快速“勘探”生态系统恢复机会。在这里，我们描述了创建这样的地图的方法，该地图源于基于先前发布的文献综述数据训练的机器学习模型。随后，我们展示了如何在 Google Earth Engine 应用中实施该地图。
数据与方法
我们使用了一份最近发布的数据集，该数据集包含森林群体生物量测量值、年龄和地理位置（Cook-Patton et al. 2020），用于训练一个机器学习模型，以预测常用的查普曼-理查兹（CR）生长函数的一个参数。
在清理了异常值和不现实观测数据后，我们剩下了大约 2000 个观测值，如下图所示，符号大小与每个地点的观测数量成正比：
基于地点的全球数据分布；符号大小与每个地点的测量数量成正比。作者提供的图片。
观测数据分布在 390 个地点。大多数地点（64%）只有一个测量点，而有一个地点有 274 个测量点。
Cook-Patton 等（2020）将这些数据与来自美国和瑞典的其他库存数据结合，创建了全球碳积累速率地图。然而，该工作假设了一个线性碳积累模型，而有更多生物学上更为现实的替代模型。这里我们展示了如何为曲线拟合方程（CR 函数）创建全球参数层。我们的方法类似于 Chazdon 等（2016）的工作，创建了拉丁美洲热带地区 Michaelis-Menten 方程的参数值地图。然而，这里我们没有将曲线拟合参数模型限制为环境协变量的线性组合，而是将曲线拟合参数作为统计学习方法（XGBoost）的响应变量，以捕捉特征之间的任何非线性或交互行为。我们选择用 CR 方程表示生长，因为它是一个灵活的曲线，可以呈现 S 形或对数形状，是林业行业中一种简单而生物学上现实的树木生长模型的标准方法（Bukoski 等，2022）：
Chapman-Richards（CR）方程。作者提供的图片。
其中y是时间t的生物量，y_max是森林成熟时的生物量上限，b控制时间 0 时的生物量，m是形状参数，k是我们估计的生长参数。在这项工作中，我们假设b = 1，将生物量限制为零，m = ⅔，与类似工作中的做法一致（Bukoski 等，2022）。这就只剩下y、y_max、t和k作为未知量。以下是几个示例 CR 曲线，m = ⅔，b = 1，y_max = 1，和不同的k值：
作者提供的图片
参数k明显影响达到接近最大潜在生物量水平所需的时间：
作者提供的图片
我们从 Walker 等（2022）的潜在生物量图中提取了相关地点的y_max值。与这篇出版物一起提供的地图包括了当前和最大潜在生物量的预测，涵盖了各种假设和条件。特别是，我们使用了这里提供的 Base_Pot_AGB_MgCha_500m.tif 地图。通过为每个测量点估算最大生物量，我们将 Cook-Patton 等（2020）中的生物量和森林立地年龄的配对值分别代入y和t，并计算了该数据集中每个测量值的k值。以下是我们获得的k值分布，显示出一个广泛的范围和较长的右尾：
由作者提供的图片
在解读这些数据时，重要的是要注意，由于我们使用了树木生长数据源，我们的模型将适用于自然再生森林。最近已经完成了对单一物种种植园的文献回顾（Bukoski 等，2022），并且正在进行对农林系统（Cook-Patton等，正在准备中）和多样化树种种植项目（Werden等，正在准备中——请联系 Leland Werden（leland@crowtherlab.com）参与公民科学文献回顾，特别是如果你会讲除英语外的其他语言）。一个类似的相关工作涉及对自然再生和单一栽培在缓解气候变化方面相对有效性的经济计量比较（Busch等，正在准备中）。
我们为我们的模型使用了 61 个空间显式特征，包括生物群落（在进行一次性编码后 12 个特征）、来自 SoilGrids 项目的土壤属性（11 个特征）、Terraclimate提供的月度气候数据（1960–1990 年期间的 14 个特征）以匹配Bioclim数据（19 个特征），以及高程、坡度、方位和阴影的地形特征（4 个特征）。这些特征类似于其他研究中应用机器学习来绘制碳积累速率（Cook-Patton 等，2020）或其他相对时间不变的生态系统属性（如某地点的最大潜在生物量）（Walker 等，2022）所使用的特征。
我们使用了XGBoost来构建多个回归模型以预测k并探索响应变量与潜在特征列表之间的关系。为了帮助解释，我们的目标之一是选择一个特征尽可能少但性能良好的模型。我们使用SHAP（SHapley Additive exPlanation）值来确定特征的重要性，并发现模型性能和选定特征的显著差异，取决于用于训练与测试的数据。因此，我们对所有模型开发使用了 10 折交叉验证，而不是留出单独的测试集。我们使用GroupKFold以确保来自同一地点的测量不会在折叠之间被分割；换句话说，没有任何一个折叠在训练数据和测试数据中都有来自同一地点的数据，从而减少模型评估中的空间自相关效应。
-
修剪相关预测变量： 为了开始选择特征，我们最初根据平均绝对 SHAP 值对所有特征进行排序。为此，我们为每个 10 个训练折叠训练了一个模型，然后计算相应验证折叠的平均绝对 SHAP 值，以查看特征在训练样本之外进行预测时的重要性。我们对所有 10 个折叠中每个特征的平均绝对 SHAP 值进行了汇总，并按降序排序以进行“排名选择投票”特征选择程序。从列表顶部的特征开始，按顺序进行，我们丢弃了所有 Pearson’s r > 0.8 且排名较低的预测变量，以减少多重共线性，并开始修剪特征集的过程。此步骤对最终模型性能的影响微乎其微，并将特征集减少到 41 个特征。
-
每折训练单独模型并结合见解： 使用这个较小的特征集，我们对 10 个折叠中的每一个分别进行了反向选择：在每次迭代中，我们根据验证集上的平均绝对 SHAP 值对特征进行排序。根据这个排序，我们确定了特征最少的模型，其生物量估计的 RMSE（均方根误差）在 1 Mg AGB/ha（每公顷地上生物量）内，约为 1%的差异——换句话说，最简单的模型几乎“与最佳模型一样好”。由于这个步骤，不同折叠的最佳模型具有不同的特征和特征数量。为了结合各折叠的见解，我们执行了类似于前一步的排名选择投票，以确定哪些特征应该考虑用于最终模型及其排序（27 个特征）。
-
使用每次折叠的共同特征集进行最终特征选择： 再次使用向后选择程序，我们检查了 10 个折叠的验证集上的均值 RMSE 和 R²，但对于每个折叠使用相同的特征集。我们发现模型验证性能在折叠间通常很嘈杂，这突显了需要额外的数据收集以开发更稳健的模型。我们还发现，相对较少特征的模型在验证得分上几乎与更多特征的模型一样好。
作者提供的图像
我们选择了一个具有 5 个特征的模型（最高验证 R²，RMSE 略低于第二低）。最终的向后选择（步骤 3）的结果如上所示，而最终模型的 SHAP 值（在验证集上收集）如下所示。像这里展示的 SHAP 汇总图，也称为“蜜蜂群”图，有一行代表每个特征，在每行中有一个数据集（此处为验证数据）中的每个预测点。点的颜色表示特征值，垂直偏移表示该 SHAP 值的数据密度，x 轴坐标表示 SHAP 值，或该样本对模型预测的有符号影响（Lundberg 和 Lee，2017）。SHAP 汇总图显示了特征的高值或低值是否对模型预测产生了正面或负面的影响及其大小。
作者提供的图像
为了实现这一模型以创建预测的* k *值的全图，我们利用了我们的开源库earthshot_model_utils来导出空间平铺预测因子的 CSV 文件，执行模型推断，然后将结果导出并合并成一个单一的 geotiff（见下一节）。
讨论与结论
与相关研究工作的比较，我们的模型表现似乎相似。交叉验证 R²的平均值（±标准误差）为 0.485 ± 0.04，而使用这些及其他数据进行训练的线性碳积累率模型的 20%保留测试集的 R²为 0.445（Cook-Patton 等，2020）。
在模型的五个特征中，最重要的特征是等温性，它表示某个区域昼夜温差相对于夏冬温差的大小。其值范围从 0 到 100，其中昼夜温差表示为年温差的百分比。这种昼夜温差与年温差的关系被认为影响物种分布，并且对“热带、岛屿和海洋环境”有用（O’Donnell 和 Ignizio，2012）。它可能是环境支持树木生长的一个指标，较大的特征值（相对于年温差较大的昼夜温差）会导致较高的k预测（相对于最大值的较快增长，如上图的 SHAP 值图所示）。等温性在赤道附近相对较大，年温差最小，因此也一般代表较暖的温度。从剩余的特征来看，SHAP 值清楚地表明，树木在风速较小的地方以及土壤湿度较高的地方（无论是物理上还是生物上）预测生长速度较快，这些都是模型行为上合理的。与 Palmer 干旱严重性指数和土壤碳的关系则更加复杂和互动，我们在这里不再详细讨论。
使用 earthshot_model_utils，我们使用训练好的模型在 1 公里空间分辨率下对非洲进行了推断，生成了* k* 值的大陆地图。将这些值代入 CR 公式，并结合潜在的生物量y_max，我们可以生成森林再生开始后的任何时间点的 AGB 地图。我们的 Google Earth Engine 应用 显示了生长开始后 10 年、20 年和 30 年的生物量；30 年的生物量如下图所示。这个数据产品可以用于几乎即时的自然森林再生预测，具有高空间分辨率和大范围地理尺度，能够快速评估碳项目。
这是对自然再生开始 30 年后的生物量的估算。图片由作者提供。
虽然我们选择使用地面数据来创建这项工作的生物量积累模型，但我们团队也在投入研究其他可能的方法。这些方法包括更多依赖遥感观测而非地面观测的统计建模方法，以及基于过程的模型。遥感指标如冠层高度和生物量使用 LiDAR（激光雷达）或 SAR（合成孔径雷达）技术进行测量，这些技术通常没有多时相组件，因此无法跟踪生物量的增长轨迹。学术研究人员和公司正在利用多时相遥感数据（如 Landsat）作为特征构建机器学习模型，以有效地在空间（例如 Walker 等人 2020 年）和时间上（如CTREEs、Chloris Geospatial 和Sylvera的努力）外推 LiDAR 测量。这开辟了广阔的可能性，但需要注意从太空评估树木大小的局限性和必要的地面验证。基于过程的建模将提供另一种有价值的视角，特别是其在建模生物多样性和其他生态系统服务细节方面的能力（参见例如 Fisher 和 Koven, 2020），并努力与遥感数据结合以提供基准数据（Ma 等人 2022）。
为了继续这里描述的工作，我们计划创建模型预测的不确定性地图，使用多种基于自助抽样的模型，以及蒙特卡洛方法来考虑最大潜在生物量的不确定性。我们计划在未来的工作中探索和比较多种建模技术，利用多个基于地面数据集的方法。
参考文献
Bukoski, J.J., Cook-Patton, S.C., Melikov, C. 等人. 全球单一栽培森林中地上碳积累的速率及驱动因素. 《自然通讯》13, 4206 (2022). doi.org/10.1038/s41467-022-31380-7
Chazdon, R.L., Broadbent, E.N., Rozendaal, D.M.A. 等人. 拉丁美洲热带地区二次生长森林再生的碳封存潜力. 《科学进展》2, 5 (2016). doi.org/10.1126/sciadv.1501639
Cook-Patton, S.C., Leavitt, S.M., Gibbs, D. 等人. 全球自然森林再生的碳积累潜力地图. 《自然》585, 545–550 (2020). doi.org/10.1038/s41586-020-2686-x
Fisher, R.A., 和 Koven, C.D. 关于陆地表面模型未来发展及代表复杂陆地系统的挑战的观点. 《先进地球系统模型杂志》12, 4 (2020). doi.org/10.1029/2018MS001453
Lundberg, S.M. 和 Lee, Su-In。《统一模型预测解释方法》。Advances in Neural Information Processing Systems 30 (2017)。 proceedings.neurips.cc/paper_files/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf
Ma, L., Hurtt, G., Ott, L. 等。《生态系统人口模型 (ED v3.0) 的全球评估》。Geosci Model Dev 15, 1971–1994 (2022)。 doi.org/10.5194/gmd-15-1971-2022
O’Donnell, M.S. 和 Ignizio, D.A. 《支持美国大陆生态应用的生物气候预测因子：美国地质调查局数据系列 691(2012)》。 pubs.usgs.gov/ds/691/ds691.pdf
Walker, W.S., Gorelik, S.R., Cook-Patton, S.C. 等。《陆地上碳存储增加的全球潜力》。Proc Natl Acad Sci USA 119, 23 (2022)。 doi.org/10.1073/pnas.2111312119
最初发布于 https://www.earthshot.eco。
《交通拥堵分析：使用图论》
学习如何使用图论找到城市基础设施中的潜在关键点
·
关注 发布于 Towards Data Science ·10 分钟阅读·2023 年 8 月 19 日
–
图论在现实问题中有很多应用，例如社交网络、分子生物学或地理空间数据。今天，我将展示最后一种应用，即分析城市的道路布局，以预测关键街道、交叉口，以及基础设施的变化如何影响这些因素。但首先，让我们从基础知识开始。
图及其中心性度量
图是顶点及其边的集合：
集合 E 是无序元组 (x, y) 的子集，其中 x 和 y 是图的顶点，且 x 不等于 y。[图像由作者提供]
边表示节点之间的连接。如果边没有方向，我们称之为无向图。无向图的一个现实例子可以是化学分子，其中顶点是原子，化学键表示为边。
血清素分子是一个简单无向图的例子。 [来源]
然而，有时我们需要了解边是从 u 到 v，从 v 到 u，还是双向的。例如，如果马克喜欢爱丽丝，并不一定意味着这是双向的（ ☹ ）。在这些情况下，我们可以将边定义为有序元组而不是无序元组。
方括号表示公式中的无序元组，而圆括号表示有序元组。[图片由作者提供]
人际互动可以使用有向图来描述。[图片由作者提供]
使用图结构，我们可以定义一个中心性度量。这是一个用于回答以下问题的指标：
这个顶点/边在图中的重要性有多大？
而且有许多方法可以回答这个问题。
评估图组件重要性的不同方法
根据任务的不同，我们可以从不同的角度来评估中心性。最常见的度量指标有：度、紧密性和中介性。我们将使用扎卡里·卡拉泰俱乐部图来讨论这些指标 [更多信息]。它展示了不同空手道俱乐部成员之间的联系。你可以在这里找到用于生成图片的代码。
度中心性
最基本的中心性。它仅定义在顶点上，等于顶点的度数（即邻接顶点的数量）。例如，我们可以回想一下人际关系图，在人际友谊的情况下，这个指标将回答以下问题：
“这个人有多受欢迎？”
空手道俱乐部图的节点度中心性。中心性度量标准按图的最大度数（即节点数减一）进行了归一化。[图片由作者提供]
图中的路径
对于接下来的两个中心性，我们需要将一些概念引入到图论知识中。所有这些都是非常直观的，从边的权重开始。我们可以给边添加权重，以标记它们之间的差异。例如，在交通图的情况下，这可以是道路长度。
在图中，我们可以定义路径，即从 A 到 B 需要遍历的顶点列表。路径中的连续顶点是邻居，第一个顶点是 A，最后一个是 B。路径距离是沿途边的权重之和。A 和 B 之间的最短路径是距离最小的路径。
A 和 F 之间的最短路径是 [A, C, E, D, F]，距离为 20。 [source]
近似中心性
拥有所有这些新知识后，我们可以回到我们的指标。下一个是近似中心性，它告诉我们一个节点离图中其余部分的距离。它对特定顶点的定义是图中所有其他顶点的最短路径的平均值的倒数。这样，较短的平均路径转化为更高的近似中心性。
Karate Club 图中的节点近似中心性。 [图片由作者提供]
介数中心性
介数中心性提供了图中哪些节点对交通流量至关重要的信息。想象一个拥有广泛道路网络的城市，其中每个交汇点都是一个节点。这些节点中的一些在日常通勤中作为关键连接点，而其他的可能是交通流量几乎没有影响的死胡同。前者的介数中心性得分较高，计算方法是通过交汇点的最短路径的比例。
Karate Club 图中的节点介数中心性。 [图片由作者提供]
城市规划作为图
现在，我们拥有了描述和分析图的工具，可以开始将城市规划提取为图形形式。为此，我们可以使用 Open Street Maps (OSM)，通过 osmnx 库将其导入 Python 作为 NX 图。我们将从一个较小的示例开始，讨论为了提高工作效率和时间的处理过程。
Grzegórzki 是克拉科夫市的十八个区之一，有两个复杂的环形交叉口——Mogilskie 和 Grzegórzeckie，以及许多交汇点。因此，我们将能够看到数据工程中的大多数潜在陷阱。
Grzegórzki 的行政边界。 [©Google]
让我们从将 OSM 仓库中的数据导入到 Python 图中，并绘制结果开始：
原始 OSM 数据导入。白点是节点，代表道路交汇点。 [图片由作者提供]
这个图表有些问题——你能发现是什么吗？
我们为单一道路部分获取了多个边缘，结果图形中有近 3000 个“交汇点”。这并没有提供正确的表示（我们不能在路中间掉头，并且每个节点都会使计算变得更慢）。为了解决这种情况，我们将通过删除两个交汇点之间道路上的所有节点来进行图拓扑简化。在 OSMnx 中，我们有一个名为*ox.simplify_graph()*的函数来实现这一点。
拓扑简化后的道路布局。现在每个节点表示道路交叉点。[作者提供的图片]
还有一个问题——如你所见，我们大多数道路都有两个边缘，每个方向一个。因此，每个交汇点都有多个节点，这是一种不希望出现的行为。想象一下我们在一个交汇点上，我们要左转，而没有专用的左转车道（或者已经满了）。只要我们不能完成转弯，其他车辆就会被阻挡。在我们当前的图形中，这不是真实的。左转由两个独立的节点组成，一个是左转节点，另一个是跨越对面车道的节点。这会表示这些是两个独立的操作，但实际上并非如此。
这就是为什么我们要合并交汇点，这意味着我们将把相互靠近的多个节点合并为一个。我们会选择一个足够大的合并半径，以将交汇点的多个部分合并为一个，但另一方面保持环形交汇点为多个节点结构，因为它们只能部分被阻塞。为此，我们将使用 osmnx 函数ox.consolidate_intersections()。
交汇点合并后的道路布局。[作者提供的图片]
交汇点的比较。之前和之后。[作者提供的图片]
在这些操作之后，我们几乎准备好进行分析了。最后一个问题是克拉科夫的市政边界——由于许多人从邻近城镇旅行，而图形分析仅包括图形中的数据，我们需要包含这些区域。我将在下一章中介绍不这样做的影响。这里是我们的图形：
颜色表示最大速度。颜色越亮，值越高。我们可以看到 A4 高速公路用黄色标记。大多数道路，用蓝色标记，为 50 km/h。[作者提供的图片]
你可以在这个jupyter notebook中找到用于生成此地图的源代码以及下一章中使用的所有图形。
道路布局的介数中心性
在这个案例研究中，我们将专注于使用介数中心性测量来估计道路交通。未来，这可能会扩展到图论中的其他技术，包括GNN（图神经网络）使用。
我们将从计算道路布局表示中所有节点和边的介数中心性开始。为此，我们将使用NetworkX库。
克拉科夫每个道路段的介数中心性。 [作者提供的图片]
由于图中道路数量众多，很难看出哪些组件在交通中最关键。让我们查看图的中心性度量分布。
克拉科夫道路布局中街道和交叉口的中心性度量分布。 [作者提供的图片]
我们可以使用这些分布来过滤掉不太重要的交叉口和街道。我们将选择每个前 2% 的部分，其阈值如下：
-
节点的中心性为 0.047，
-
边的中心性为 0.021。
中心性度量在阈值处理后的图。 [作者提供的图片]
我们可以看到，介数中心性最高的道路段是：
-
A4 高速公路和 S7 作为克拉科夫的环城高速（注意克拉科夫没有北部环城道路），
-
第二环路的西部以及其与 A4 的连接，
-
第三环路的北部（替代缺失的北部环城道路），
-
Nowohucka 街道连接第二环路和城市的东北部，
-
Wielicka 道路从市中心通向东南部高速公路部分。
让我们将这些信息与 Google Maps 上克拉科夫的实际交通地图进行比较：
克拉科夫周一通勤的典型交通 [©2023 Google, source]
我们可以看到我们的见解与交通雷达的结果相符。其背后的机制很简单——那些在图中具有高介数中心性的组件是最常用来通行最短路径的。如果驾驶员选择最佳路线，那么交通量最高的街道和交叉口将具有最高的介数中心性。
让我们回到图形工程的最后部分——扩展图形边界。我们可以检查如果仅将城市边界纳入分析会发生什么：
克拉科夫道路的介数中心性，未考虑邻近城镇。 [作者提供的图片]
A4 高速公路，由于其环城性质，是最重要的组件之一，但在整个图中的中心性度量却是最低的！这是因为 A4 位于城市的边缘，大部分交通来自外部，我们不能将这一因素纳入介数中心性。
如何使用介数中心性分析布局变化对交通的影响，
让我们来看看图分析的不同场景。假设我们想预测道路封闭（例如由于事故）如何影响交通。我们可以使用中心性测量来比较两个图之间的差异，从而检查中心性的变化。
在本研究中，我们将模拟 A4–7 高速公路段上的汽车事故，这是一个常见的情况。事故将导致该段完全封闭。
我们将通过从图中去除 A4–7 段并重新计算中心性测量来创建一个新的道路网络。
新布局的中心性测量。红色的 A4 部分代表缺失部分。[作者提供的图片]
让我们来看看中心性分布：
去除 A4–7 高速公路段后的克拉科夫道路布局中街道和交叉口的中心性测量分布。[作者提供的图片]
我们可以看到它仍然与原始情况非常相似。为了检查中心性测量的变化，我们将计算残差图，其中中心性测量是原始道路布局与事故后之间的差异。正值将表示事故后的更高中心性。在其中一个图中缺失的节点和交叉口（如 A4–7）将不会被包含在残差图中。以下是残差的测量分布：
去除 A4–7 高速公路段后的中心性变化分布。[作者提供的图片]
再次，我们将筛选出受影响的前 2%的街道和节点。这次的阈值是：
-
节点的测量值为 0.018，
-
边的测量值为 0.017。
去除 A4–7 高速公路段后，介于中介中心性增加最大的街道和交叉口。[作者提供的图片]
我们可以看到，连接环路分段到市中心的道路有所增加，其中第二环路位于城市的西侧，包含两座横跨维斯瓦河的桥梁之一的变化最大。
图的中心性分析无法在道路网络上实现的内容
在图分析中，有一些我们无法考虑的因素。在本分析中我们可以看到的两个最重要的因素是：
- 图的中心性分析假设节点之间的流量分布是均匀的。
这在大多数情况下是错误的，因为乡村和城市的人口密度不同。然而，还有其他因素可以减少这种影响，例如，相比于生活在城市中心的人，居住在邻近乡村的人更倾向于选择汽车作为通勤方式。
- 图分析仅考虑图中存在的事物。
在提供的示例中，这一点不易察觉，尤其是对克拉科夫以外的人来说。我们来看一下Zakopianka。它是连接市中心和克拉科夫南部大多数市镇的主要交通干道，也是贯穿全国的 DK7（国家公路 7 号）的组成部分。
DK7 公路 — 绿色部分表示快速路。 [source]
如果我们将克拉科夫的 DK7 典型交通情况与我们的中心性测量进行比较，它们完全不同。平均介数中心性约为 0.01，这比前 2%的阈值小两倍。然而在现实中，它是最拥堵的路段之一。
Zakopianka 的平均拥堵情况与介数中心性进行比较。 [©2023 Google, source]
总结
图论及其分析在多个场景中都有应用，例如本研究中展示的交通分析。通过对图进行基本操作和度量，我们可以在比建立整个模拟模型更短的时间内获得有价值的见解。
这个完整的分析可以通过几十行 Python 代码来完成，并且不限于一种道路布局。我们也可以很容易地过渡到图论的其他分析工具。
像所有事物一样，这种方法也有其缺点。主要缺点是对均匀交通分布的假设以及范围仅限于图结构。
包含本研究中使用的代码的 Github 仓库可以在这里找到。
使用 MapReduce 进行大规模数据处理
深入探讨 MapReduce 和并行化
·发布于 Towards Data Science ·4 分钟阅读·2023 年 7 月 19 日
–
图片由 Luca Nicoletti 提供，发布在 Unsplash
在当前的市场环境中，组织必须进行数据驱动的决策，以保持竞争力并促进创新。因此，每天都会收集大量数据。
尽管数据持久性问题在很大程度上已得到解决，这要归功于云存储的广泛可用性和价格实惠，现代组织仍然面临着高效有效处理大量数据的挑战。
在过去几十年中，出现了许多编程模型来解决大规模处理大数据的挑战。毫无疑问，MapReduce 是最受欢迎和有效的方法之一。
什么是 MapReduce
MapReduce 是一个分布式编程框架，最初由 Jeffrey Dean 和 Sanjay Ghemawat 于 2004 年在 Google 开发，并受到函数式编程基本概念的启发。他们的提案涉及一个包含两个步骤的并行数据处理模型；map 和 reduce。
简单来说，map 步骤涉及将原始数据分割成小块，以便对每个数据块应用转换逻辑。因此，可以在创建的块上并行处理数据，最后，reduce 步骤将汇总/整合处理过的块，并将最终结果返回给调用者。
MapReduce 算法如何工作
尽管 MapReduce 算法通常被认为是一个两步过程，但它实际上包含三个不同的阶段。
1. Map： 在这个第一步骤中，数据被拆分成更小的块，并分布到通常属于处理单元集群的多个节点上。每个创建的块被分配给一个 mapper。Mapper 的输入是一组 <key, value> 对。数据处理执行后（依然是 <key, value> 形式），Mapper 会将生成的输出写入临时存储。
例如，我们可以考虑以下例子，其中输入文本首先被拆分到三个 Mapper 上，输入以键值对的形式提供。
MapReduce 算法的映射步骤 — 来源：作者
2. Shuffling： 在这个步骤中，算法将数据洗牌，以便具有相同键的记录被分配到相同的工作节点。这通常是整个 MapReduce 过程生命周期中最昂贵的操作。
MapReduce 中的洗牌步骤 — 来源：作者
3. Reduce： 在这最后一步中，每个 Reducer 将接受对应 Mapper 输出的 <key, value> 对作为输入。所有具有相同键的 Mapper 输出将分配给相同的 Reducer，Reducer 将汇总这些值，并将汇总结果以 <key, value> 对的形式返回。
Reduce 步骤中的 MapReduce — 来源：作者
MapReduce 和 Hadoop
MapReduce 是 Apache Hadoop 框架的一部分，用于访问存储在 Hadoop 分布式文件系统（HDFS）中的数据。Hadoop 由四个基本模块组成：
-
Hadoop Distributed File System (HDFS)：这是一个分布式文件系统，可以以容错的方式存储大型数据集
-
Yet Another Resource Negotiation (YARN)：这是一个节点管理器，监控集群和资源，同时也作为作业调度器。
-
MapReduce
-
Hadoop Common：这是一个提供常用 Java 库的模块
之前我们提到过，Mapper 和 Reducer 在计算机集群的独立节点上运行。实际上，这些工作节点是 Hadoop 框架的一部分，该框架决定了每种情况下所需的 Mapper 数量，这取决于输入大小的体积。
Hadoop 设计时提供了容错功能。在节点发生故障时，Hadoop 会在另一个映射节点上重新运行任务并生成所需的输出。
最后的思考
MapReduce 在分布式计算中是一个突破性的概念，使许多组织能够处理大量数据并提取有价值的洞察。
熟悉这个概念至关重要，特别是在利用如 Spark 等依赖于 MapReduce 框架的技术时。
👉 成为会员 ，在 Medium 上阅读所有故事。您的会员费直接支持我和您阅读的其他作者。您还将获得对 Medium 上每个故事的全面访问权限。
[## 通过我的推荐链接加入 Medium — Giorgos Myrianthous
成为 Medium 会员后，您的会员费的一部分将分配给您阅读的作者，并且您可以全面访问每个故事…
gmyrianthous.medium.com](https://gmyrianthous.medium.com/membership?source=post_page-----f0d8776d0fcf--------------------------------)
3 月版：数据与因果关系
每月版
数据科学家如何处理因果推断
·
关注 发表在 Towards Data Science ·4 分钟阅读·2023 年 3 月 2 日
–
照片由 Joey Genovese 提供，刊登在 Unsplash
在最近的作者焦点问答中，Matteo Courthoud 反映了在工业界或学术界工作时，做出可靠预测的重要性日益增加：
我认为在未来，因果推断将变得越来越重要，我们将看到社会科学的理论方法与计算机科学的数据驱动方法之间的趋同。
我们希望你能阅读我们生动对话的其余部分；与此同时，Matteo 的观察激励我们深入档案，寻找关于因果推断及更广泛的因果主题的其他有见地的文章。我们在本月版中分享的选集内容从入门到更高级别，展示了数据科学和机器学习从业者在工作中每天使用的一些不同方法。
我们希望你喜欢探索这些推荐阅读！一如既往，我们很感激你将 TDS 作为你学习旅程的一部分；如果你希望以其他方式支持我们的工作（并在此过程中获得我们整个档案的访问权限），请考虑成为 Medium 会员。
[TDS 编辑]（https://medium.com/u/7e12c71dfa81?source=post_page-----396b2881aea9--------------------------------）
TDS 编辑亮点
-
因果关系的科学与艺术（第一部分）（2023 年 1 月，11 分钟）
Quentin Gallea, PhD对因果推断基础知识的易懂介绍旨在回答两个关键问题：理解因果关系为何如此重要，因果关系为何如此难以评估？
-
使用合成控制的因果推断（2022 年 11 月，9 分钟）
在 A/B 测试不是建立因果关系的好选择的情况下，准实验技术可能是解决方案。diksha tiwari提出合成控制作为一种替代方案。
-
通过回归的因果效应（2023 年 1 月，8 分钟）
Shawhin Talebi已覆盖因果推断的理论与实践超过一年。你可以回到这个系列的起点，或直接跳转到这篇关于回归技术及如何利用它们建立变量之间关系的最新文章。
-
因果推断的事件研究：应该做与不应该做的（2022 年 12 月，17 分钟）
事件研究在准实验背景中是另一种有帮助的方法；正如Nazlı Alagöz在这篇深入说明的文章中指出的，有许多陷阱需要避免，以便我们不从数据中得出错误的见解。
-
使用 SciPy 进行两组独立样本均值的统计显著性检验（2022 年 11 月，8 分钟）
对于那些热衷于动手实验数据的读者，Zolzaya Luvsandorj的教程是理想的起点：它提供了进入假设检验的实用入口，并包括所有所需的 Python 代码。
-
识别：可信因果推断的关键（2023 年 2 月，8 分钟）
正如Murat Unal在一篇深刻的新文章中解释的那样，“没有明确的识别，没有任何复杂的建模或估计能够帮助我们从数据中确定因果关系。” 阅读 Murat 的概述，以更好地理解识别及其重要性。
原创特点
探索我们最新的问答和阅读推荐。
-
“我写作的主要驱动力一直是学习**”** 我们与Matteo Courthoud的问答，他反思了离开学术界、对因果推断的兴趣以及公开写作的价值。
-
面对数据偏差依然困难且必要 我们分享了在机器学习和人工智能领域中持续主导对话的话题的核心阅读资料。
-
如何作为数据科学家培养良好习惯 在最近的一次总结中，我们汇总了经验丰富的数据从业者提供的技巧和策略，以便更流畅和可靠的工作流程。
热门文章
如果你错过了，这里是上个月 TDS 中最受欢迎的一些文章。
-
ChatGPT 的工作原理：模型背后的技术 由Molly Ruby撰写
-
你可能在不知不觉中成为高级 Pythonista 的 5 个迹象 由Bex T.撰写
-
使用 OpenAI 和 Python 提升你的简历：一步一步的指南 由Piero Paialunga撰写
-
如何用 Python 构建 ELT 作者：玛丽·陈
-
如何创建有效的自学计划，以成功自学数据科学 作者：马迪逊·亨特
-
你还在使用肘部法则吗？ 作者：萨穆埃尔·马赞提
-
如何为你的数据找到最佳理论分布 作者：厄尔多安·塔斯克森
我们很高兴在二月迎来了全新的 TDS 作者团队——他们包括 萨曼莎·霍德、阿尔瓦罗·佩尼亚、泰米托普·索博杜、弗雷德里克·霍特尔、吉尔·肖姆龙、拉斐尔·比绍夫、肖恩·史密斯、布鲁诺·阿尔维西奥、乔里斯·盖林、德米特里·埃柳塞耶夫、科里·贝克、波尔·马林、皮奥特·拉赫特、布鲁诺·波尼 和 诺布尔·阿克森 等。如果你有有趣的项目或想法要与我们分享，我们非常乐意听取你的意见!
下个月见。
以规模化方式掌握语义搜索：使用 FAISS 和 Sentence Transformers 在闪电般的推理时间内索引数百万份文档
深入了解一个高性能的语义搜索引擎的端到端演示，利用 GPU 加速、高效的索引技术和强大的句子编码器处理多达 100 万份文档的数据集，实现 50 毫秒的推理时间
·发表于 Towards Data Science ·阅读时间 15 分钟·2023 年 3 月 31 日
–
介绍
在搜索和信息检索领域，语义搜索已经成为一场变革。它使我们能够根据文档的意义或概念而不仅仅是关键字匹配来进行搜索和检索。与传统的基于关键字的搜索方法相比，语义搜索能够提供更复杂、更相关的结果。然而，挑战在于将语义搜索扩展到处理大量文档的语料库，而不会因分析每个文档的语义内容的计算复杂性而不堪重负。
在这篇文章中，我们迎接了通过利用两种前沿技术来实现可扩展语义搜索的挑战：FAISS 用于高效的语义向量索引，Sentence Transformers 用于将句子编码为这些向量。FAISS 是一个出色的库，旨在快速检索高维空间中的最近邻，使得即使在大规模下也能迅速进行语义最近邻搜索。Sentence Transformers，一个深度学习模型，生成句子的密集向量表示，有效捕捉其语义含义。
本文展示了如何利用 FAISS 和句子变换器的协同作用，构建一个具有卓越性能的可扩展语义搜索引擎。通过将 FAISS 和句子变换器集成，我们可以对来自大量文档的语义向量进行索引，从而实现快速准确的语义搜索体验。我们的方法可以实现新的应用，如上下文化问答和先进的推荐系统，当搜索 1M 文档的语料库时推理时间低至 50 毫秒。我们将指导你实现这一最先进的端到端解决方案，并在基准数据集上展示其性能。
图 1：搜索是你在这个世界上导航所需的一切（来源）
本文属于“大语言模型纪实：导航 NLP 前沿”系列，这是一系列每周更新的文章，将探讨如何利用大模型的力量来完成各种 NLP 任务。通过深入这些前沿技术，我们旨在赋能开发者、研究人员和爱好者，充分利用 NLP 的潜力，开启新的可能性。
迄今已发布的文章：
如往常一样，代码可在我的Github上找到。
句子变换器用于语义编码
深度学习带来了句子变换器的力量，它们制作密集的向量表示，捕捉句子意义的本质。这些模型在大量数据上进行训练，生成上下文化的词嵌入，旨在准确重建输入句子，并将语义相似的句子对拉近。
为了利用句子变换器在语义编码中的潜力，你需要首先选择一个合适的模型架构，如 BERT、RoBERTa 或 XLNet。确定模型后，我们将把文档语料库输入其中，为每个句子生成固定长度的语义向量。这些向量是句子核心主题和话题的紧凑数值表示。
以两个句子为例：“狗追逐猫”和“猫追逐狗”。通过句子变换器处理后，它们的语义向量将紧密相关，即使词序不同，因为基本意义相似。另一方面，像“天空是蓝色的”这样的句子会产生更远的向量，因为其含义不同。
使用句子变换器对整个语料库进行编码，我们得到一组语义向量，这些向量概括了文档的总体含义。为了使这种变换后的表示准备好进行高效检索，我们使用 FAISS 对其进行索引。敬请关注，我们将在下一节中深入探讨这个话题。
FAISS 用于高效索引
FAISS 支持多种索引结构，以优化不同的使用场景。它是一个设计用于在大量向量中快速找到与给定查询向量最接近的匹配项的库。
-
倒排文件（IVF）：对相似向量的簇进行索引。适用于中维向量。
-
产品量化（PQ）：将向量编码到量化子空间中。适用于高维向量。
-
基于簇的策略：将向量组织成分层的簇集以进行多级搜索。适用于非常大的数据集。
要使用 FAISS 进行语义搜索，我们首先加载我们的向量数据集（来自句子变换器编码的语义向量）并构建 FAISS 索引。我们选择的具体索引结构取决于语义向量的维度和期望的效率等因素。然后，我们通过将语义向量传递到 FAISS 索引中对其进行索引，FAISS 将高效地组织它们以实现快速检索。
对于搜索，我们将新的句子编码成一个语义向量查询并将其传递给 FAISS 索引。FAISS 将检索最接近的语义向量并返回最相似的句子。与线性搜索相比，后者会将查询向量与每个索引向量进行比对，FAISS 能够提供更快的检索时间，通常与索引向量的数量呈对数级缩放。此外，这些索引具有高度的内存效率，因为它们压缩了原始密集向量。
倒排文件索引
FAISS 中的倒排文件（IVF）索引将相似的向量聚集到“倒排文件”中，适用于中维向量（例如，100–1000 维）。每个倒排文件包含相互接近的向量子集。在搜索时，FAISS 仅搜索与查询向量最接近的倒排文件，而不是搜索所有向量，即使有许多向量，也能实现高效搜索。
要构建 IVF 索引，我们需要指定倒排文件（簇）的数量以及每个倒排文件的最大向量数量。然后，FAISS 将每个向量分配给最近的倒排文件，直到没有倒排文件超过最大数量。倒排文件包含代表性点，这些点总结了它们内部的向量。在查询时，FAISS 计算查询向量与每个倒排文件代表性点之间的距离，并仅搜索与查询向量最接近的倒排文件以找到最匹配的向量。
例如，如果我们有 1024 维的图像特征向量并且想要在 100 万个向量中进行快速搜索，我们可以创建一个包含 1024 个倒排文件（簇）的 IVF 索引，每个倒排文件最多包含 1000 个向量。在这种方法中，FAISS 只会搜索与查询最接近的倒排文件，从而比线性搜索更快。
汇总所有内容
在本节中，我们将使用 FAISS 和 Sentence Transformers 构建一个可扩展的语义搜索引擎。我们将展示如何评估这种方法的性能基准，并讨论进一步的改进和应用。
可扩展的语义搜索引擎
为了构建一个可扩展的语义搜索引擎，我们首先初始化ScalableSemanticSearch类。该类负责使用 Sentence Transformers 对句子进行编码，并使用 FAISS 对其进行索引，以实现高效搜索。它还提供了保存和加载索引、测量时间和内存使用的实用方法。
semantic_search = ScalableSemanticSearch(device="cuda")
接下来，我们使用编码方法对大量文档进行编码，该方法返回一个语义向量的 numpy 数组。该方法还创建了一个索引和句子之间的映射，这在检索前几个结果时会很有用。
embeddings = semantic_search.encode(corpus)
现在，我们使用 build_index 方法构建 FAISS 索引，该方法以嵌入向量作为输入。该方法根据嵌入中的数据点数量创建 IndexIVFPQ 或 IndexFlatL2 索引。
semantic_search.build_index(embeddings)
基于数据集大小选择索引方法
我们定义了两种索引方法：L2 距离的精确搜索和产品量化与 L2 距离的近似搜索。我们还将讨论选择第一种方法用于较小数据集（少于 1500 个文档）和第二种方法用于较大数据集的原因。
1. L2 距离的精确搜索
L2 距离的精确搜索是一种精确搜索方法，它计算查询向量与数据集中每个向量之间的 L2（欧几里得）距离。该方法可以保证找到精确的最近邻，但对于大型数据集可能比较慢，因为它执行线性扫描。
使用案例： 这种方法适用于需要精确最近邻的小型数据集，并且计算成本不是问题。
2. 产品量化与 L2 距离的近似搜索
近似搜索与产品量化和 L2 距离是一种近似最近邻搜索方法，它结合了倒排文件结构、产品量化和 L2 距离，以高效地在大数据集中搜索相似向量。该方法首先使用 k-means（faiss.IndexFlatL2 作为量化器）对数据集进行聚类，然后应用产品量化来压缩残差向量。这种方法比起蛮力方法使用更少的内存，从而实现更快的搜索速度。
使用案例： 这种方法适用于不严格要求精确最近邻的大型数据集，主要关注搜索速度和内存效率。
选择不同方法的原因依据数据集大小
对于包含少于 1500 个文档的数据集，我们设置了 L2 距离的精确搜索方法，因为在这种情况下计算成本不是一个重大问题。此外，这种方法可以保证找到最近邻，这对于较小的数据集是可取的。
对于较大的数据集，我们倾向于使用近似搜索与产品量化和 L2 距离方法，因为它比精确搜索方法更高效，消耗的内存也更少。当优先考虑搜索速度和内存效率而非找到精确的最近邻时，近似搜索方法更适合大数据集。
搜索过程
在构建索引后，我们可以通过提供输入查询和返回的结果数量来执行语义搜索。search 方法计算输入句子与已索引嵌入之间的余弦相似度，并返回最匹配句子的索引和分数。
query = "What is the meaning of life?"
top = 5
top_indices, top_scores = semantic_search.search(query, top)
最后，我们可以使用 get_top_sentences 方法检索最顶级的句子，该方法接受索引到句子的映射和顶级索引作为输入，并返回顶级句子的列表。
top_sentences = ScalableSemanticSearch.get_top_sentences(semantic_search.hashmap_index_sentence, top_indices)
完整模型
我们模型的完整类如下：
class ScalableSemanticSearch:
"""Vector similarity using product quantization with sentence transformers embeddings and cosine similarity."""
def __init__(self, device="cpu"):
self.device = device
self.model = SentenceTransformer(
"sentence-transformers/all-mpnet-base-v2", device=self.device
)
self.dimension = self.model.get_sentence_embedding_dimension()
self.quantizer = None
self.index = None
self.hashmap_index_sentence = None
log_directory = "log"
if not os.path.exists(log_directory):
os.makedirs(log_directory)
log_file_path = os.path.join(log_directory, "scalable_semantic_search.log")
logging.basicConfig(
filename=log_file_path,
level=logging.INFO,
format="%(asctime)s %(levelname)s: %(message)s",
)
logging.info("ScalableSemanticSearch initialized with device: %s", self.device)
@staticmethod
def calculate_clusters(n_data_points: int) -> int:
return max(2, min(n_data_points, int(np.sqrt(n_data_points))))
def encode(self, data: List[str]) -> np.ndarray:
"""Encode input data using sentence transformer model.
Args:
data: List of input sentences.
Returns:
Numpy array of encoded sentences.
"""
embeddings = self.model.encode(data)
self.hashmap_index_sentence = self.index_to_sentence_map(data)
return embeddings.astype("float32")
def build_index(self, embeddings: np.ndarray) -> None:
"""Build the index for FAISS search.
Args:
embeddings: Numpy array of encoded sentences.
"""
n_data_points = len(embeddings)
if (
n_data_points >= 1500
): # Adjust this value based on the minimum number of data points required for IndexIVFPQ
self.quantizer = faiss.IndexFlatL2(self.dimension)
n_clusters = self.calculate_clusters(n_data_points)
self.index = faiss.IndexIVFPQ(
self.quantizer, self.dimension, n_clusters, 8, 4
)
logging.info("IndexIVFPQ created with %d clusters", n_clusters)
else:
self.index = faiss.IndexFlatL2(self.dimension)
logging.info("IndexFlatL2 created")
if isinstance(self.index, faiss.IndexIVFPQ):
self.index.train(embeddings)
self.index.add(embeddings)
logging.info("Index built on device: %s", self.device)
@staticmethod
def index_to_sentence_map(data: List[str]) -> Dict[int, str]:
"""Create a mapping between index and sentence.
Args:
data: List of sentences.
Returns:
Dictionary mapping index to the corresponding sentence.
"""
return {index: sentence for index, sentence in enumerate(data)}
@staticmethod
def get_top_sentences(
index_map: Dict[int, str], top_indices: np.ndarray
) -> List[str]:
"""Get the top sentences based on the indices.
Args:
index_map: Dictionary mapping index to the corresponding sentence.
top_indices: Numpy array of top indices.
Returns:
List of top sentences.
"""
return [index_map[i] for i in top_indices]
def search(self, input_sentence: str, top: int) -> Tuple[np.ndarray, np.ndarray]:
"""Compute cosine similarity between an input sentence and a collection of sentence embeddings.
Args:
input_sentence: The input sentence to compute similarity against.
top: The number of results to return.
Returns:
A tuple containing two numpy arrays. The first array contains the cosine similarities between the input
sentence and the embeddings, ordered in descending order. The second array contains the indices of the
corresponding embeddings in the original array, also ordered by descending similarity.
"""
vectorized_input = self.model.encode(
[input_sentence], device=self.device
).astype("float32")
D, I = self.index.search(vectorized_input, top)
return I[0], 1 - D[0]
def save_index(self, file_path: str) -> None:
"""Save the FAISS index to disk.
Args:
file_path: The path where the index will be saved.
"""
if hasattr(self, "index"):
faiss.write_index(self.index, file_path)
else:
raise AttributeError(
"The index has not been built yet. Build the index using `build_index` method first."
)
def load_index(self, file_path: str) -> None:
"""Load a previously saved FAISS index from disk.
Args:
file_path: The path where the index is stored.
"""
if os.path.exists(file_path):
self.index = faiss.read_index(file_path)
else:
raise FileNotFoundError(f"The specified file '{file_path}' does not exist.")
@staticmethod
def measure_time(func: Callable, *args, **kwargs) -> Tuple[float, Any]:
start_time = time.time()
result = func(*args, **kwargs)
end_time = time.time()
elapsed_time = end_time - start_time
return elapsed_time, result
@staticmethod
def measure_memory_usage() -> float:
process = psutil.Process(os.getpid())
ram = process.memory_info().rss
return ram / (1024**2)
def timed_train(self, data: List[str]) -> Tuple[float, float]:
start_time = time.time()
embeddings = self.encode(data)
self.build_index(embeddings)
end_time = time.time()
elapsed_time = end_time - start_time
memory_usage = self.measure_memory_usage()
logging.info(
"Training time: %.2f seconds on device: %s", elapsed_time, self.device
)
logging.info("Training memory usage: %.2f MB", memory_usage)
return elapsed_time, memory_usage
def timed_infer(self, query: str, top: int) -> Tuple[float, float]:
start_time = time.time()
_, _ = self.search(query, top)
end_time = time.time()
elapsed_time = end_time - start_time
memory_usage = self.measure_memory_usage()
logging.info(
"Inference time: %.2f seconds on device: %s", elapsed_time, self.device
)
logging.info("Inference memory usage: %.2f MB", memory_usage)
return elapsed_time, memory_usage
def timed_load_index(self, file_path: str) -> float:
start_time = time.time()
self.load_index(file_path)
end_time = time.time()
elapsed_time = end_time - start_time
logging.info(
"Index loading time: %.2f seconds on device: %s", elapsed_time, self.device
)
return elapsed_time
端到端演示
本节将提供使用 SemanticSearchDemo 类的可扩展语义搜索引擎的端到端演示。
上述代码中的主函数。目标是理解不同概念和组件如何结合创建一个实用的应用程序。
初始化 SemanticSearchDemo 类：要初始化 SemanticSearchDemo 类，需要提供数据集路径、ScalableSemanticSearch 模型、可选的索引路径和可选的子集大小。这种灵活性使得可以使用不同的数据集、模型和子集大小。
demo = SemanticSearchDemo(
dataset_path, model, index_path=index_path, subset_size=subset_size
)
加载数据：load_data 函数主动读取和处理文件中的数据，然后返回一个句子列表。系统使用这些数据来训练语义搜索模型。
sentences = demo.load_data(file_name)
subset_sentences = sentences[:subset_size]
训练模型：train 函数在数据集上训练语义搜索模型，并返回训练过程的消耗时间和内存使用情况。
training_time, training_memory_usage = demo.train(subset_sentences)
进行推理：infer 函数接受一个查询、一组待搜索的句子和返回的结果数量。它对模型进行推理，并返回最匹配的句子、推理过程的消耗时间和内存使用情况。
top_sentences, inference_time, inference_memory_usage = demo.infer(
query, subset_sentences, top=3
)
演示的完整类如下：
class SemanticSearchDemo:
"""A demo class for semantic search using the ScalableSemanticSearch model."""
def __init__(
self,
dataset_path: str,
model: ScalableSemanticSearch,
index_path: Optional[str] = None,
subset_size: Optional[int] = None,
):
self.dataset_path = dataset_path
self.model = model
self.index_path = index_path
self.subset_size = subset_size
if self.index_path is not None and os.path.exists(self.index_path):
self.loading_time = self.model.timed_load_index(self.index_path)
else:
self.train()
def load_data(self, file_name: str) -> List[str]:
"""Load data from a file.
Args:
file_name: The name of the file containing the data.
Returns:
A list of sentences loaded from the file.
"""
with open(f"{self.dataset_path}/{file_name}", "r") as f:
reader = csv.reader(f, delimiter="\t")
next(reader) # Skip the header
sentences = [row[3] for row in reader] # Extract the sentences
return sentences
def train(self, data: Optional[List[str]] = None) -> Tuple[float, float]:
"""Train the semantic search model and measure time and memory usage.
Args:
data: A list of sentences to train the model on. If not provided, the data is loaded from file.
Returns:
A tuple containing the elapsed time in seconds and the memory usage in megabytes.
"""
if data is None:
file_name = "GenericsKB-Best.tsv"
data = self.load_data(file_name)
if self.subset_size is not None:
data = data[: self.subset_size]
elapsed_time, memory_usage = self.model.timed_train(data)
if self.index_path is not None:
self.model.save_index(self.index_path)
return elapsed_time, memory_usage
def infer(
self, query: str, data: List[str], top: int
) -> Tuple[List[str], float, float]:
"""Perform inference on the semantic search model and measure time and memory usage.
Args:
query: The input query to search for.
data: A list of sentences to search in.
top: The number of top results to return.
Returns:
A tuple containing the list of top sentences that match the input query, elapsed time in seconds, and memory usage in megabytes.
"""
elapsed_time, memory_usage = self.model.timed_infer(query, top)
top_indices, _ = self.model.search(query, top)
index_map = self.model.index_to_sentence_map(data)
top_sentences = self.model.get_top_sentences(index_map, top_indices)
return top_sentences, elapsed_time, memory_usage
我们可扩展语义搜索引擎的性能评估
为了评估我们可扩展的语义搜索引擎的性能，我们可以测量各种操作的时间和内存使用情况，如训练、推理和加载索引。ScalableSemanticSearch 类提供 timed_train、timed_infer 和 timed_load_index 方法来测量这些基准。
train_time, train_memory = semantic_search.timed_train(corpus)
infer_time, infer_memory = semantic_search.timed_infer(query, top)
关于执行时间和内存使用情况的训练和推理性能图表如下。我们将讨论和解释这些结果，基于我们使用的语料库大小选择的算法。
图 2：不同数据集大小的训练时间和内存使用情况
图 3：不同数据集大小的推断时间和内存使用情况
使用 L2 距离的精确搜索
使用 L2 距离的精确搜索是一种穷举搜索方法，它执行线性扫描以找到最近邻。
理论复杂度
-
时间复杂度：O(n) — 因为它需要将查询向量与数据集中每个向量进行比较。
-
内存复杂度：O(n) — 它存储数据集中的所有向量。
观察到的复杂度
- 从上述图表中可以观察到，对于少于 1500 个文档，训练时间和内存使用量都随着文档数量线性增加，这与预期的理论复杂度相符。
使用产品量化和 L2 距离的近似搜索
使用产品量化和 L2 距离的近似搜索方法是一种近似最近邻搜索方法，它采用产品量化和倒排文件结构来提高效率。聚类数量（k）是此方法中的一个重要因素，计算公式为：max(2, min(n_data_points, int(np.sqrt(n_data_points))))。
简单来说，这个公式确保了：
-
至少有 2 个聚类，提供了最低水平的分区。
-
聚类的数量不会超过数据点的数量。
-
作为启发式方法，使用数据点数量的平方根来平衡搜索精度和计算效率。
理论复杂度
-
时间复杂度（训练）：O(n * k) — k-means 聚类算法在训练阶段的复杂度。
-
内存复杂度（训练）：O(n + k) — 它存储聚类的质心和残差码。
-
时间复杂度（推断）：O(k + m) — 其中 m 是要搜索的最近聚类的数量。由于层次结构和近似，它比线性搜索更快。
-
内存复杂度（推断）：O(n + k) — 需要存储倒排文件和质心。
观察到的复杂度
从上述图表中可以观察到，对于超过 1500 个文档：
-
训练时间复杂度：增长速度快于线性，这与预期的理论复杂度 O(n * k) 相符，因为 k 随着数据点（n）的增加而增长。
-
训练内存复杂度：内存使用量随着文档数量的增加呈非线性增长，这与预期的理论复杂度 O(n + k) 相符。
-
推断时间复杂度：执行时间几乎保持不变，与期望的理论复杂度 O(k + m) 一致，因为 m 通常远小于 n。
-
推断内存复杂度：内存使用量随着文档数量线性增加，这与预期的理论复杂度 O(n + k) 相符。
我们还可以通过将搜索引擎的前几个结果与手动整理的真实结果集进行比较，来评估搜索引擎的准确性和召回率。我们可以通过迭代各种查询并比较结果来计算整个数据集的平均准确性和召回率。
结论
展示的方法突显了使用 FAISS 和 Sentence Transformers 进行语义搜索的可扩展性，同时揭示了改进机会。例如，集成先进的 Transformer 模型来对句子进行编码，或测试替代的 FAISS 配置，可能会加快搜索过程。此外，研究最先进的模型，如 GPT-4 或 BERT 变体，可能会提高语义搜索任务的性能和准确性。
可扩展的语义搜索引擎的几个潜在应用包括：
-
在广泛的知识库中检索文档
-
在自动化系统中回答问题
-
提供个性化推荐
-
生成聊天机器人回应
利用 FAISS 和 Sentence Transformers，我们开发了一种可扩展的语义搜索引擎，能够高效处理数十亿份文档并提供准确的搜索结果。这种创新的方法可能会显著影响语义搜索的未来以及其在各个行业和应用中的影响。
随着数字数据的增长，对高效且准确的语义搜索引擎的需求变得越来越重要。基于 FAISS 和 Sentence Transformers，可扩展的语义搜索引擎为克服这些挑战奠定了坚实的基础，并革新了我们搜索和获取相关信息的方式。
未来的发展涉及整合更先进的自然语言处理和机器学习技术，以增强搜索引擎的能力。这些改进可能包括用于更好理解上下文、意图以及查询词和短语之间关系的无监督学习方法，以及处理语言使用中的歧义和变体的技术。
关于我
连续创业者和 AI 领域的领军人物。我为企业开发 AI 产品，并投资于 AI 相关的初创公司。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
![[source]](https://en.wikipedia.org/wiki/Shortest_path_problem#/media/File:Shortest_path_with_direct_weights.svg){kind=link}
