如何用Python编写更好的科学代码?

2022-02-20 18:33:26

任何科学努力的很大一部分都在于编写代码。无论是典型的机器学习建模、分析,还是参与数据项目,大部分时间都花在了新功能的原型设计上。由于是探索性的,预计该计划的几个部分将被替换或修改,通常会超出最初的计划。

与“消费者软件”不同,变化通常不是由客户的需求引起的,而是由过程的结果引起的。因此,如果实验证据表明存在不同的路径,以不需要“完全重建”的方式设计它是非常有价值的。

编写科学代码有两个(额外的)具体挑战:第一个与数学错误有关。计算中的错误通常很难追踪,尤其是当代码语义正确时。没有发现漏洞。没有提出任何例外。一切看起来都不错,但(数字)结果是错误的。特别是,通过实施概率模型,结果有时可能看起来不错,这取决于一些初始条件或随机因素。

第二个来自前面描述的事实。总会有实验性的部分,而且它们一直在变化。因此,关键是要设计每一个部件,使大部分工作能够为下一阶段的发展奠定坚实的基础。

在本文中,我们将重点介绍一些模式,这些模式可以使代码更加健壮、易于理解,并且总体上更易于处理。您将看到简单的改进可以带来更多可重用的代码和更少的bug。

为了演示,我们的任务是计算随机过程结果的预期值。从数学上讲,可以归结为以下公式:

正如你可能想知道的,这里的挑战是,你希望它能与任何分布:连续或离散。或者,如果不合理,至少要认识到问题的本质,这样才能顺利地调整代码。

让我们从一个不太好的代码开始。假设你想掷六面骰子。由于每个结果的概率都相等,因此需要计算抽样结果的平均值。

导入随机导入统计数据def die(sides=6):返回随机。randint(1,sides)def expected_value(n_samples):samples=[范围内(n_samples)的die()]返回统计信息。平均值(样本)

首先,die函数一次返回一个样本。需要调用N次才能获得N个样本,这很慢。

其次,期望_值函数强烈依赖于生成样本的模函数。显然,一旦你考虑使用不同的骰子,比如说12面骰子。在这种设计中,您需要“打开”预期的_值,以接受额外的参数,然后将其传递到模具,以扩展到更一般的情况。虽然这会起作用,但它使得ExpEdTyx值的接口不直观,但是该解决方案仍然依赖于使用骰子作为采样源,从而使得很难考虑其他分布。

这是非常明显的,但你只是把问题转移到其他地方……现在,样本变成了一个新的实体,用于存储数据(甚至非常大),而且它是相当匿名的。expected_value函数希望接收它,但准备它的工作由您自己负责。

另一种方法是将模具保持在内部,将其作为对象传递给预期的_值。

从functools导入部分十二面_die=partial(die,12)def预期_值(die,n_samples):samples=[die()用于范围内(n_samples)]返回统计信息。平均值(样本数)ev=预期值(十二面,10000)

该想法使用了一个准备好的“版本”模具,并使预期的_值将其用作样本来源。然而,出现了一个新问题:预期的_值仅与die兼容。它不能用任何其他“样本生成器”计算结果,或者至少不能保证它能正确地计算结果。

第三个想法是在更抽象的层面上认识问题,并设计更好的接口。

存在一个概率分布,我们可以从中取样。(它可以是一个骰子,一枚硬币,正态分布——没关系)。

有一种数学运算可以消耗和转换数据。(例如,计算平均值、方差等)。

让我们更加关注如何构建正确的抽象,以及如何控制它们的相互兼容性。

从数学上讲,概率分布可以是函数——连续或离散,有限或无限,我们可以从中抽取样本。根据问题的不同,此功能的“处方”可能会非常不同。我们可以使用“现有”公式,如高斯分布或泊松分布,但它也可以是从直方图衍生的“自定义”创建。

由于@abstractmethod的存在,我们的发行版强制我们对从这个抽象派生的任何子类实现sample方法。对于我们的死亡,这可以是:

将numpy作为np类模具导入(分发):def__init__(self,sides=6):self。侧面=侧面def样本(自身,n_样本):返回np。随机的randint(1,自边+1,大小=n_个样本)

这里,通过调用特定于抛出公平骰子的方法sample来交付样本:骰子(12)。样本(10000)。此外,得益于numpy,我们可以通过使用np替换列表理解来快速创建大量数据。恩达雷。

事实上,事情还可以进一步改善。目前,调用Die()会返回如下内容<__主要的。模具位于0x7F43F448400>;,这并不具有信息性。对于Pythonthy,它们是同一类的两个不同对象实例。要解决这个问题,我们还需要实现另外两种方法:

def__repr__(self):返回f";{self._class__._name__}(sides={self.sides})和#34;def__eq__(self,other):如果是instance(other,self._class___):返回self。两边。双方返回错误

__repr______________________________________________。

每次实现这四种方法可能会很乏味。此外,当前的Die实现并不能通过将属性指定给像这个Die这样的现有对象来防止对象的更改,即使是意外更改。侧面=20。

从dataclasses导入dataclass@dataclass(冻结=True)类死亡(分布):侧面:int=6 def样本(self,n):返回np。随机的randint(1,自边+1,大小=n)

这个例子的行为和之前的一样。此外,设置freezed=True,为die指定一个新值。双方将提出一个例外。如果我们想要一个新的模具,我们应该创建一个新的对象。

现在,我们的expected_value函数可能会将die作为一个分布对象,并通过调用它的sample方法进行计算。

上面的例子很简洁。我们确切地知道预期_值的作用,并且很容易测试。然而,n面模具并不是我们可能想要计算的期望值的唯一分布。例如,抛硬币的结果不是数字的(除非我们建立一个惯例并坚持下去)。当然,提供一些关于哪些接口可以一起使用以及如何使用的提示是有意义的。

对于python这样的动态类型化语言,您不必拘泥于变量的类型。然而,使用各种IDE和工具(如mypy),键入可以帮助您发现潜在的故障点,并使代码更加透明。

输入import Generic、Sequence、TypeVar D=TypeVar(";D";)类分布(ABC,泛型[D]):@abstractmethod def sample(self,n:int)->;序列[D]:@数据类(冻结=True)类模具(分布[int]):侧面:int=6 def样本(自身,n:int)->;序列[int]:返回np。随机的randint(1,self.sides+1,size=n)@dataclass(冻结=True)类硬币(分布[str]):结果:元组=(";H";,";T";)公平性:浮动=0.5 def样本(自身,n:int)->;Sequence[str]:p=(self.fairity,1.0-self.fairity)返回np。随机的选择(self.outcocts,size=n,p=p)@dataclass(冻结=True)class Gaussian(分布[float]):mu:float=0.0 sigma:float=1.0 def样本(self,n:int)->;序列[float]:np。随机的正常(loc=self.mu,scale=self.sigma,size=n)

这里发生了几件事。多亏了D=TypeVar(";D";),我们现在可以定义一个新的变量类型,通过它我们可以参数化每个分布的类型。您可以注意到,Distribution不仅在抽象基类之后继承,而且在泛型[D]之后继承,这也将它转换为一个新的(参数化的)类型。现在,它变成了一种身份,构成了一种新的数据类型。

每个版本的示例都会返回一个特定类型的序列,该序列对每个发行版的上下文都有意义。这样,我们就有了一个统一的接口,它也是一个参数化的接口。我们可以使用它来确保预期_值的正确行为:

当将例如die=die()或gaussian=gaussian()传递到预期的_值时(因为int和float都是数字),传递coin=coin()将被例如mypy标记出来,说明

正如您所见,使用类型化设计接口有助于将意图形式化,并尽早捕获错误。您甚至可以利用numpy的数据类型将其提升到下一个级别。这样,您不仅可以确保不同的元素组合在一起,还可以更加注意数据的内存占用。

进口numpy。键入npt类分布(ABC,泛型[D]):@abstractmethod def sample(self,n:int)->;NPNDArray[np.泛型]:。。。级模具(分配[int]):侧面:int=6 def样本(自身,n:int)->;《不扩散条约》。NDArray[np.uint8]:返回np。随机的randint(1,self.sides+1,size=n,dtype=np.uint8)

这样,你甚至会被告知,如果死亡。sample方法返回与严格无符号8位整数不同的数字。问题是你是否想深入到那一步?这是一件值得思考的事情。

让我们回到设计计算部分。到目前为止,我们已经预期了_值,它可以用于数值分布。当然,我们可以计算骰子和高斯的期望值,但不能计算硬币的期望值。与当前的设计不同。

我们可以通过映射来创建代理分布,例如(";H";,";T";)->;(0,1),或

第一种方法创建了一个人造身体,其想法依赖于惯例。它不会阻止任何人用(";H";,";t";)定义另一个代理->;(1,0),导致错误难以检测。

def预期_值(d:分布[d],f:可调用[[d],任意]=lambda x:x,n:int=1000)->;float:返回np。平均值(np.沿_轴应用_(f,轴=0,arr=d.样本(n)))

expected_value的第二个参数是一个可调用的(函数),我们可以选择使用它来转换Coin()分布的结果。然而,默认情况下,它将保持结果不变。

die=die(12)预期_值(die,n=10000)gaussian=gaussian(mu=4.0,sigma=2.0)预期_值(gaussian,n=100000)#但coin=coin(公平性=0.75)预期_值(coin,f=lambda x:np.其中(x==";H";,1.0,0.0))

在这里,我们不仅避免创建代理分布,而且还设法避免将预期的_值转换为任何特定的数据转换方式。expected_value函数只执行它承诺的操作:计算期望值。如果需要任何调整或转换,则由外部提供。请注意,这里我们还有一个选项:我们可以定义一个命名函数(例如coin_转换),以防我们计划重用它,或者在单独的定义没有增加价值时坚持使用lambda。

事实证明,抽象数学计算本身非常有用,尤其是在设计迭代算法时。通常,除了主要的计算之外,我们还必须关注一些附带的结果,比如收敛、提前停止(最大迭代)、度量等等。

让我们以常数为例。从数学上讲,我们可以通过以下极限得到它的值:

定义近似值(初始值:浮点,最大值:int=10,ε:浮点=0.01)->;浮动:e_值=范围内n的初始值(1,最大值+1):新的e_值=(1.0+1.0/n)**n如果abs(新的e_值-e_值)<;epsilon:return new_e_value e_value=new_e_value return new_e_value

首先,函数做三件事而不是一件事。第8行。是计算的绝对本质。然而,由于早期的停止和收敛条件,我们留下了大量的代码开销,这与实际计算紧密耦合。虽然这两个条件看起来更一般,但如果我们选择替换近似的主题(例如平方根),我们将不得不复制粘贴这段附加代码,并确保它不会破坏新算法。

第二,关于参数化这两个条件,我们唯一的选择是对max_iter和epsilon的值进行硬编码,或者允许用户提供它们作为参数。它破坏了界面,使测试更加困难。

最后,该算法“急切地”生成数据。它没有把重点放在数学上,并在“被要求时”提供值,而是把数据扔给你。对于大量数据,这可能会导致内存问题。

现在,让我们通过分离不同部分之间的责任来同时解决这三个问题。我们有三件事:

输入import Iterator import itertools def approximate_e()->;迭代器[float]:n=1,而True:yield(1.0+1.0/n)**n+=1 def early_stop(值:迭代器[float],最大值:int=100)->;迭代器[float]:返回itertools。islice(值,最大值)定义收敛(值:迭代器[float],ε:float=0.01)>;迭代器[float]:对于itertools中的a,b。成对(值):如果(a-b)<;伊普西隆:休息

这个设计使用迭代器,它实现“惰性”加载。数据项仅在请求时逐个返回(因此使用关键字yield)。多亏了这一点,我们(几乎)不必担心记忆力。

此外,这三个功能中的每一个都可以单独存在。它们有特定的接口,可以单独进行单元测试。

值=近似值(e()值=提前停止(值,最大值=50)值=收敛(值,ε=0.001)值中值:打印(";e变成:";,值)

def pairwise(值:迭代器[D])->;迭代器[Tuple[D,D]]:a=next(values,None)如果a为None:返回值中的b:yield a,ba=b

科学编程带来了额外的挑战。这是一门令人兴奋的学科,但陷阱的数量似乎随着问题的复杂性至少呈二次增长。

在本文中,我们讨论了在每项科学编程任务中似乎反复出现的两个主要组成部分:数据和计算。希望通过以正确的方式对它们进行抽象,您不仅可以提高代码的效率,减少错误,还可以让编码成为更好的体验。为了你和你的团队!