0

*2010年6月17日編集

私は自分のコードを改善する方法を理解しようとしています(よりPythonicにします)。また、生化学で一般的なシナリオを説明する、より直感的な「条件文」を書くことに興味があります。以下のプログラムの条件付き基準は、回答#2で説明しましたが、コードに満足していません。正常に機能しますが、明確ではなく、より複雑な条件付きシナリオに実装するのは簡単ではありません。アイデアは大歓迎です。コメント/批評を歓迎します。最初の投稿体験@stackoverflow-必要に応じてエチケットについてコメントしてください。

このコードは、次の演習の解決策となる値のリストを生成します。

「選択したプログラミング言語で、GillespieのFirst Reaction Algorithmを実装して、反応A ---> Bの時間的動作を研究します。この場合、AからBへの遷移は、別の化合物Cが存在する場合にのみ発生します。ここで、以下のペトリネットでモデル化されているように、Cは動的にDと相互変換します。100分子のA、1のCがあり、反応の開始時にBまたはDが存在しないと仮定します。kABを0.1s-1に設定します。 kCDとkDCの両方を1.0秒に-1。100秒以上のシステムの動作をシミュレートします。」

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の単純な確率シミュレーションの背後にある数学に関する情報:

通常、このようなプロセスは離散イベントとしてシミュレートされ、各イベントは特定の速度定数「k」と時間間隔「dt」で発生する可能性のあるイベント「n」の数が与えられた確率「P」で発生します。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の数値はイベントが発生したことを示し、.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

ギレスピーアルゴリズムはわかりませんが、プログラムが正しい平衡状態に収束していることを確認したと思います。したがって、私はあなたの質問を次のように解釈します

「これが実用的な物理プログラムです。どうすればもっとPythonicにすることができますか?」

次のようなことをする方がおそらくもっとPython的でしょう

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

Noneを返す代わりにStopIteration例外を発生させる方が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はおそらくスパース配列である必要がありますが、これらのアイデアが出発点になることを願っています。

于 2010-06-17T10:30:56.640 に答える
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 []の3つの値のため、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 []の2番目の値である遷移1の場合、基準は、遷移0のdtよりも小さくなければならないことを示しています...つまり、まったく起こったのはもっと早く起こったに違いない、それは化学論理と一致している。私の最大の懸念は、「t + =間隔」の線で義務付けられている累積t値が完全に公平ではない可能性があることです...t1の発火はt0の発火とは独立しているため、t0の発火とたとえば.1秒は同じの使用からt1を除外します。発砲まで1秒...しかし、私はこれの修正に取り組んでいます...提案を歓迎します!これは、トランジション1と2の起動を含む、スクリプトからの詳細な出力です。

時間遷移ABCD

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

0.0 0100 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 に答える