Python / Biomolecular Physics-试图对表现出条件行为的系统进行简单的随机模拟！

*编辑于2010年6月17日 我正在努力了解如何改进我的代码（使其更加pythonic）。此外，我有兴趣编写更直观的“条件”，描述生物化学中常见的场景。我在答案＃2中解释了以下程序中的条件标准，但我对代码不满意 - 它工作正常，但不是很明显，并且不容易实现更复杂的条件场景。想法欢迎。评论/批评欢迎。首次发布体验@ stackoverflow-如果需要，请评论礼仪。 该代码生成一个值列表，这些值是以下练习的解决方案： “在您选择的编程语言中，实施Gillespie的第一反应算法来研究反应的时间行为A ---> B，其中只有当存在另一种化合物C时才能从A到B的转变，并且其中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 "Timet", "Transitiont", "At", "Bt", "Ct", "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()
    
2019-05-10 5 条评论
分享

没有找到相关结果

    已邀请:

    窝头菊

    不确定你是否看过这个。 http://stompy.sourceforge.net/html/userguide_doc.html 我从事类似的工作，最近我偶然发现了这一点。     
    2019-05-10 0 0
    分享

    春驹晴陪

    关于化学rxns的简单随机模拟背后的数学信息： 通常，这样的过程被模拟为离散事件，每个事件发生概率'P'给定特定速率常数'k'，并且在时间间隔'dt'中有多个可能事件'n'：P = 1-e **（ - KDTN）。在这里，我们忽略了每个事件的实际时间（~0），而是关注事件发生的时间间隔。任何熟悉N选择K问题/ bernouli试验的人都会理解1 / e的存在，例如当N = K且N-> oo时，不从N中选择特定元素的概率接近1 / e。因此，在随机化学反应（一阶）中，分子不会发生反应（未被选择）的概率是1 / e的一些幂...功率取决于时间间隔和速率常数以及有问题的分子数和速率常数。相反，1-（1 / e）^ xyz给出了任何特定分子反应（被选择）的概率。 在模拟方面，将我们的总时间间隔划分为更小的间隔并使用随机数发生器来预测事件是否在给定的时间间隔内发生是合乎逻辑的。如果我们将单个偶数的dt除以10个较小的间隔，则0到0.1之间的数字表示发生了事件，而.1和1.0之间的数字表示它没有。然而，事件发生的确切时间存在不确定性 - 因此我们必须使我们的间隔变小 - 这很快就会变成一场失败的战斗，因为这种方法存在不确定性。 这个问题的解决方案是采用上面等式两边的自然对数（'ln'，默认为log（）），并求解dt，得到dt =（-ln（1-P）） /（K * N）。然后随机生成概率P，为每个事件给出确定的dt。     
    2019-05-10 0 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会更加pythonic，但这是一个细节。 您还应将浓度存储在单个数据结构中。如果你使用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可能应该是一个稀疏数组，但我希望这些想法可以作为一个开始。     
    2019-05-10 0 0
    分享

    肺鬼耙扮群

    ****编辑****我原来解释错了!!!!以下是正确的 - 贾斯汀，这个程序使用一个聪明的标准来“加权”每个事件。通过将kAB，kCD和kDC乘以C或D，RAB，RCD和RDC值都被赋予真/假参数，在这种情况下，它们可以是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且小于interval。对于过渡0，间隔为1e36（应该代表无穷大）。摩擦是然后将该间隔设置为转换0，因此对于dt []中的第二个值，转换1，标准指出它必须小于转换0的dt才能发生...或者换句话说事情发生的必然发生得更快，这与化学逻辑一致。我最关心的是“t + =间隔”线所规定的累积t值可能不完全公平...因为t1点火与t0点火无关，因此t0点火并且比如说.1秒，不应该排除t1使用相同的.1秒来解雇...但我正在努力解决这个问题...欢迎提出建议！这是从脚本中打印出来的详细信息，包括过渡1和2的触发： 时间转换A B C D. 
    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非常合法。     
    2019-05-10 0 0
    分享
    为什么被折叠? 0 个回复被折叠

    要回复问题请先登录或注册

    或代码 OrCode.com 备案号：粤ICP备15020848号-1
    Escape time: 0.035156011581421