Drkcore

27 06 2009 Python DMPK Tweet

pythonでdeconvolution

演習本にデコンボリューション法が載っていたのだけど、実装までは載ってなかったので、調べたら、論文が見つかった。

ProductName ファーマコキネティクス―演習による理解
杉山 雄一,山下 伸二,加藤 基浩
南山堂 / ¥ 6,300 ()
在庫あり。

GOTOとIFの組み合わせがwhileなのかforなのか悩むとこが幾つかあったが、pythonで書きなおした。結構、fortranっぽさが残っている気もするが。

#!/usr/bin/env python

def calc_eps():
    eps = 1.0
    while eps + 1.0 > 1.0: eps = eps/2.0
    return eps

def sign(a1,a2):
    if a2 >= 0: return abs(a1)
    else:       return -abs(a1)

def deconv(a,div,b,dose):
    n = len(a); m = len(b); l = n + m -1
    eps = calc_eps()
    u = [0.0 for i in range(l)]
    v = [0.0 for i in range(l)]
    gamma = [0.0 for i in range(20)]
    g = [0.0 for i in range(20)]

    a.sort(lambda x,y:cmp(x[1],y[1]))

    for II in range(n-1):
        x1 = a[II][1] * -1
        x2 = a[II+1][1] * -1
        qx1 = 0.0
        qx2 = 0.0

        for i in range(n):
            s1 = 1.0
            s2 = 1.0
            for j in range(n):
                if j != i:
                    s1 = s1*(a[j][1] + x1)
                    s2 = s2*(a[j][1] + x2)
            qx1 = qx1 + a[i][0] * s1
            qx2 = qx2 + a[i][0] * s2

        x3 = x1
        qx3 = qx1
        d = x2 -x1
        e = d

        while(1):

            if abs(qx3) < abs(qx2):
                x1  = x2; x2  = x3; x3  = x1
                qx1 = qx2; qx2 = qx3; qx3 = qx1

            tol = 4.0*eps*abs(x2)
            xm  = 0.5*(x3-x2)

            if abs(xm) <= tol or qx2 == 0.0: break

            if abs(e) >= tol and abs(qx1) > abs(qx2):
                if x1 == x3:
                    s = qx2/qx1; p = 2.0*xm*s; q = 1.0-s
                else:
                    q = qx1/qx3; r = qx2/qx3; s = qx2/qx1
                    p = s*(2.0*xm*q*(q-r) - (x2-x1)*(r-1.0))
                    q = (q-1.0)*(r-1.0)*(s-1.0)

                if p > 0.0: q = -q
                p = abs(p)

                if 2.0*p < 3.0*xm*q - abs(tol*q) and p < abs(0.5*e*q):
                    e = d
                    d = p/q
                else:
                    d = xm
                    e = d
            else:
                d = xm
                e = d

            x1 = x2
            qx1 =qx2
            x2 = x2 + d
            if abs(d) <= tol: x2 = x2 - d + sign(tol,xm)
            qx2 = 0.0
            for i in range(n):
                s2 = 1.0
                for j in range(n):
                    if j != i: s2 = s2*(a[j][1] +x2)
                qx2 = qx2 + a[i][0] * s2
            if qx2*(qx3/abs(qx3)) > 0.0:
                x3 = x1
                qx3 = qx1
                d = x2 - x1
                e = d

        gamma[II] = x2
        s2 = 0.0
        for j in range(n):
            s1 = 0.0
            for k in range(n):
                if k != j:
                    s1 = s1 + 1.0/(gamma[II] + a[k][1])
            s2 = s2 + a[j][0]*s1/(gamma[II] + a[j][1])
        g[II] = 1.0/s2

    ak1 = 0.0
    ak2 = 0.0
    ak3 = 0.0
    for i in range(n):
        ak1 = ak1 + a[i][0];
        ak2 = ak2 + a[i][0]*a[i][1]
        if n > 1 and i != n-1: ak3 = ak3 + g[i]/gamma[i]
    ak1 = 1.0/ak1
    ak2 = ak1*ak1*ak2
    ak3 = ak2 - ak3
    ak4 = 100.0*div/dose
    uz = 0.0

    for i in range(l):
        s = 0.0
        if i > m-1:
            v[i] = -gamma[i-m]
            for j in range(m): 
                s = s + b[j][0]/(gamma[i-m] + b[j][1])
            u[i] = ak4*g[i-m]*s/gamma[i-m]
        else:
            v[i] = b[i][1]
            u[i] = ak4*b[i][0]*(ak1 - ak3/b[i][1])
            if n != 1:
                for j in range(n-1):
                    s = s + g[j]/(gamma[j]*(gamma[j] + b[i][1]))
                u[i] = u[i] - ak4*b[i][0]*s
        uz = uz - u[i]

    return uz,u,v

if __name__ =='__main__':
    print "=== Data set 1 ==="
    a     = [[1.0427,5.0526],[0.97617,0.97567]]
    div   = 1.0
    b     = [[-0.80547,3.6373],[0.80547,0.83033]]
    dose  = 0.6

    print deconv(a,div,b,dose)

    print "=== Data set 2 ==="
    a     = [[1.3377,3.1851],[0.52381,0.61065]]
    div   = 1.0
    b     = [[-2.2637,2.0800],[2.2637,1.2513]]
    dose  = 0.6

    print deconv(a,div,b,dose)

    print "=== Data set 3 ==="
    a     = [[1.0308,5.9131],[1.0487,1.0262]]
    div   = 1.0
    b     = [[-1.6103,3.3563],[1.6103,1.3917]]
    dose  = 0.6

    print deconv(a,div,b,dose)

    print "=== Data set 4 ==="
    a     = [[1.3377,3.1851],[0.52381,0.61065]]
    div   = 1.0
    b     = [[-5.5440,2.2488],[5.5440,1.7521]]
    dose  = 0.6

    print deconv(a,div,b,dose)

Table.1 のテスト計算を実行。論文のデータセット4のa2はタイポでしょう。

=== Data set 1 ===
(103.38116460095375, 
[99.756988719460651, -23.218031509886707, -179.92012181052769], 
[3.6373000000000002, 0.83033000000000001, 2.9469592648362699])
=== Data set 2 ===
(94.012343334076377, 
[212.39909284817011, 2395.320119577972, -2701.7315557602187], 
[2.0800000000000001, 1.2513000000000001, 1.3350740721242431])
=== Data set 3 ===
(94.362883301636657, 
[-1704.7606134762168, 73.013373899587776, 1537.3843562749923], 
[3.3563000000000001, 1.3916999999999999, 3.4906828227939406])
=== Data set 4 ===
(91.15932635164495, 
[370.51731307178676, -1111.1884882257646, 649.51184880233291], 
[2.2488000000000001, 1.7521, 1.3350740721242431])

About

  • もう5年目(wishlistありマス♡)
  • 最近はPythonとDeepLearning
  • 日本酒自粛中
  • ドラムンベースからミニマルまで
  • ポケモンGOゆるめ

Tag

Python Deep Learning javascript chemoinformatics Emacs sake and more...

Ad

© kzfm 2003-2021