ただのメモ

他の人に見せても良い方のメモ

from numba import jit で高速に

numbaのjitというものを使うと普通のpythonの計算がとても速くやってくれるという噂を聞きました。

 

試しに使ってみることにします。

 

始めましょうか。

 

 

===========================

from numba import jit
import time

@jit
def function():
    number = 10 ** 8
    x = [0] * number
    total_time = 10
    time_mesh = 10 ** 3
    dt = total_time / time_mesh
    a = 0.01
    b = 0.02
    for i in range(number - 1):
        x[i + 1] = x[i] + dt * (a * x[i] - b * x[i] ** 2)

t0 = time.time()
function()
t1 = time.time()
print('Calculation time=', float(t1 - t0), '[s]')

 

==========================

計算にかかった時間が、

numba jit 1.5秒

negative control 12.6秒

となった。

 

さらに、計算回数を10倍にして、増やすと、

numba jit 18.1秒

negative control 147.6秒

となった。

 

やはり、多少コードの変更は必要だが、numbaを使った方が圧倒的に速いことがわかった。

 

参考文献はこちら

 

 

生存関数とハザード関数

生存時間解析の話で、生存関数とハザード関数の二つが登場する。

両者の関係についてざっくり理解したい。こちら。

 

生存関数というのは、生存時間を変数にとって、その変数以上、生きられる確率はいくらかを表す。

例えば、仮に1)0歳から100歳までしか生きられないし、2)死亡率が一定の生き物がいたら、

その生存関数は0歳の時は1で、100歳の時0となる、一次関数となる。

 

ハザード関数は、生存時間は変数にとって、その変数時点まで生きたとして、その瞬間に死ぬ確率、確率密度関数を表す。

例えば、さっきの生き物なら、ハザード関数に50歳を代入したら、50歳まで生きて、50歳で死ぬ確率、となる。

50歳まで生きられる確率は、0.5。50歳から50+Δt歳の間に死ぬのは、Δt / 100 

よって、求める値は、 \lim_{\Delta t \rightarrow 0} 0.01 \Delta t / 0.5 \Delta t = 0.02

この生き物だと、年を追うごとに、ハザード関数が上昇していく。

 

気持ち的には、

大域的に、集団全体を描写するときには、生存関数みたいな、何歳以上生きてるのがどれくらい、というのが使えて、

局所的な、今何歳なんだけどどれくらい死にやすいかな、というのが使える、

というイメージ

heteroclinic orbit

studyhacker.net

言語化能力を高めるための手法として、アウトプットする、ということがある。

ここでいう、アウトプットする、とは、自分が聞いたことや調べたことを、自分なりの言葉にして相手に伝える、という一連の動作である。

 

言語化することのメリットは、3つ。

自分の思考の一貫性を保つことができる。

自分の思考の客観視ができる。

自分の思考の要約力が高められる。

 

このアウトプットの際に役に立つものとして、5W1Hというものがある。

これを使いながら、アウトプットの精度を上げていきたい。

 

そういうわけで、今日も調べたことをアウトプットする。

 

ヘテロクリニック軌道

en.wikipedia.org

 

 

ODEの解を軌道を見た時に、過去に極限と飛ばしたら出てくる点と、

未来に極限に飛ばしたら出てくる点がある。

そして、その2点を結ぶ線、それこそが、ヘテロクリニック軌道。

 

記号力学系では、ある周期的な配列から、別の周期的な配列に落ち着くまでの、記号列を指しているようだ。(?)

Brusselator

Brusselator model 

Brusselator - Wikipedia

 

 b = 1 + a^2を超えるか否かでダイナミクスが変化するらしい

hopf 分岐の現象の例を理解する上で良い。

Euler methodでモデルする。

 

======================

import numpy as np
import matplotlib.pyplot as plt
import random
import time
import matplotlib.animation as animation


def main():
    time_start = time.time()
    fig, ax = plt.subplots(figsize=(10,10))

    def f(x,y):
        a = 1;
        b = 2.5
        dt = T / N
        x_diff = a + x ** 2 * y - b * x - x
        y_diff = b * x - x ** 2 * y
        return (dt * x_diff, dt * y_diff)

    ims = []
    x = np.random.random(100)
    y = np.random.random(100)

    T = 10 * 20
    N = 10000 * 5

    for i in range(N):
        dx, dy = f(x, y)
        x, y = x + dx, y + dy
        if i % 100 == 0:
            im1 = ax.plot(x,y,color = "red",linestyle = "None",marker = "o")
            ims.append(im1)

    ani = animation.ArtistAnimation(fig,ims,interval = 10)

    ax.set_aspect("equal")
    time_end = time.time()
    plt.show()
    print(time_end - time_start)


if __name__ == '__main__':
    main()

========================

 

 

 

Sylvester's law of inertia

対称行列の固有値の正、零、負の個数の保存則

シルヴェスターの慣性法則 - Wikipedia

 

確認してみる。

=====================

import numpy as np

A = np.array([[1,2,3],[1,2,-1],[3,-1,3]])

w,v = np.linalg.eig(A)
print(w)

S = np.array([[1,0,1],[0, 4, 1],[-1,-1,4]])
D = S @ A @ S.T
w_d, v_d = np.linalg.eig(D)
print(w_d)

Henon Map

エノン写像

パラメータa,bを変化させると、ダイナミクスが変化する。

それを体験してみることにした。

 

=================================

import numpy as np
import matplotlib.pyplot as plt
import random
import time
import matplotlib.animation as animation


def main():
    time_start = time.time()
    fig, ax = plt.subplots(figsize=(10,10))

    def f(x,y):
        a = 1.4
        b = 0.01
        return (y, a - b * x - y ** 2)

    ims = []
    x = np.random.random(100)
    y = np.random.random(100)

    N = 1000
    for i in range(N):
        x1, y1 = f(x, y)
        im1 = ax.plot(x1,y1,color = "red",linestyle = "None",marker = "o")
        x, y = x1, y1
        ims.append(im1)

    ani = animation.ArtistAnimation(fig,ims,interval = 100)

    ax.set_aspect("equal")
    time_end = time.time()
    plt.show()
    print(time_end - time_start)


if __name__ == '__main__':
    main()

誤り訂正についてのリンク

誤り訂正の考え方

誤り検出符号の考え方

誤り訂正符号の考え方

 

0か1の2値データを扱っている時に、情報通信ないしは、貯蓄の時に、データが書き変わってしまうことがある。

その時に、これはおかしい、と誤りを検出して、誤りを訂正する作業は情報通信において非常に重要である。

 

あらかじめ、情報に余分な仕組みで、誤りが発生したときに起動するようなものを用意しておく必要がある。

 

なので、元の情報に、何かしら余剰を持たせるための、符号化を行う。

ここでは、生成行列を用意して、それをもとに符号を作る。行列をかけたら、実際の情報より、長い列ができる。それを伝達する。

 

そこから、情報伝達後に、パリティ検査行列という、行列を噛ませる。すると、誤りがない場合は0、誤りがある場合は0でない何かしらを返す。

今回は、そこから誤りが発生した位置を取り出している。

 

もっと学ぼうとすると、線形代数整数論が絡んでくるらしい。非常に面白そうである。

numeron

import math
"""
https://y.honkakuha.success-games.net/game/su-115-hitandblow/
"""

N = int(input("what is N?"))
number = [i for i in range(10)]
#number = [i for i in range(26)]
#alphabet = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z']

def nlp(x):
    integer_string =
    for i in range(len(x)):
        integer_string.append(number[alphabet.index(x[i])])

    return integer_string

def inv_nlp(x):
    str_string =
    for i in range(len(x)):
        str_string.append(alphabet[number.index(x[i])])

    return str_string


def making_numeron(N):
    temp = [[i] for i in range(len(number))]
    for i in range(N-1):
        next_temp =
        for j in range(len(number)):
            for k in range(len(temp)):
                if j not in temp[k]:
                    next_temp.append(temp[k] + [j])

        temp = next_temp
    return temp


numeron = making_numeron(N)
#print(len(numeron))

def check(a,b,c,d):
    eat = 0
    bite = 0
    #a equal to candidate number, b equal to guess
    for i in range(N):
        if a[i] in b:
            if a[i] == b[i]:
                eat += 1
            else:
                bite += 1
    if c == eat and d == bite:
        return 1
    else:
        return 0

 

def candidate(temp, b, c, d):
    next =
    for i in temp:
        if check(i, b, c, d):
            next.append(i)

    return next
    #calculate probable candidate

 

def hint(a,b):
    eat = 0
    bite = 0
    for i in range(N):
        if a[i] in b:
            if a[i] == b[i]:
                eat += 1
            else:
                bite += 1
    return (eat, bite)

def predict(temp):
    #calculate least damage solution
    #temp equal to set of possible solutions
    minimum = 10**8
    get_answer = 0
    predicted_answer = [-1] * N
    #answer = [0] * (N ** 2)
    #print(len(temp))

    for i in numeron:
        #if get_answer == 1:
        #    break
        if i in temp:
            priority = 1
        else:
            priority = 0

        answer = [0] * (N ** 2+ 2 * N+1)

        for j in temp:
            hint_eat, hint_bite = hint(i,j)
            #print("ok")
            #if hint_eat == N:
            #    get_answer = 1
            #    true_answer = j
            #    break
            #else:
            answer[N*hint_eat + hint_bite] = answer[N*hint_eat + hint_bite]+ 1
        else:
            if sum(answer) == 1:
                if answer[N*N] == 1:
                    get_answer = 1
                    true_answer = i
        if minimum > max(answer):
            predicted_answer = i
            #print(answer,get_answer,i,)
            minimum = min(minimum, max(answer))
        elif minimum == max(answer) - priority:
            if priority == 1:
                predicted_answer = i

    #print(max(answer))
    if get_answer == 1:
        return true_answer
        pass
    else:
        return predicted_answer

 

def main():
    now = numeron
    print(len(now),"start")
    guess = number[:N]
    while True:
        #q = input("game continue? (y/n):")
        #if q == "y":
        #    pass
        #else:
        #    print("see you again!")
        #    break

        #guess_str = list(str(input("what is your guess?")))
        #guess = nlp(guess_str)
        #guess = [int(element) for element in guess_integers]

        try:
            eat_now = int(input("how many eat?"))
            bite_now = int(input("how many bite?"))
        except:
            print("See you!")
            break
        """
        try:
            eat_now = nlp(input("how many eat?"))
            bite_now = nlp(input("how many bite?"))
        except:
            print("See you!")
            break
        """
        if eat_now == N:
            print("Clear!")
            break
        now = candidate(now, guess, eat_now, bite_now)
        #print(len(now),now)

        prediction = predict(now)
        guess = prediction
        print(f"guess is {prediction}")
        print(f"nubmer of candidate answer :{len(now)}")
        #print(inv_nlp(prediction))

main()

メモ(施設配置の最適化など)

施設配置を最適化する問題がある。

 

https://www.jstage.jst.go.jp/article/bjsiam/23/4/23_KJ00008992858/_pdf

 

いくつか代表的な手法があって、

p-median

p-center

集合カバー

最大カバー

の4つを紹介している。

 

需要(人)とその方のいる場所が組として与えられる。

そして、施設を置く場所の候補も与えられる。

人は必ず一つの施設しか利用できない。

その時、何かしらの尺度を最適化するような、施設の置き方を考える。

 

p-medianでは、利用者と利用施設の距離の総和の最小化を考える。

p-centerでは、利用者と利用施設の距離の最大値の最小化を考える。

集合カバーでは、利用者を利用施設の距離が一定の値以下にできるような最小の施設数を考える。

最大カバーでは、施設の数が固定されていて、それで利用者と利用施設の距離が一定以下という制約のもとで、利用できる人の数(需要)を最大化することを考える。

 

 

メモ(Burgers Equationなど)

Burgers方程式なるものがある。これは非線形偏微分方程式だが、

Cole-hopf変換で、単純な(熱)の拡散方程式に帰着する。

バーガース方程式 - Wikipedia

 

PDE関連の話で、Homotopyが登場したので、立ち寄る

そもそも、Homotopyは、ある位相空間Xから別の位相空間Yへと渡す、「2種類」の連続関数の間の関係を扱い、

位相空間Xと単位区間の直積を、位相空間Yへ移す写像で、

単位区間が0の時、連続関数f、単位区間における1の時、連続関数gの役割を果たすようなものだ。

要は、単位区間を移動するにつれて、連続的に別の関数に切り替わっていくものだろう。

Wikiのドーナツからコップの写像のイメージが直感的に良い。

 

ここで、Nullhomotopicというのは、定数関数とhomotopicであるということだ。

Nullhomotopicに関連して、Lusternik-Schnirelmann categoryというものがある。

 

位相空間XのLusternik-Schnirelmann categoryとは、homotopy不変量で、Xの開被覆で、それぞれの包含写像がNullhomotopicであるとき、そのような開被覆の要素数の最小整数である。

 

Brezis Lieb lemmaというものがある。

これは何かというと、(X,μ)が測度空間で、FnがX上の可測複素関数列でほとんど至る所で関数Fに収束する。

このときLemmaは、ある正の数pについて、|f|^pと、|fn|^p と、|f- fn|^pの関係を等式で繋いだもので、変分問題を解くときに使えるらしい。

 

Klein-Gordon equationなるものがある。これは、波動方程式に一次の項が追加されたもので、相対論的粒子に量子論演算子を噛ませると、導出される式でもある。

クライン-ゴルドン方程式 - Wikipedia

 

https://en.wikipedia.org/wiki/Sine-Gordon_equation

 

 

メモ(脳動脈瘤と神経症状、尿道炎)

前交通動脈と、中大脳動脈の分岐部、内頸動脈ー後交通動脈分岐部、の3つで脳動脈瘤が好発する。

 

内頸ー後交通動脈分岐部動脈瘤では、動眼神経が、

内頸ー眼動脈分岐部では、視神経、視交叉部圧迫症状が起きる

内頸動脈海綿静脈洞部では、Ⅲ,Ⅳ,Ⅴ1、Ⅴ2、Ⅵ脳神経圧迫

前交通動脈瘤では、視交叉圧迫症状が起こる

 

淋菌性尿道炎(淋菌)、潜伏期が短い、排尿痛・外尿道口からの発赤、尿道口からの膿性分泌物、Gram陰性双球菌、淋菌培養、尿道狭窄、セフェム系、スペクチノマイシン

非淋菌性尿道炎(クラミジア最多)、潜伏期長い、排尿痛、漿液性分泌物、クラミジアDNA診断、精巣上体炎、テトラサイクリン系、マクロライド系、ニューキノロン

 

悪性腫瘍か慢性炎症は体重減少を来しうる

大腸癌では血便、慢性膵炎では脂肪の吸収不良による白色便

肝機能障害では血中ビリルビン増加で尿濃染

手掌紅斑は、肝疾患や肝硬変の所見

Osler結節は塞栓

頸静脈虚脱は循環血液量減少

うっ血乳頭は、視神経乳頭の腫れ、頭蓋内圧亢進が原因

項部硬直は、髄膜炎くも膜下出血

 

羽ばたき振戦では、CTかMRIはダメ。ぶれるから

 

硬膜穿刺後頭痛の原因は、硬膜の穿刺孔から出る髄液漏出

それによる、脳圧低下や脳支持組織牽引。

輸液や安静、カフェイン、自己血パッチ。

 

手術後、疼痛、腫脹を伴う突然の下腿浮腫となれば、深部静脈血栓症

チオペンタールは脳活動抑制、脳圧低下。気管支痙攣、舌根沈下、血圧低下、嘔吐。

 

揮発性吸入麻酔薬

エーテル、はろたん、セボフルラン

吸入麻酔薬

亜酸化窒素

静脈麻酔薬

プロポフォールチオペンタールケタミンミダゾラム

オピオイド

モルヒネフェンタニル、レミフェンタニル

筋弛緩薬

サクシニルコリン(脱分極性)ロクロニウム(非脱分極性)ベクロニウム(非脱分極性)、ツボクラリン(非脱分極)

 

局所麻酔薬

メピバカイン、ブピバカイン、ロピバカイン

 

モルヒネヒスタミン遊離、喘息には禁忌

コデインは鎮咳作用のあるオピオイド

セボフルラン、揮発性吸入麻酔薬、気管支拡張作用

チオペンタール全身麻酔の導入に使用。しゃっくり、咽頭痙攣、気管支痙攣。呼吸抑制。気管支喘息に禁忌