0

* 2010 年 6 月 17 日编辑

我正在尝试了解如何改进我的代码(使其更 Pythonic)。此外,我有兴趣编写更直观的“条件”来描述生物化学中常见的场景。下面程序中的条件标准我在答案#2 中解释过,但我对代码不满意——它工作得很好,但并不明显,而且对于更复杂的条件场景也不容易实现。欢迎提出想法。欢迎评论/批评。第一次发帖经验@stackoverflow-如果需要,请评论礼仪。

该代码会生成一个值列表,这些值是以下练习的解决方案:

“在您选择的编程语言中,实施 Gillespie 的第一反应算法来研究反应 A--->B 的时间行为,其中从 A 到 B 的转变只有在存在另一种化合物 C 的情况下才会发生,并且其中 C 与 D 动态相互转换,如下面的 Petri 网模型所示。假设在反应开始时有 100 个 A 分子,1 个 C 分子,并且没有 B 或 D 存在。将 kAB 设置为 0.1 s-1 和kCD 和 kDC 均达到 1.0 s-1。模拟系统在 100 s 内的行为。

def sim():
    # Set the rate constants for all transitions
    kAB = 0.1
    kCD = 1.0
    kDC = 1.0

    # Set up the initial state
    A = 100
    B = 0
    C = 1
    D = 0

    # Set the start and end times
    t = 0.0
    tEnd = 100.0

    print "Time\t", "Transition\t", "A\t", "B\t", "C\t", "D"

    # Compute the first interval
    transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)
    # Loop until the end time is exceded or no transition can fire any more
    while t <= tEnd and transition >= 0:
        print t, '\t', transition, '\t', A, '\t', B, '\t', C, '\t', D
        t += interval
        if transition == 0:
            A -= 1
            B += 1
        if transition == 1:
            C -= 1
            D += 1
        if transition == 2:
            C += 1
            D -= 1
        transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)


def transitionData(A, B, C, D, kAB, kCD, kDC):
    """ Returns nTransition, the number of the firing transition (0: A->B,
    1: C->D, 2: D->C), and interval, the interval between the time of
    the previous transition and that of the current one. """
    RAB = kAB * A * C
    RCD = kCD * C
    RDC = kDC * D
    dt = [-1.0, -1.0, -1.0]
    if RAB > 0.0:
        dt[0] = -math.log(1.0 - random.random())/RAB
    if RCD > 0.0:
        dt[1] = -math.log(1.0 - random.random())/RCD
    if RDC > 0.0:
        dt[2] = -math.log(1.0 - random.random())/RDC
    interval = 1e36
    transition = -1
    for n in range(len(dt)):
        if dt[n] > 0.0 and dt[n] < interval:
            interval = dt[n]
            transition = n
    return transition, interval       


if __name__ == '__main__':
    sim()
4

4 回答 4

1

有关化学 rxns 的简单随机模拟背后的数学信息:

通常,像这样的过程被模拟为离散事件,每个事件发生的概率为“P”,给定特定的速率常数“k”和时间间隔“dt”内的许多可能事件“n”:P=1-e **(-k dtn)。在这里,我们忽略了每个事件的实际时间(~0),而是关注事件发生的时间间隔。任何熟悉 N 选择 K 问题/伯努利试验的人都会理解 1/e 的存在,例如,当 N=K 且 N->oo 时,不从 N 中选择特定元素的概率接近 1/e。因此,在随机化学反应(一阶)中,分子不会发生反应(未被选择)的概率是 1/e 的某个幂...该幂取决于时间间隔和速率常数以及分子数和速率常数有问题。相反,1-(1/e)^xyz 给出了任何特定分子将发生反应(被选择)的概率。

在模拟方面,将我们的总时间间隔分成更小的间隔并使用随机数生成器来预测事件是否在给定的时间间隔内发生是合乎逻辑的——例如,如果我们将单个的 dt 分成更小的 10 个间隔,0 到 0.1 之间的数字表示发生了事件,而 0.1 到 1.0 之间的数字表示没有发生。然而,事件发生的确切时间存在不确定性——因此我们必须缩短间隔——这很快就会成为一场失败的战斗,因为这种方法仍然存在不确定性。

这个问题的解决方法是取上式两边的自然对数(这里是'ln',py中默认是log()),求解dt,得到dt= (-ln(1-P)) /(k*n)。然后随机生成概率 P,为每个事件给出一个确定的 dt。

于 2010-06-17T23:46:33.413 回答
1

不确定你是否看过这个。

http://stompy.sourceforge.net/html/userguide_doc.html

我从事类似的工作,最近碰巧发现了这一点。

于 2011-03-04T06:10:16.383 回答
0

****编辑****我最初解释错了!!!!以下是正确的——贾斯汀,这个程序使用一个聪明的标准来“加权”每个事件。RAB、RCD 和 RDC 值都通过将 kAB、kCD 和 kDC 乘以 C 或 D(在这种情况下可以是 1 或 0)来给出真/假参数。D 的值为零,因此 RDC 将阻止 dt[2] 在

对于范围内的 n(len(dt)):如果 dt[n] > 0.0 并且 dt[n] < 间隔:

陈述。此外,以下——

    if transition == 1:
        C -= 1
        D += 1
    if transition == 2:
        C += 1
        D -= 1

规定当事件 C->D 发生时(转换 1),下一个事件必然是 D->C(转换 2),因为在 dt[] 中的三个值中,只有 dt[1] 不为零,因此满足上述标准。那么,我们如何衡量转换 0 或转换 1 发生的可能性呢?这有点棘手,但在以下几行中是固有的:

interval = 1e36
transition = -1
for n in range(len(dt)):  
    if dt[n] > 0.0 and dt[n] < interval:            
        interval = dt[n]            
        transition = n            
return transition, interval   

"for n in range (len(dt)):" 返回列表 dt[] 的所有值。下一行指定必须满足的条件,即每个值必须大于 0 且小于间隔。对于过渡 0,区间是 1e36(应该代表无穷大)。问题是间隔然后设置为过渡 0,因此对于 dt[] 中的第二个值,过渡 1,标准规定它必须小于过渡 0 的 dt 才能发生......或者换句话说它必须发生得更快才能发生,这与化学逻辑一致。我最担心的是,“t += 间隔”行所规定的累积 t 值可能并不完全公平……因为由于 t1 触发独立于 t0 触发,因此 t0 触发和占用 0.1 秒不应该排除 t1 使用相同的 . 1 秒开火......但我正在努力解决这个问题......欢迎提出建议!这是脚本的详细打印输出,包括转换 1 和 2 的触发:

时间转换 ABCD

dt0= 0.0350693547214
dt1= 2.26710773787
interval= 1e+36
dt= 0.0350693547214
transition= 0

0.0 0 100 0 1 0

dt0= 0.000339596342313
dt1= 0.21083283004
interval= 1e+36
dt= 0.000339596342313 
transition= 0

0.0350693547214 0 99 1 1 0

dt0= 0.0510125874767
dt1= 1.26127048627
interval= 1e+36 
dt= 0.0510125874767
transition= 0

0.0354089510637 0 98 2 1 0

dt0= 0.0809691957218
dt1= 0.593246425076
interval= 1e+36
dt= 0.0809691957218
transition= 0

0.0864215385404 0 97 3 1 0

dt0= 0.00205040633531
dt1= 1.70623338677
interval= 1e+36
dt= 0.00205040633531
transition= 0

0.167390734262 0 96 4 1 0

dt0= 0.106140534256
dt1= 0.0915160378053
interval= 1e+36
dt= 0.106140534256
transition= 0
interval= 0.106140534256
dt= 0.0915160378053
transition= 1

0.169441140598 1 95 5 1 0

dt2= 0.0482892532952
interval= 1e+36
dt= 0.0482892532952
transition= 2

0.260957178403 2 95 5 0 1

dt0= 0.112545351421
dt1= 1.84936696832
interval= 1e+36
dt= 0.112545351421
transition= 0

0.309246431698 0 95 5 1 0

贾斯汀,我不确定你所说的 dt[2] 小于 1e36 使其在过渡 2 中“停留”是什么意思?这不会发生,因为

if transition == 2:
            C += 1
            D -= 1

陈述。任何人都知道实现此目的的更直接方法

哈哈,让燃烧开始吧——你们太棒了——我真的很感谢你的反馈!Stackoverflow 是非常合法的。

于 2010-06-17T23:37:55.723 回答
0

我不知道 Gillespie 算法,但我假设您已经检查过程序是否收敛到正确的平衡。因此,我将您的问题解释为

“这是一个有效的物理程序,我怎样才能让它更 Pythonic”

做类似以下的事情可能会更pythonic

R = [ kAB * A * C,  kCD * C, kAB * A * C]
dt = [(-math.log(1-random.random())/x,i) for i,x in enumerate(R) if x > 0]
if dt:
     interval,transition = min(dt)
else:
     transition = None

如果你想在物理学中使用 python,那么我建议你学习 numpy。因为 numpy 对于许多问题来说更快。所以这里是一个 numpy 解决方案的一些未经测试的部分。将以下内容添加到程序的标题中

from numpy import log, array, isinf, zeros
from numpy.random import rand

然后你可以用类似下面的东西替换里面的TransitionData

R = array([ kAB * A * C,  kCD * C, kAB * A * C])
dt = -log(1-rand(3))/R
transition = dt.argmin()
interval = dt[transition]
if isinf(interval):
    transition = None

我不知道引发 StopIteration 异常而不是返回 None 是否更符合 Python 风格,但这是一个细节。

您还应该将浓度存储在单个数据结构中。如果你使用numpy,那么我建议你使用数组。同样 yoy 可以使用数组 dABCD 来存储浓度的变化(您可能会想出更好的变量名称)。在循环外添加以下代码

ABCD = array([A,B,C,D])

dABCD = zeros(3,4)
dABCD[0,0] = -1#A decreases when transition == 0
dABCD[0,1] = 1 #B increases when transition == 0
dABCD[1,2] = -1#C decreases when transition == 1
dABCD[1,3] = 1 #D increases when transition == 1
.....      etc

现在您可以将主循环替换为以下内容

while t <= tEnd:
       print t, '\t', transition, '\t', ABCD
       transition, interval = transitionData(ABCD, kAB, kCD, kDC)
       if transition != None:
            t += interval
           ABCD += dABCD[transition,:]
       else:
           break;
else:
       print "Warning: Stopping due to timeout. The system did not equilibrate"

可能还有更多工作要做。例如,dABCD 应该是一个稀疏数组,但我希望这些想法可以作为一个开始。

于 2010-06-17T10:30:56.640 回答