探索性数据分析-Python-中的伽马光谱-第三部分-
探索性数据分析：Python 中的伽马光谱（第三部分）
原文：
towardsdatascience.com/exploratory-data-analysis-gamma-spectroscopy-in-python-part-3/
在这个故事中，我将测试两种方法：
-
我将创建一个公共 Streamlit 应用程序，该应用程序将免费托管在Streamlit Cloud上（应用程序链接已添加到文章末尾）。
-
作为一种更灵活和通用的解决方案，我将创建一个基于 Python HTMX 的应用程序，它可以与真实硬件通信并在实时进行预测。
与前一部分相同，我将使用一个Radiacode闪烁探测器来获取数据（免责声明：本测试中使用的设备由制造商提供；我从中不获得任何商业利润，并且我没有获得关于所有测试的任何编辑意见）。没有 Radiacode 硬件的读者将能够使用 Kaggle 上可用的文件测试应用程序和模型。
让我们开始吧！
1. 同位素分类模型
这个模型在前一部分中进行了描述。它基于 XGBoost，我使用不同的放射性样品训练了这个模型。我使用了可以合法购买的样品，比如复古铀玻璃或 20 世纪 50 年代制造的带有镭表盘的旧手表。正如之前提到的，我还使用了一个 Radiacode 闪烁探测器，这使我能够获取物体的伽马光谱。仅在 10-20 年前，这类探测器只能在大型实验室中找到；如今，它们的价格与中档智能手机相当。
该模型包含三个组件：
-
基于 XGBoost 的模型本身。
-
一系列放射性同位素（如铅-214 或锕-228），模型就是基于这些同位素训练的。Radiacode 闪烁探测器返回 1024 个光谱值，其中 23 个被用于模型。
-
一个标签编码器，用于将列表索引转换为可读名称。
让我们将所有这些封装到一个单一的 Python 类中：
from xgboost import XGBClassifier
from sklearn.preprocessing import LabelEncoder
class IsotopesClassificationModel:
""" Gamma Spectrum Classification Model """
def __init__(self):
""" Load models """
path = self._get_models_path()
self._classifier = self._load_model(path + "/XGBClassifier.json")
self._isotopes = self._load_isotopes(path + "/isotopes.json")
self._labels_encoder = self._load_labels_encoder(path + "/LabelEncoder.npy")
def predict(self, spectrum: Spectrum) -> str:
""" Predict the isotope """
features = SpectrumPreprocessing.convert_to_features(
spectrum, self._isotopes
)
preds = self._classifier.predict([features])
preds = self._labels_encoder.inverse_transform(preds)
return preds[0]
@staticmethod
def _load_model(filename: str) -> XGBClassifier:
""" Load model from file """
bst = XGBClassifier()
bst.load_model(filename)
return bst
@staticmethod
def _load_isotopes(filename: str) -> List:
with open(filename, "r") as f_in:
return json.load(f_in)
@staticmethod
def _load_labels_encoder(filename: str) -> LabelEncoder:
le = LabelEncoder()
le.classes_ = np.load(filename)
return le
@staticmethod
def _get_models_path() -> str:
""" Get path to models. Model files are stored in
'models/V1/' folder """
parent_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
return parent_dir + f"/models/{IsotopesClassificationModel.VERSION}"
一个 Spectrum 类包含我们从辐射探测器获得的光谱数据：
@dataclass
class Spectrum:
""" Radiation spectrum data """
duration: datetime.timedelta
a0: float
a1: float
a2: float
counts: list[int]
在这里，counts 是一个伽马光谱，它由 1024 个通道值组成的列表表示。光谱数据可以使用官方的 Radiacode Android 应用程序导出，或者直接从设备使用 radiacode Python 库检索。
要将光谱加载到模型中，我创建了一个 SpectrumPreprocessing 类：
class SpectrumPreprocessing:
""" Gamma Spectrum Preprocessing """
@staticmethod
def convert_to_features(spectrum: Spectrum, isotopes: List) -> np.array:
""" Convert the spectrum to the list of features for prediction """
sp_norm = SpectrumPreprocessing._normalize(spectrum)
energies = [energy for _, energy in isotopes]
channels = [SpectrumPreprocessing.energy_to_channel(spectrum, energy) for energy in energies]
return np.array([sp_norm.counts[ch] for ch in channels])
@staticmethod
def load_from_xml_file(file_path: str) -> Spectrum:
""" Load spectrum from a Radiacode Android app file """
在这里，我跳过了一些已经在 前一部分 中发布的代码块。伽马光谱的特征提取也在那里进行了说明，我强烈建议首先阅读那部分内容。
现在，让我们测试一下模型！我使用了一个 Radiacode 探测器，在 10 分钟内收集了一个伽马光谱：
辐射探测器，图片由作者提供
这款中国吊坠被宣传为“离子生成”，它有轻微的放射性。使用官方 Radiacode Android 应用程序收集的伽马光谱如下所示：
屏幕截图由作者提供
等待大约 10 分钟后，我将光谱导出为 XML 文件。现在，我们可以运行模型：
from spectrum import SpectrumPreprocessing
from ml_models import IsotopesClassificationModel
sp = SpectrumPreprocessing.load_from_file("spectrum.xml")
model = IsotopesClassificationModel()
result = model.predict(sp)
print(result)
#> Thorium
如我们所见，模型运行良好。我们可以将峰值与已知同位素的光谱（例如，这里 或 这里）进行比较，并确认该光谱属于钍。
2. Streamlit
模型运行良好；然而，我们生活在 21 世纪，几乎没有人会运行控制台应用程序来获取结果。相反，我们可以将应用程序放在网上，这样所有的 Radiacode 用户都可以运行它。
有许多 Python 框架可以用来制作基于浏览器的应用程序，其中 Streamlit 可能是数据科学社区中最受欢迎的。对我们来说重要的是，Streamlit 社区云 平台允许每个人完全免费发布他们的应用程序。为此，让我们首先制作这个应用程序。
2.1 Streamlit 应用程序
Streamlit 框架相对容易使用，至少如果我们需要一个看起来标准的应用程序的话。我个人不是这个方法的粉丝。这些框架将所有底层实现细节隐藏起来。制作原型很简单，但 UI 逻辑将与一个非常狭窄的框架紧密耦合，并且无法在其他地方重用。做所有非标准的事情，这通常不受框架支持，可能几乎不可能实现，或者在没有深入研究大量抽象和代码页的情况下很难实现。然而，在我们的情况下，原型就是我们所需要的。
通常，Streamlit 代码很简单，我们只需要描述我们页面的逻辑层次结构：
import streamlit as st
import logging
logger = logging.getLogger(__name__)
def is_xml_valid(xml_data: str) -> bool:
""" Check if the XML has valid size and data """
return len(xml_data) < 65535 and xml_data.startswith("<?xml")
def get_spectrum(stringio: StringIO) -> Optional[Spectrum]:
""" Load spectrum from the StringIO stream """
xml_data = stringio.read()
if is_xml_valid(xml_data):
return SpectrumPreprocessing.load_from_xml(xml_data)
return None
def main():
""" Main app """
st.set_page_config(page_title="Gamma Spectrum")
st.title("Radiacode Spectrum Detection")
st.text(
"Export the spectrum to XML using the Radiacode app, and "
"upload it to see the results."
)
# File Upload
uploaded_file = st.file_uploader(
"Choose the XML file", type="xml", key="uploader",
)
if uploaded_file is not None:
stringio = StringIO(uploaded_file.getvalue().decode("utf-8"))
if sp := get_spectrum(stringio):
# Prediction
model = IsotopesClassificationModel()
result = model.predict(sp)
logger.info(f"Spectrum prediction: {result}")
# Show result
st.success(f"Prediction Result: {result}")
# Draw
fig = get_spectrum_barchart(sp)
st.pyplot(fig)
if __name__ == "__main__":
logger.setLevel(logging.INFO)
main()
如我们所见，完整的应用程序只需要最少的 Python 代码。Streamlit 将为我们渲染所有的 HTML，包括标题、文件上传和结果。作为额外的好处，我还会使用 Matplotlib 显示一个光谱：
def get_spectrum_barchart(sp: Spectrum) -> plt.Figure:
""" Get Matplotlib's barchart """
counts = SpectrumPreprocessing.get_counts(sp)
energy = [
SpectrumPreprocessing.channel_to_energy(sp, x) for x in range(len(counts))
]
fig, ax = plt.subplots(figsize=(9, 6))
ax.spines["top"].set_color("lightgray")
ax.spines["right"].set_color("lightgray")
# Bars
ax.bar(energy, counts, width=3.0, label="Counts")
# X values
ticks_x = [SpectrumPreprocessing.channel_to_energy(sp, ch) for ch in range(0, len(counts), len(counts) // 20)]
labels_x = [f"{int(ch)}" for ch in ticks_x]
ax.set_xticks(ticks_x, labels=labels_x, rotation=45)
ax.set_xlim(energy[0], energy[-1])
ax.set_ylim(0, None)
ax.set_title("Gamma spectrum")
ax.set_xlabel("Energy, keV")
ax.set_ylabel("Counts")
return fig
现在，我们可以在本地运行应用程序：
streamlit run st-app.py
之后，我们的应用程序完全可用，可以在浏览器中进行测试：
屏幕截图由作者提供
如前所述，我不是非常高级框架的粉丝，更喜欢更好地理解事物“底层”的工作原理。然而，考虑到我只用了大约 100 行代码就制作了一个功能齐全的 Web 应用程序，我无法抱怨——对于原型设计，它工作得很好。
2.2 Streamlit 社区云
当应用程序在本地进行测试时，是时候使其公开了！Streamlit Cloud 是一个免费服务，显然，它有很多限制：
-
该应用程序在类似 Docker 的容器中运行。您的 GitHub 账户必须与 Streamlit 关联。当容器启动时，它会从 GitHub 拉取您的代码并运行。
-
在撰写本文时，容器资源限制为 2 个核心和高达 2.7GB 的 RAM。运行 70B 大小的 LLM 将过于受限，但对于一个小型的 XGBoost 模型来说，这已经足够了。
-
Streamlit 不提供任何永久存储。在关闭或重启后，所有日志和临时文件都将丢失（如果需要，您可以使用 API 密钥并从 Python 代码连接到任何其他云存储）。
-
在一段时间的非活动状态（大约 30 分钟）后，容器将被停止，所有临时文件也将丢失。如果有人打开应用程序链接，它将再次运行。
如读者所猜，一个不活跃的应用程序几乎不花费 Streamlit 任何费用，因为它只存储了一个小的配置文件。这对于免费服务来说是一个很好的解决方案——它允许我们免费发布应用程序，并提供一个链接供人们运行。
要在 Streamlit 中发布应用程序，我们需要执行三个简单的步骤。
首先，我们需要将我们的 Python 应用程序提交到 GitHub。一个requirements.txt文件也是强制性的。Streamlit 容器使用它来安装所需的 Python 依赖项。在我的情况下，它看起来像这样：
xgboost==3.0.2
scikit-learn==1.6.1
numpy==1.26.4
streamlit==1.47.0
pillow==11.1.0
matplotlib==3.10.3
xmltodict==0.14.2
可以使用.streamlit/config.toml文件更改服务器设置。在我的情况下，我将上传文件大小限制为 1MB，因为所有光谱文件都较小：
[server]
# Max size, in megabytes, for files uploaded with the file_uploader.
# Default: 200
maxUploadSize = 1
第二，我们需要使用 GitHub 账户登录到share.streamlit.io并授予访问源代码的权限。
最后，我们可以创建一个新的 Streamlit 项目。在项目设置中，我们还可以选择所需的 URL 和环境：
图片由作者提供
如果一切操作正确，我们可以看到我们的应用程序正在运行：
图片由作者提供
在这个时候，全世界的用户也可以访问我们的应用程序！在我的情况下，我选择了gammaspectrumdetection这个名字，应用程序可以通过这个链接访问。
3. FastAPI + HTMX 应用程序
正如读者所见，Streamlit 是一个很好的简单原型解决方案。然而，在辐射探测器的案例中，我希望看到来自真实 Radiacode 硬件的数据。这在 Streamlit 中是不可能的；这个库根本不是为了这个目的而设计的。相反，我将使用几个生产级框架：
如前所述，那些没有 Radiacode 硬件的读者将能够使用从真实设备保存的原始日志文件重放数据。文章末尾提供了应用程序和所有文件的链接。
让我们开始吧！
3.1 HTML/HTMX
应用程序连接到 Radiacode 探测器，我决定在页面上显示连接状态、辐射水平和频谱图。在底部，将显示频谱收集时间和机器学习模型预测。
这个布局的index.html文件看起来是这样的：
<!DOCTYPE html>
<head>
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Gamma Spectrum & Monitoring</title>
<script src="https://cdn.jsdelivr.net/npm/[[email protected]](/cdn-cgi/l/email-protection)/dist/htmx.min.js" crossorigin="anonymous"></script>
<script src="https://cdn.jsdelivr.net/npm/chart.js"></script>
<link rel="stylesheet" href="{{ url_for('static', path='styles.css') }}">
</head>
<body>
<script type="text/javascript">
function createChart() {
const ctx = document.getElementById('gammaChart').getContext('2d');
window.chart = new Chart(ctx, {
...
});
}
function updateChartData(labels, data) {
window.chart.data.labels = labels;
window.chart.data.datasets[0].data = data;
window.chart.update();
}
document.addEventListener('DOMContentLoaded', function() {
// console.log('Page loaded');
createChart();
}, false);
</script>
<div class="dashboard">
<div id="reload-status" hx-get="/device/poll" hx-target="#device-info" hx-swap="outerHTML" hx-trigger="load, every 1s"/>
<div class="info-row-top">
<div class="info-box-left">
<span id="device-info">Device: n/a</span><br/>
</div>
<div class="info-box-right">
<span id="device-cps">CPS: n/a</span>
</div>
</div>
<canvas id="gammaChart" class="gamma-chart"></canvas>
<div id="gammaChartData"></div>
<div class="info-row-bottom">
<div class="info-box-left">
<span>Collection Time:</span><br/>
<span id="collection-time">n/a</span>
</div>
<div class="info-box-right">
<span>Prediction:</span><br/>
<span id="prediction-result">n/a</span>
</div>
</div>
<button class="reset-button" hx-post="/device/reset_spectrum" hx-swap="none">Reset Spectrum</button>
</div>
</body>
</html>
在这里，我使用了Chart.js来制作图表，所有的 UI 控件都位于div部分。HTMX 允许我们不用 JavaScript 就能制作一个交互式 Web 应用。这个 HTMX 字符串执行实际的更新：
hx-get="/device/poll" hx-target="#device-info" hx-swap="outerHTML" hx-trigger="load, every 1s"
在这里，HTMX 被配置为每秒调用/device/poll端点，并将结果放入具有#device-info ID 的控件中。
如果需要，读者可以更改 HTML。我不是前端开发者，也许这个 UI 可以做得更有效。然而，页面是工作的，并且完成了它的任务：
作者截图
3.2 FastAPI
如前所述，我决定使用FastAPI作为后端部分。
首先，我们需要创建一个 FastAPI 实例并添加所需的端点：
from fastapi import FastAPI, Request, Response
from fastapi.responses import FileResponse
from contextlib import asynccontextmanager
@asynccontextmanager
async def lifespan(app: FastAPI):
logging.info("FastAPI::Started")
app.state.device = RadiaCodeDevice()
app.state.device.start_in_background()
yield
app.state.device.stop()
logging.info("FastAPI::Ended")
app = FastAPI(lifespan=lifespan)
app.mount("/static", StaticFiles(directory="static"), name="static")
@app.get("/device/poll")
async def device_poll():
""" Get device data """
# The code is placed below
@app.post("/device/reset_spectrum")
async def device_reset_spectrum():
""" Reset gamma spectrum data """
app.state.device.reset_spectrum()
return "OK"
@app.get("/")
def read_index(request: Request):
context = {"request": request}
return templates.TemplateResponse("index.html", context)
@app.get("/favicon.ico", include_in_schema=False)
async def favicon():
return FileResponse("static/favicon.ico")
if __name__ == '__main__':
import uvicorn
uvicorn.run("app:app", host='0.0.0.0', port=8000, reload=True)
正如读者所见，我使用了一个lifespan回调在 FastAPI 服务器启动时运行 Radiacode 连接。app.state是 FastAPI 中的一个全局状态对象。我们只有一个设备实例，因此使用全局状态是个不错的选择。FastAPI 本身基于 asyncio，Radiacode 连接线程在后台工作。在 asyncio 中，我不能阻塞主线程，所以我使用了两个队列（一个用于数据，一个用于命令）进行数据交换。
在 HTMX 中更新页面的方式很有趣。假设我想更新 CPS（每秒计数）值。作为一个提醒，在 HTML 中，我有这些控件：
<span id="device-info">Device: n/a</span><br/>
<div class="info-box-left">
<span id="device-info">Device: n/a</span><br/>
</div>
<div class="info-box-right">
<span id="device-cps">CPS: n/a</span>
</div>
<div id="reload-status" hx-get="/device/poll" hx-target="#device-info"
hx-swap="outerHTML" hx-trigger="load, every 1s"
/>
在对/device/poll的响应中，我需要为相应的 ID 返回适当的 HTML：
def get_connection_data(self) -> str:
""" Get device data for the HTML poll request """
status = "Connected"
return f'<span id="device-info">{status}</span>'
def get_cps_data(self) -> str:
""" Get CPS data for the HTML poll request """
return f'<span id="device-cps" hx-swap-oob="true">CPS: {self.radiation_cps:.2f}</span>'
def get_spectrum_data(self) -> str:
""" Get spectrum data for the HTML poll request """
return f'<script>updateChartData([], [])</script>'
def get_spectrum_prediction(self) -> str:
...
@app.get("/device/poll")
async def device_poll():
""" Get device data """
poll = app.state.device.get_connection_data()
cps = app.state.device.get_cps_data()
spectrum = app.state.device.get_spectrum_data()
prediction = app.state.device.get_spectrum_prediction()
return Response(
content=poll + cps + prediction + spectrum,
media_type="application/text"
)
HTMX 会自动触发/device/poll端点，并用新数据替换网页上的控件。hx-swap-oob键允许我们在单个请求中更新多个控件。正如我们所见，HTMX 允许我们无缝地将浏览器中的前端页面与 Python 后端连接起来。
get_spectrum_data方法略有不同。HTMX 不支持直接更新 Chart.js 条形图。然而，我可以返回一个将在页面上执行的 JavaScript 代码块。我使用这个特性通过调用放在index.html中的updateChartData方法来更新图表。
完整的源代码更长，这里我只展示了关键部分。文章末尾有一个完整项目的链接。
4. 测试
最后，让我们看看这两个应用的实际运行情况！我的模型可以分类同位素，我想通过一个模型之前未“看到”的物体来测试它。为此，我参观了市中心的一家矿物店：
图片由作者提供
我知道一些矿物可能略微具有放射性，所以我使用辐射探测器选择了最好的（或者根据你的喜好，最差的）一个。它是一种蓝色磷灰石。我在家里进行了测试，我的 GMC 辐射计显示其放射性水平为 0,71 µSv/h，该设备算法将其视为“高”：
图片由作者提供
在这里，“高”并不意味着“危险”——这个值只是对于平均背景辐射水平来说较高，通常约为 0,1 µSv/h。然而，在这种情况下，较高的水平仅在矿物周围 1-2cm 处。它并不危险，也不会造成任何伤害（作为比较，飞机中的辐射，由宇宙射线引起，大约比这高 1.5 倍）。而且，如前一部分所述，盖革计数器可以显示数值，但不能告诉我们为什么这个矿物是放射性的。所以，我将使用 Radiacode 伽马探测器来找出原因。
作为提醒，我首先创建了一个Streamlit 应用并将其发布在 Streamlit 社区云平台上。该应用在线可用，但如果没有人在使用它，应用就处于“睡眠模式”：
截图由作者提供
我使用 Radiacode 安卓应用保存了一个光谱。当 Streamlit 应用启动时（容器启动需要 2-3 分钟），我上传了文件：
分类结果，截图由作者提供
作为第二次测试，让我们运行本地应用。它功能更强大，可以通过 USB 与 Radiacode 探测器通信，也可以作为 UI 来控制硬件。另一个关键的区别是我们可以实时看到探测器中的数据。
我按下了“重置光谱”按钮以重新启动累积，在收集数据 2-3 分钟后，矿物被鉴定为钍：
作者截图
这是一个有趣的发现，因为根据Wikipedia，磷灰石可能含有铀，但那里没有提到钍。然而，根据USGS（美国地质调查局）页面，钍通常可以在磷灰石中找到。在这种情况下，Wikipedia 文章只是不够完整。
结论
在这篇文章中，我测试了检测放射性同位素的机器学习模型。我测试了两种使用模型的方式——一个发布在 Streamlit 社区云中的 Streamlit 应用程序，以及一个功能齐全的 FastAPI 应用程序，它可以从硬件获取数据，并可作为真实辐射探测器的网络界面使用。
如我们所见，该模型运行良好，我能够在蓝色磷灰石中找到微量的钍。这个事实甚至没有在 Wikipedia 中提到，学到新东西总是令人愉快的。然而，正如前一部分所提到的，模型的性能有限。我不是核机构，该模型仅基于可以合法购买的放射性物体进行训练，如带有镭表盘的复古手表。我没有铯或钚等测试源。尽管如此，训练模型并在应用程序中使用它仍然很有趣。可能这样的项目没有商业价值，对同位素分类模型的市场需求为 0.00%。然而，在这样做的时候，我能够使用 FastAPI、HTMX 或 CSS 等工具，这些工具在其他项目中可能很有用。最后但同样重要的是，我玩得很开心，这是最重要的部分 😉
所有数据文件都是使用一个 Radiacode 涡轮探测器收集的。这是一个既好又紧凑的设备，价格与中档智能手机相当，正如我们所见，它允许我们进行一些有趣的实验（免责声明：我从中没有任何利润或其他商业利益）。对于那些没有 Radiacode 硬件的读者，所有收集到的数据都可以在 Kaggle 上免费获取。该存档包含用于训练模型的多种物体的光谱。XML 格式的光谱可以用来测试 Streamlit 应用程序。"replay"文件夹包含从 Radiacode 设备保存的日志，可以用来测试 HTMX 应用程序。
这篇文章的完整源代码可在我的Patreon 页面上找到。这种支持帮助我购买未来类似测试所需的设备和电子产品。读者也可以通过LinkedIn与我联系，我在那里定期发布一些不够完整，不足以成为完整文章的小帖子。
感谢阅读。
浙公网安备 33010602011771号