首页 > 编程语言 >基于Python程序模拟核酸检测寻找最优化方案

基于Python程序模拟核酸检测寻找最优化方案

时间：2023-04-17 19:33:10浏览次数：44  
标签：方案 group Python 检测 range 核酸 分组 print 最优化

本文中的数学建模问题来源于NKU的数学建模第二次实战演练，由于本次是我来进行程序的编写，故将代码与笔记记录在这里。

问题提要

现有800万市民报名参与核酸检测，如果对每人逐一进行检测，所需时间和检测能力都超过现实情况，所以拟采用混样检测(group testing)方式进行。先考虑混样规模为\(k\)人份一组（一份的咽拭子样本不宜过多，限定\(k\leq10\)）。\(k\)个人的样本混在一起进行检测，如果是阴性，那就全部通过，如果是阳性，可以将该组再进一步分组，或者逐一进行检测。假设人群的阳性率为\(\alpha\) (不超过10%)，每个人最多采集\(n\) 个样本（限定\(n\)不超过3，太多，对于被测试者不友好，也不宜保存和管理）。在满足上述条件的情况下，请设计一个方案（可以是理论证明和公式，也可以是模拟程序），提供给核酸检测团队，帮助核酸检测人员尽量高效的完成大规模（假设固定为800万）核酸检测。

问题分析

根据问题的实际情况，考虑利用Phython程序来进行模拟。为了计算结果的速度，我们考虑将样本压缩到80000人，当然接下来的思路实际上可以套用到任意总数上，只不过到了800000程序运行时间就肉眼可见的慢了。
下面我们考虑进行多次分组。

第一次分组

根据问题条件，首次混样分组以\(k\)人份为一组，且限定\(k\leq10\).则此时可以知道第一次分组混检次数为\(N_1=\lceil \frac{8\times10^7}{k} \rceil\).第一次分组混检次数与分组在\(k\)确定时已经确定。

第二次分组

下面考虑第二次分组，第一次混检结果为阴性的组此时已无需考虑，结果为阳性的组仍需进一步检测以确定阳性病例。为使问题研究简便，下面我们考虑单独一个人数为\(k\)的阳性小组在不同分组情形下的检测次数。

下面我们将第二次分组的问题简化抽象为以下问题：

  1. 对于4人小组，我们可将该小组拆分为\(1+1+1+1,1+2+1,2+2,1+3\)的正整数和的形式，针对本问题不考虑其本身，对于任意的正整数\(k\),如何罗列出\(k\)人的小组拆分成多个正整数和的方案？

  2. 我们将1.中拆分之后的小组称为细分小组，如何计算或模拟不同细分方案下核酸检测所需次数？

针对以上问题，考虑到实际细分小组中的阳性病例的随机性，考虑利用Python程序模拟，通过产生随机样本，经过多次模拟试验，估测出各方案的可能取值范围，从而帮助我们做出决策。

第三次核酸检测

在经历以上步骤后，第三次核酸检测由于\(n\leq3\)的限制只能进行每人的单检（在仍然不能确定阳性与阴性的前提下）。由此我们仅需考虑两次分组问题。

模型假设

为方便有效地对问题进行研究, 我们不妨作如下假设:

  1. 假设在新冠病毒核酸检测中, 所用的试剂判断准确率为 100%.

  2. 假设在新冠病毒核酸检测中, 所用的试剂同样可以准确地检测出混合样本, 且判断准确率为 100%. 在混合样本均为阴性时, 显阴性特性, 否则显阳性特性。

  3. 假设每个人的取样数足够多, 在样本分组后, 各组之间相互独立, 而且样本分组或稀释不会影响最终检测的结果。

程序建构

生成样本

下面是生成总样本的代码：

def generate_mass(n):
    mass = []
    #alpha = random.uniform(0.01,0.1)
    alpha = 0.05
    for i in range(1,n+1):
        if i <= n*alpha:
            mass.append(1)
        else:
            mass.append(0)
    random.shuffle(mass)    # 洗牌模拟真实情形
    return mass

这个函数在传入总人数\(n\)后生成并返回一个长度为\(n\)的列表，其中内容为整数0，1.0代表实际阴性人员，1表示实际阳性病例。我们利用random.shuffle()方法将生成的列表洗牌打乱以模拟真实情形。如果想要让\(\alpha\)也随机生成，可以选用注释掉的那一句代码。

第一次分组

在生成样本之后我们考虑将样本列表进行第一次分组，即生成一个包含\(N_1\)个长度为\(k\)的列表：

def group1(mass,k,n):
    global N,N1
    N=0
    groups1 = []
    group = []
    for i in range(1,n+1,k):
        for j in range(i,i+k):
            if j<=n:
                group.append(mass[j-1])
            else:
                pass
        groups1.append(group)
        N +=1
        group = []
    #print(N)
    N1=N
    return groups1

注意这里的\(N\)与\(N1\)为程序开始时定义的全局变量，用于存储总检测次数和第一次分组混检的检测次数。返回的groups1即为我们需要的列表。

第二次分组

分组方案生成

第二次分组由于既可单检，也可混检，因此我们需要考虑\(k\)的加法分解。
例如\(k=4\)，我们用列表表示可能的正整数分组就有\([1,1,1,1],[1,1,2],[2,2],[1,3]\)，我们不考虑\([4]\)本身，因为这对第二次检测是没有意义的。
因此我们需要将任意\(k\)的方案输出：

# 考虑找出k的所有正整数加法分解，且每个数最大为k-1

def find_all_partitions(k):
    # 初始化结果列表
    partitions = []
    # 递归函数，用于生成所有的分解
    def find_partitions(n, lst):
        if n == 0:
            # 如果n为0，说明已经找到了一组分解
            lst.sort() # 对分解进行排序
            if lst not in partitions: # 如果分解不在结果列表中，就添加到结果列表
                partitions.append(lst)
        else:
            # 否则从1开始尝试，递归查找所有的分解
            for i in range(1, n + 1):
                find_partitions(n - i, lst + [i])

    # 调用递归函数，生成所有的分解
    find_partitions(k, [])
    # 返回结果列表
    return partitions

利用递归函数我们可以将任意正整数\(k\)的加法分解求出。注意代码中的排序作用是除去仅仅只是顺序不同的方案。我们利用循环从\(k=2\)至\(k=10\)可以输出我们需要的全部方案。

分组方案实施

def group2(groups1):
    
        global N,solutions,N1
        N=N1
        positives = []
        positives2=[]   #需要第三次检测的组
        groups2 = []
        # 先踢掉全阴性的组
        for group in groups1:
            for i in range(0,len(group)):
                if group[i] == 1:
                    positives.append(group)
        #print(positives)

    # 再考虑第二次分组检测，从所有可能分组中进行循环
    
        for positive in positives: 
            num=0
            for i in range(0,len(solution)):
                    # 利用切片取出细分小组
                    groups2.append(positive[num:num+solution[i]])
                    N+=1
                    num += solution[i]
        print('在{0}的情况下第二次分组为'.format(solution))
        #print(groups2)
        for group in groups2:
            for i in range(0,len(group)):
                if group[i] == 1 and len(group)>=2:
                    positives2.append(group)
        #print('这里面需要参加第三次的组为：')
        #print(positives2)
        for i in range(0,len(positives2)):
            for j in range(0,len(positives2[i])):
                N+=1
        print(str(N))

上面的代码将会输出对应分组方案的总检测次数，如果需要多次试验可以将第一个print注释掉（在明确结果对应的分组之后），从而方便我们将数据复制到表格当中进行分析。

结果分析

在运行程序20次试验之后我们将结果汇总并将检测次数最少的八个方案列表如下（总数为\(80000\)，\(\alpha=0.05\)）：

分组方案 平均次数 标准差
[2, 2, 2] 35339 78.56818695
[1, 2, 2] 35661 80.75215167
[1, 1, 1, 2] 35833 64.58018272
[1, 1, 1, 1] 36000 0
[1, 1, 1, 1, 1] 36000 0
[1, 1, 2, 2] 36007 84.90959899
[2, 2, 3] 36038 114.1926442
[1, 2, 3] 36184 106.2002236

根据表中的结果来看分组方案中平均次数较少的方案有如下特点：

  • 第二次分组中少量均分方案检测次数相对较少，其中单检和最大组为2人的分组方案尤其较少。

  • 第一次分组时\(k=5,6\)的情况检测次数相对较少。

  • 第二次分组单检的情况次数稳定，而混检次数可能根据实际情况出现接近上百次的浮动。

我们因此可以得到如下的结论：

  • 第一次分组令\(k=5,6\)对于减少检测次数是有效的，考虑到整除问题选用\(k=5\)较为合理且便于操作。

  • 最具有实操性且检测次数在较少的同时还能保持稳定的方案是第二次分组时直接采用单检方案。

综上我们可以认为\(k=5\)且第二次对所有阳性小组进行单检的方案是最适宜于实际操作的。

注意事项

模型还能从以下角度进行改进：

  • 在试验时可以在\(\alpha = 0.01,0.02,\cdots,0.1\)的条件下分别进行20次试验，可以找到\(\alpha\)不同时分别的最优方案。

  • 代码中仍有\(\frac{n}{k}\)不整除时带来的误差，这源于将分组方案直接套用到不足\(k\)人小组造成的误差。

完整代码

完整代码如下：

# -*- coding: utf-8 -*-
"""
Created on Fri Apr 14 22:54:04 2023
核酸检测模拟程序

@author: xzqbear
"""
import random
N,N1=0,0
solution = []
solutions = []
# 生成总人数列表
def generate_mass(n):
    mass = []
    #alpha = random.uniform(0.01,0.1)
    alpha = 0.05
    for i in range(1,n+1):
        if i <= n*alpha:
            mass.append(1)
        else:
            mass.append(0)
    random.shuffle(mass)    # 洗牌模拟真实情形
    return mass

# 第一次分组
def group1(mass,k,n):
    global N,N1
    N=0
    groups1 = []
    group = []
    for i in range(1,n+1,k):
        for j in range(i,i+k):
            if j<=n:
                group.append(mass[j-1])
            else:
                pass
        groups1.append(group)
        N +=1
        group = []
    #print(N)
    N1=N
    return groups1

'''
在进行第二次分组之前我们需要考虑k人小组可能的分组情形
我们考虑先寻找k这个正整数可能的加法分解情况，例如对于4，可能的分解有：
1+1+1+1，1+2+1，1+3，2+2，4
根据问题，我们可以认为其本身作为加法分解是没有意义的，所以上述情况
符合提议的分解需要去掉4.
下面考虑找出所有可能的分解情况
'''
# 考虑找出k的所有正整数加法分解，且每个数最大为k-1

def find_all_partitions(k):
    # 初始化结果列表
    partitions = []
    
    # 递归函数，用于生成所有的分解
    def find_partitions(n, lst):
        if n == 0:
            # 如果n为0，说明已经找到了一组分解
            lst.sort() # 对分解进行排序
            if lst not in partitions: # 如果分解不在结果列表中，就添加到结果列表
                partitions.append(lst)
        else:
            # 否则从1开始尝试，递归查找所有的分解
            for i in range(1, n + 1):
                find_partitions(n - i, lst + [i])

    # 调用递归函数，生成所有的分解
    find_partitions(k, [])
    
    # 返回结果列表
    return partitions


# 第二次分组

def group2(groups1):
    
        global N,solutions,N1
        N=N1
        positives = []
        positives2=[]   #需要第三次检测的组
        groups2 = []
        # 先踢掉全阴性的组
        for group in groups1:
            for i in range(0,len(group)):
                if group[i] == 1:
                    positives.append(group)
        #print(positives)
    

        
    # 再考虑第二次分组检测，从所有可能分组中进行循环
    
        for positive in positives: 
            num=0
            for i in range(0,len(solution)):
                    # 利用切片取出细分小组
                    groups2.append(positive[num:num+solution[i]])
                    N+=1
                    num += solution[i]
        #print('在{0}的情况下第二次分组为'.format(solution))
        #print(groups2)
        for group in groups2:
            for i in range(0,len(group)):
                if group[i] == 1 and len(group)>=2:
                    positives2.append(group)
        #print('这里面需要参加第三次的组为：')
        #print(positives2)
        for i in range(0,len(positives2)):
            for j in range(0,len(positives2[i])):
                N+=1
        print(str(N))
    

if __name__ == '__main__':

    n = 80000    # 人数
    mass = generate_mass(n)
    #print(mass)
    for k in range(2,11):
        groups1 = group1(mass, k, n)
        #print(groups1)
        solutions = find_all_partitions(k)
        del solutions[-1]
        for solution in solutions:
            group2(groups1)

标签：方案,group,Python,检测,range,核酸,分组,print,最优化
From： https://www.cnblogs.com/han-son-xiong/p/17323416.html

相关文章

  • Python替换文件内容
    文件部分内容如下：<mappingcell="A1">request.aaPriceChangeDesc</mapping><mappingcell="B1">request.aaStartDate</mapping><mappingcell="C1">request.aaSumCode</mapping><mappingcell=......
  • 自动化脚本：一键安装python自定义版本
     1：环境：centos7python2.72：脚本内容：#!/usr/bin/envpythonimportosimportsysimportrequestsimporttarfileimportshutilimportsubprocess#Installnecessarypackagestry:subprocess.check_call(["yum","install","-y&qu......
  • [oeasy]python0132_变量含义_meaning_声明_declaration_赋值_assignment
    变量定义回忆上次内容上次回顾了一下历史python是如何从无到有的看到Guido长期的坚持和努力 编程语言的基础都是变量声明python是如何声明变量的呢？ 变量想要定义变量首先明确什么是变量变量就是数值能变的量英文名称varia......
  • [oeasy]python0132_变量含义_meaning_声明_declaration_赋值_assignment
    变量定义回忆上次内容上次回顾了一下历史python是如何从无到有的看到Guido长期的坚持和努力编程语言的基础都是变量声明python是如何声明变量的呢？变量想要定义变量首先明确什么是变量变量就是数值能变的量英文名称variable计算机在内存中分配出空间用来存储这些能变的量那......
  • random模块&string模块谈python中随机数
    一、概述随机数在程序设计中的属于比较基础的内容，主要用于验证场景（如验证码，生成账号对应的密码等），今天结合random模块和string模块来谈谈python中随机数那些事儿。二、随机数实现相关模块2.1random模块random.random()返回一个随机浮点数。>>>importrandom>>>print(ran......
  • PYTHON学习路径计划图整理
    PYTHON学习路径计划图Python工作环境及基础语法知识了解对于Python基础语法学习部分，学习周期大概为4周，需要的相关资源在网络上都能找到免费的资源，而且质量都不错。相关中文资源如下：1.python工作集成环境包Python(x,y): 下载地址Pycharm： 下载地址2.python数据分析相关库（Pa......
  • 功能不够用？使用C++编写通达信插件及接入Python（二）
    参考：https://zhuanlan.zhihu.com/p/613157262一、准备工作（参考上一篇）安装VS2019 安装pycharm下载 http://help.tdx.com.cn/book.asp《通达信DLL函数编程规范.rar》二、下载python3.x的32位版本，http://www.python.org，随便找个32位版就行了。我准备下载Windowsembeddabl......
  • Python数据分析方向第三方库推荐！
    在Python的应用范围中，数据分析是一个非常不错的就业方向，发展前景好、需求市场大，也成为很多人的首选方向，因此Python也提供了许多数据分析的第三方库。那么Python数据分析第三方库有哪些?这篇文章为大家详细介绍一下。1、PandasPandas是Python强大、灵活的数据分析和探索......
  • 【python-docx】添加表格、合并单元格、设置样式
    1.表格1.1添加表格添加表格很简单，只需要调用一下add_table()即可，返回一个Table对象，参数可以指定行、列、样式fromdocximportDocumentdoc=Document()#添加一个5行3列的表格，样式是网格实线table=doc.add_table(5,3,style="TableGrid")doc.save('./test.docx')1.2添......
  • 盘点Python内置函数sorted()高级用法实战
    今日鸡汤清川带长薄，车马去闲闲。大家好，我是Python进阶者。一、前言前几天在Python钻石交流群有个叫【emerson】的粉丝问了一个Python排序的问题，这里拿出来给大家分享下，一起学习下。其实这里【瑜亮老师】、【布达佩斯的永恒】等人讲了很多，只不过对于基础不太好的小伙伴们来说，还是有......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99