IBM Qやってみた

こちらの動画をまとめてみる

  • 量子コンピュータとは、量子性を積極的に利用して計算する計算機。
  • 量子性とは重ね合わせの原理や、量子もつれ。
  • 古典コンピュータとは、データのもたせ方やルールが違う。
  • 使いたい理由としては、1）計算量の大きいものを早く解くため、2）量子化学、量子機械学習、だが、特定の問題に対してのみ適応。
  • 出現しうる全データの重ね合わせで表現される。
  • 測定とは、量子ビットが0か1かを読み出す処理。重ね合わせの原理により、確率的に0か1が出現する状態を作れる。
  • 量子計算では、最後に測定のタイミングまでに、求める結果の出現確率が1に近づくように実装する。（実際、これがGroverのアルゴリズムに対応）
  • 量子状態はベクトル、ゲートは行列になる。
  • ベクトルは足し引きが同じ成分同士では出来る。定数倍、掛け算も出来る（列を行に一個ずつかけていく）掛け算は非可換！
  • 1つのビットは2つの値。
  • 量子ゲートでは、X（0と1入れ替える）、Y、Z（1の状態の確率振幅符号を逆にする）、H（0か1の状態から重ね合わせ状態を作る）、CNOT（上の状態が1なら下の状態にXゲートをする。それ以外は何もしない）の他に、S、Tがある。何もしないのも、ある意味ゲート。
  • IBM Q、Qiskit、IBM Q Experienceがある。
  • IMB Qは実際に使えるので、やってみる。
  • 最初に量子回路を作る。

    f:id:medical-science:20211010005928p:plain

  • この理論的な状態の確率分布はこのようになる。

    f:id:medical-science:20211010010227p:plain

  • 実際にIBMの量子コンピュータに投げると、結果が帰ってくる（待ち時間はそこそこ）

    f:id:medical-science:20211010005757p:plain

  • ノイズがある（431，53，47，493）けれども、まあ大体は予測出来そう。
  • もうちょっと複雑な回路でやってみる価値があるが、あまり期待はしない。

ゲート方式は諦めて、アニーリング方式を採用した方が良さそう。

Sparse matrix with diffusion equation

疎行列を使って一気に解く

 

def RDE():
    dtnum = 4801 #時間をいくつに分割して計算するか 一の位を1とした方が良い 時刻tのfor文の終わりの数字に関係あり
    dxnum = 101 #xとyをいくつに分割して計算するか
    dynum = 101

    thick = 10 #x方向の大きさ
    width = 10 #y方向の大きさ
    time_calc = 48 #計算する時間
    #beta = 0.1 #上の式中ではラムダ　温度拡散率 lambda = 熱伝導率/(密度*比熱)
    
    #constant for PDE
    D = 0.001
    
    """計算準備"""
    #空の解を用意
    solution_k1 = np.zeros([dxnum,dynum])
    solution_k2 = np.zeros([dxnum,dynum]) + 100 * numpy.random.random*1 #solution_k2の初期値が初期条件になる
    
    #境界条件
    #今回はNeumann boundary conditionを作りたい
    
    
    dt = time_calc / (dtnum - 1)
    dx = thick / (dxnum - 1)
    dy = width / (dynum - 1)
    r_x = (dt / (2 * dx ** 2))
    r_y = (dt / (2 * dy ** 2))
    
    
    """Ax=bのAを用意"""
    a_matrix_1 = np.identity(dxnum) * (1 + 2 * D * r_x) \
        + np.eye(dxnum, k=1) * (- D * r_x) \
            + np.eye(dxnum, k=-1) * (- D * r_x)
    a_matrix_1[0,0] /= 2
    a_matrix_1[-1,-1] /= 2
    
    a_matrix_2 = np.identity(dynum) * (1 + 2 * D * r_y) \
        + np.eye(dynum, k=1) * (- D * r_y) \
            + np.eye(dynum, k=-1) * (- D * r_y)
    a_matrix_2[0,0] /= 2
    a_matrix_2[-1,-1] /= 2
    
    #疎行列を格納　csr方式
    a_matrix_csr1 = csr_matrix(a_matrix_1)
    a_matrix_csr2 = csr_matrix(a_matrix_2)
    
    a_matrix_arr1 = (a_matrix_csr1 for i in range(dynum))
    a_matrix_arr2 = (a_matrix_csr2 for i in range(dxnum))
    
    #さらに大きな行列を作る
    a_matrix_csr_agg1 = csr_matrix(block_diag(a_matrix_arr1))
    a_matrix_csr_agg2 = csr_matrix(block_diag(a_matrix_arr2))
    
    #ADI法ではk+1時刻とk+2時刻を計算するので、for文の回数は半分になる
    #number = Decimal(str(dtnum/2)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)
    number = int(dtnum/2)
    
    solution = np.zeros([dxnum,dynum,number]) #解を代入する行列を作成
    
    #ADI法の計算
    
    for k in range(number): #時刻tについて
        #大きな行列で一気に計算を済ませてしまう。
        #for loopが減るので、計算時間が１/3になった。
        #一回目の拡散
        template_1 = np.zeros*2
        template_1[1 : dxnum + 1, 1 : dynum + 1] = solution_k2
        template_1[:, 0] = template_1[:, 2]
        template_1[:, -1] = template_1[:, -3]
        
        b_array_x = solution_k2 \
            + D * r_y * (template_1[1: dxnum + 1, 0: dynum] \
                - 2 * solution_k2 \
                    + template_1[1: dxnum + 1, 2: dynum + 2]) \
                + dt * (1 / 2) * function(solution_k2)
        b_array_x[0, :] /= 2
        b_array_x[-1,:] /= 2
        
        #一気に解くために、1次元配列にする。
        #ただし、足し合わせる方向性を変えるために、一度転置行列にしている
        b_array_x_1d = b_array_x.T.reshape(1, dxnum * dynum)[0]
        temp_solve1 = spla.spsolve(a_matrix_csr_agg1, b_array_x_1d)
        #得られた解を2次元配列に戻す。
        solution_k1 = temp_solve1.reshape(dynum,dxnum).T
        
        
        template_2 = np.zeros*3
        template_2[1:dxnum+1, 1:dynum+1] = solution_k1
        template_2[0,:] = template_2[2,:]
        template_2[-1,:] = template_2[-3,:]
        
        b_array_y = solution_k1 \
            + D * r_x * (template_2[0: dxnum, 1: dynum + 1] \
                - 2 * solution_k1 \
                    + template_2[2: dxnum + 2, 1: dynum + 1]) \
                + dt * (1 / 2) * function(solution_k1)
        b_array_y[:, 0] /= 2
        b_array_y[:,-1] /= 2
        
        #一気に解くために1次元配列にする。
        
        b_array_y_1d = b_array_y.reshape(1 , dxnum * dynum)[0]
        temp_solve2 = spla.spsolve(a_matrix_csr_agg2, b_array_y_1d)
        #得られた解を2次元配列に戻す。
        solution_k2 = temp_solve2.reshape(dxnum,dynum)
        
        solution[:,:,k] = solution_k2
        
    return solution

 

 

 

 

 

*1:dxnum, dynum

*2:dxnum + 2, dynum + 2

*3:dxnum + 2, dynum + 2

block_diagonal_matrix

import numpy as np 
from scipy.sparse import csr_matrix, block_diag

def main():
    a = np.array([[1,1,0],[1,2,0],[3,0,4]])
    z = np.zeros*1
    b = np.array([[1,11,0],[1,21,0],[33,0,0]])
    mat_array = (a for i in range(10))
    c_sp = block_diag(mat_array).toarray()
    print(c_sp)
    
main()

 

 

scipy.sparceのblock_diagを使えば、行列をブロックとして扱い、それを対角方向に並べた疎行列を作ることが出来る。これが、大規模行列計算に役に立つとか。

 

詳しいdocumentはこちらを参照。

*1:3,3

Matrix reshaping (2d to 1d, 1d to 2d)

import numpy as np 


def convert_1d_array(arr, x, y):
    return arr.reshape(1,x * y)

def convert_2d_array(arr, x, y):
    return arr.reshape(x,y)

def main():
    a = np.ones*1
    a_1d = convert_1d_array(a, 5, 3)
    a_2d = convert_2d_array(a_1d, 5, 3)
    print(a_1d)
    print(a_2d)
    
    
main()

 

大規模の行列計算をするときに、こういった2次元配列を1次元配列に落とし込んだ方が計算が早く済むことがある。

そのメモ書きである。

*1:5,3

Bernoulli distribution with multiple probability

import numpy as np 

p = [0.1, 0.2, 0.9, 0.1]#異なる確率
n = 1

#10000回ベルヌーイ分布でサンプルを引く
x = np.array([np.random.binomial(n, p) for i in range(10000)])
y = np.array([np.sum(x[:,i]) for i in range(4)])/10000
print(p, "original")
print(y, "experiment")#確率が元の確率と似ているか

 

ベルヌーイ分布のpythonの記事を見てみると、異なるベルヌーイ分布からサンプリングするものについてのまとまった記事が見当たらなかったので、備忘録として残しておく。

 

そもそもベルヌーイ分布とは、当たりの確率の決まったガチャガチャ見たいなもの。

０か１かを返す。

詳しいこと（分散、期待値、エントロピー、特性関数など）はこちらのWikiを参照して欲しい。

 

 

乳歯について

乳歯は上下２０本

切歯から生え始める。

乳歯は、生後６−８ヶ月から生え始める。２−３歳で生え終わる。

6歳頃から永久歯に生え替わり始める。12歳に生え終わる。

２８本と、親知らず４本である。

 

基本的に、乳歯は月齢−６が本数の数え方になる。

 

そのほかの、以前の小児の記事はこちら。

 

medical-science.hatenablog.com

medical-science.hatenablog.com

medical-science.hatenablog.com

 

頭蓋の門について、大泉門、小泉門

門とは、wikiで調べてみる。

建造物の世界では、”敷地と外部を区切るものに、（通るために）開けられた出入り口のことらしい

 

頭蓋における門とは、

頭蓋骨が頭蓋の外と中を区切っている。

そこに開けられた穴のことを言っている感じだろう。

 

大泉門、小泉門は、その門のお名前のことだろう。

 

何が泉なのか、それは気になる。

大、小は、穴の”大きさ”を持って名付けたのか。

 

穴の閉鎖時期が、それぞれ、生後１歳半、生後１−２ヶ月だから、穴の閉じるスピードが同じだとすると、穴の面積の比が割り出せそうだ。（実際はそんな簡単じゃないだろう）

 

もう少し周りの構造物との、位置関係を学ぶ。

 

ここで言う位置関係とは、3次元空間上の座標ではなく、隣接している構造物（領域）の名前を知ろう、と言うことだ。

 

イメージ的にはネットワーク、グラフ辺りな気持ち。

 

話を戻して、大泉門は、左右の前頭骨と、左右の頭頂骨の４つの骨の間の穴。

 

小泉門は、左右の頭頂骨、後頭骨の３つの間の穴。

 

小児発育身長あれこれ

medical-science.hatenablog.com

 

 

Kaup指数：小児の発育の指標で、BMIに１０をかけたもの

15 - 18にあって欲しい

 

出生時の平均は５０cm。4頭身。

乳児は、出生後１年間で25cm/year

４歳で出生時の２倍、つまり１００cmになる。

幼児では、7cm/yearで増える。

つまり、25cm増えるタイミングが、1歳と4歳となっている

 

思春期は、7-9cm/yearとなっている。

女性の方が２年思春期が早い

 

12歳で１５０cmくらい。