suzuzusu日記

(´・ω・`)

Implementation of FitzHugh-Nagumo model using Rust

自分用のまとめ

FitzHugh-Nagumo(フィッツヒュー・南雲モデル)

HH型モデルを簡略化したモデルである.

{\displaystyle
\frac{dV(t)}{dt} = c(-\frac{V(t)^{3}}{3} + V(t) - w + I(t)) \\
\frac{dw}{dt} = V(t) - bw +a \\
}

それぞれのパラメータは以下を参考にして次のように設定した.

http://www.chem.scphys.kyoto-u.ac.jp/nonnonWWW/b8/07f/majima/majima.pdf

{ \displaystyle
a = 0.7 \\
b = 0.8 \\
c = 10 \\
}

初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード

extern crate gnuplot;

fn rk4<F: Fn(f64)->f64>(f: F, y: f64, dt: f64) -> f64 {
    let k1 = dt * f(y);
    let k2 = dt * f(y + k1*0.5);
    let k3 = dt * f(y + k2*0.5);
    let k4 = dt * f(y + k3);
    (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
}

fn main(){
    const A: f64 = 0.7;
    const B: f64 = 0.8;
    const C: f64 = 10.0;

    let mut t = -30.0;
    let dt = 0.001;
    let mut v = 0.0;
    let mut w = 0.0;
    let mut i_ext ;

    let mut x = Vec::new();
    let mut y = Vec::new();

    while t <= 100.0 {

        if t >= 30.0 && t <= 70.0 {
            i_ext = 0.35;
        } else {
            i_ext = 0.0;
        }

        let d_w = move |y:f64| v - B * y + A;
        w += rk4(d_w, w, dt);

        let d_v = move |y:f64| C*(-y.powf(3.0)/3.0 + y - w + i_ext);
        v += rk4(d_v, v, dt);

        if t >= 0.0 {
            y.push(v);
            x.push(t);
        }

        t += dt;
    }

    let mut fg = gnuplot::Figure::new();
    fg.axes2d().lines(x.iter(), y.iter(), &[gnuplot::Color("blue")]);
    fg.echo_to_file("fn.plt");
}

f:id:suzuzusu:20180327092810p:plain
実行結果

github.com

参考

Implementation of Integrate-and-fire model using Rust

自分用のまとめ

Integrate-and-fire(積分発火モデル)

入力の時間積分によって細胞の膜電位が増加し,ある閾値を超えるとスパイクが生成されるという現象を再現したモデルである. HH型モデルのような生物的なニューロンモデルとは異なった挙動をするので注意する必要がある.

{\displaystyle
\tau \frac{dV(t)}{dt} = E_{L} - V(t) + RI(t) \\
if\  V(t) \geq V_{\theta} \  then \  V(t) = V_{0}
}

それぞれの定数は以下のように設定した.

{ \displaystyle
E_{L} = -65 \\
\tau = 10 \\
V_{\theta} = -55 \\
RI = 12
}

初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード

extern crate gnuplot;

fn rk4<F: Fn(f64)->f64>(f: F, y: f64, dt: f64) -> f64 {
    let k1 = dt * f(y);
    let k2 = dt * f(y + k1*0.5);
    let k3 = dt * f(y + k2*0.5);
    let k4 = dt * f(y + k3);
    (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
}

fn main(){
    const E_L: f64 = -65.0;
    const THETA: f64 = -55.0;
    const TAU: f64 = 10.0;
    const RI: f64 = 12.0;

    let mut t = 0.0;
    let dt = 0.001;
    let mut v = E_L;

    let mut x = Vec::new();
    let mut y = Vec::new();

    while t <= 100.0 {

        if v >= THETA {
            v = E_L;
        }

        let d_v = move |y:f64| (E_L - y + RI)/TAU;
        v += rk4(d_v, v, dt);

        y.push(v);
        x.push(t);

        t += dt;
    }

    let mut fg = gnuplot::Figure::new();
    fg.axes2d().lines(x.iter(), y.iter(), &[gnuplot::Color("blue")]);
    fg.echo_to_file("if.plt");
}

f:id:suzuzusu:20180325134816p:plain
実行結果

github.com

Implementation of Izhikevich model using Rust

自分用のまとめ

Izhikevich(イジケヴィッチモデル)

Hodgkin–Huxley modelなどのコンダクタンス依存型モデルは複雑な式で計算量もそれなりにいるので,簡易な形に落とし込んだモデルがIzhikevich modelである. 以下のような2変数の常微分方程式で表すことができる.

{ \displaystyle
\frac{dv}{dt} = 0.04v^{2} + 5v + 140 - bu + I \\
\frac{du}{dt} = a(v-u)
}

ここで変数{v(t)}{u(t)}は,Izニューロンが発火するたびに以下のような規則で状態を更新する.

{v(t)}の値が30になった場合,{v(t)}の値をcにし,{u(t)}{d}を加算する.

if v = 30 mV
then v <- c, u <- u+d

それぞれの定数は以下のように設定した.

{ \displaystyle
a = 0.02 \\
b = 0.2 \\
c = -65 \\
d = 40
}

初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード

extern crate gnuplot;

fn rk4<F: Fn(f64)->f64>(f: F, y: f64, dt: f64) -> f64 {
    let k1 = dt * f(y);
    let k2 = dt * f(y + k1*0.5);
    let k3 = dt * f(y + k2*0.5);
    let k4 = dt * f(y + k3);
    (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
}

fn main(){
    const A: f64 = 0.02;
    const B: f64 = 0.2;
    const C: f64 = -65.0;
    const D: f64 = 6.0;

    let mut t = -0.0;
    let dt = 0.001;
    let mut i_ext ;
    let mut v = -70.0;
    let mut u = B*v;

    let mut x = Vec::new();
    let mut y = Vec::new();

    while t <= 100.0 {

        if t > 30.0 && t < 70.0 {
            i_ext = 14.0;
        } else {
            i_ext = 0.0;
        }

        let d_v = move |y:f64| 0.04*y*y+5.0*y+140.0-B*u+i_ext;
        v += rk4(d_v, v, dt);

        let d_u = move |y:f64| A*(v - y);
        u += rk4(d_u, u, dt);

        if v > 30.0 {
            y.push(30.0);
            v = C;
            u += D;
        }else {
            y.push(v);
        }
        x.push(t);

        t += dt;
    }

    let mut fg = gnuplot::Figure::new();
    fg.axes2d().lines(x.iter(), y.iter(), &[gnuplot::Color("blue")]);
    fg.echo_to_file("iz.plt");
}

f:id:suzuzusu:20180325123127p:plain
実行結果

github.com

参考

Implementation of Hodgkin–Huxley model using Rust

自分用のまとめ

Hodgkin–Huxley model(ホジキン・ハクスレーモデル)

神経細胞の細胞膜をコンデンサイオンチャネルを変化する抵抗素子として考えた電気回路モデル(回路図は割愛)である. ナトリウム,カリウムを通すイオンチャネルが存在し,それぞれに平衡電位を持つ.

キルヒホッフの法則の法則から以下の式が成り立つ.

{ \displaystyle
C \frac{dV(t)}{dt} = -I_{Na}(t) - I_{K}(t) - I_{L}(t) + I(t)
}

{I_{Na}}{I_{K}}はナトリウムチャンネル,カリウムチャンネルが通り流れ出す電流で,{I_{L}}は漏れ電流である.

それぞれの電流は以下の式で表される.

{ \displaystyle
I_{Na}(t) = g_{Na}m(t)^{3} h(t)(V(t) - E_{Na}) \\
I_{K}(t) = g_{K}n(t)^{4}(V(t) - E_{K}) \\
I_{L}(t) = g_{L}(V(t) - E_{L})
}

変数{h, m, n}は以下のxを読み替えた同じ式に従う.

{ \displaystyle
\alpha_{x}(V)(1-x(t)) - \beta_{x}(V)x(t)
}

{ \displaystyle
\alpha_{h}(V) = 0.07 \exp(\frac{-V}{20}) \\
\alpha_{m}(V) = \frac{0.1(25-V)}{\exp\left\{(25-V)/10\right\}-1} \\
\alpha_{n}(V) = \frac{0.01(10-V)}{\exp\left\{(10-V)/10\right\}-1} \\
\beta_{h}(V) = \frac{1}{\exp\left\{(30-V)/10\right\}+1} \\
\beta_{m}(V) = 4 \exp(\frac{-V}{18}) \\
\beta_{n}(V) = 0.125 \exp(\frac{-V}{80}) \\
}

それぞれの定数は以下のように設定した.

{ \displaystyle
E_{Na} = 115 mV \\
E_{K} = -12 mV \\
E_{L} = 10.613 mV \\
g_{Na} = 120 mV \\
g_{K} = 36 mV \\
g_{L} = 0.3 mV \\
C_{m} = 1.0 \mu F / cm^{2} \\
}

初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード

extern crate gnuplot;

fn alpha_m(v: f64) -> f64{
    (0.1*(25.0-v)) / (((25.0-v)/10.0).exp()-1.0)
}

fn alpha_h(v: f64) -> f64{
    0.07*((-v/20.0).exp())
}

fn alpha_n(v: f64) -> f64{
    (0.01*(10.0-v)) / (((10.0-v)/10.0).exp()-1.0)
}

fn beta_m(v: f64) -> f64{
    4.0*((-v/18.0).exp())
}

fn beta_h(v: f64) -> f64{
    1.0/(((30.0-v)/10.0).exp()+1.0)
}

fn beta_n(v: f64) -> f64{
    0.125*((-v/80.0).exp())
}

fn rk4<F: Fn(f64)->f64>(f: F, y: f64, dt: f64) -> f64 {
    let k1 = dt * f(y);
    let k2 = dt * f(y + k1*0.5);
    let k3 = dt * f(y + k2*0.5);
    let k4 = dt * f(y + k3);
    (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0
}

fn main(){
    const E_NA: f64 = 115.0;
    const E_K:  f64 = -12.0;
    const E_L:  f64 = 10.613;
    const G_NA: f64 = 120.0;
    const G_K:  f64 = 36.0;
    const G_L:  f64 = 0.3;
    const C_M:  f64 = 1.0;

    let mut t = -100.0;
    let dt = 0.001;
    let mut i_ext ;
    let mut v = -10.0;
    let mut x_n = 0.0;
    let mut x_m = 0.0;
    let mut x_h = 0.0;

    let mut x = Vec::new();
    let mut y = Vec::new();

    while t <= 100.0 {

        if t > 30.0 && t < 70.0 {
            i_ext = 10.0;
        } else {
            i_ext = 0.0;
        }

        let a_n = alpha_n(v);
        let a_m = alpha_m(v);
        let a_h = alpha_h(v);

        let b_n = beta_n(v);
        let b_m = beta_m(v);
        let b_h = beta_h(v);

        let d_n = |y: f64| a_n*(1.0-y) - b_n*y;
        let d_m = |y: f64| a_m*(1.0-y) - b_m*y;
        let d_h = |y: f64| a_h*(1.0-y) - b_h*y;

        x_n += rk4(d_n, x_n, dt);
        x_m += rk4(d_m, x_m, dt);
        x_h += rk4(d_h, x_h, dt);

        let i_na = G_NA *(x_m.powf(3.0))*x_h*(v- E_NA);
        let i_k = G_K *(x_n.powf(4.0))*(v- E_K);
        let i_l = G_L *(v- E_L);
        let i_sum = i_na + i_k + i_l;

        let d_v = |y:f64| (i_ext - i_sum) / C_M;
        v += rk4(d_v, v, dt);

        if t >= 0.0{
            x.push(t);
            y.push(v);
        }
        t += dt;
    }

    let mut fg = gnuplot::Figure::new();
    fg.axes2d().lines(x.iter(), y.iter(), &[gnuplot::Color("blue")]);
    fg.echo_to_file("hh.plt");
}

f:id:suzuzusu:20180325110018p:plain
実行結果

github.com

参考

  • Fundamentals of Computational Neuroscience
  • 脳の計算論

KerasのstatefulなLSTMのtrain出力を記録する

忘備録として書く.

FAQ - Keras Documentation からstatefulなRNN系だと学習や予測時に状態を更新してしまうので,単純に同じデータを入力して出力を得ることはできない.

そこでメトリック経由で出力を得る.

# 教師データ
def true(y_true, y_pred):
    return y_true[0]

# 出力データ
def pred(y_true, y_pred):
    return y_pred[0]

上記のようにカスタムメトリックから無理矢理出力を得る.あとはmetricsとして追加すればいい.

model.compile(loss='mean_squared_error', optimizer='rmsprop', metrics=[true, pred])

以下,正弦波の予測のExample

import numpy as np
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers import Dense, LSTM

def true(y_true, y_pred):
    return y_true[0]

def pred(y_true, y_pred):
    return y_pred[0]

T = 30
f = 1/T
N = 10000
hidden_unit = 30

model = Sequential()
model.add(LSTM(hidden_unit, stateful=True, batch_input_shape=(1, 1, 1), return_sequences=True))
model.add(LSTM(hidden_unit, stateful=True, return_sequences=True))
model.add(LSTM(hidden_unit, stateful=True))
model.add(Dense(1))
model.compile(loss='mean_squared_error', optimizer='rmsprop', metrics=[true, pred])

# logger
true_log = np.zeros(N)
pred_log = np.zeros(N)

# train
for i in range(N):
    x = np.sin(2*np.pi*i*f)
    t = np.sin(2*np.pi*(i+1)*f)
    history = model.train_on_batch(np.array(x).reshape(1, 1, 1), np.array(t).reshape(1,1))
    true = history[1]
    pred = history[2]
    true_log[i] = true
    pred_log[i] = pred

# plot
plt.plot(true_log[N-3*T:], label='true')
plt.plot(pred_log[N-3*T:], label='prediction')
plt.legend()
plt.show()

f:id:suzuzusu:20180317071711p:plain

すごい無理矢理感があるので,誰かもっと良い方法を知っていたら教えてほしい.

ARP tableを監視して部屋にいるか簡易通知する方法

行きたい先に目的の人がいるかどうか簡易的に知るにはどうすべきか考えたときに同じサブネットにラップトップやスマートフォンを接続しているならARP tableのMACアドレスを見て簡易的に判断できると思ったのでやってみる

  1. 事前にMACアドレスを教えてもらう
  2. pingをサブネット内の全IPに送ってARP tableを更新する
  3. ARP tableを見てMACアドレスがあれば部屋にいると判断する

以下そのスクリプト(PowerShell)

$mac_address_target="xx-xx-xx-xx-xx-xx"

while($TRUE){
    foreach ($i in @(2..254)){
        Start-Process -WindowStyle hidden -FilePath ping -ArgumentList (" -n 1 xxx.xxx.xxx." + $i.ToString())
    }
    sleep 5

    $arp_data = -split ( arp -a | Select-String $mac_address_target )

    if($arp_data.Length -gt 1){
        echo "join"
        # 通知(ex slack, etc...)
        break
    }
    
    sleep 300

}

正直pingをブロードキャストアドレスに飛ばしてなぜ返ってこないのか分からない

Mandelbrot, Intel AVX, CImg

Intel AVXのサンプルを動かしたのでそのメモ

画像出力にはCImgを使った

  • 環境
    • Visual C++ 2017
#include <immintrin.h>
#include <complex>
#include "CImg.h"

using namespace std;
using namespace cimg_library;

// simple code to compute Mandelbrot in C++
void MandelbrotCPU(float x1, float y1, float x2, float y2,
    int width, int height, int maxIters, unsigned short * image)
{
    float dx = (x2 - x1) / width, dy = (y2 - y1) / height;
    for (int j = 0; j < height; ++j)
        for (int i = 0; i < width; ++i)
        {
            complex<float> c(x1 + dx*i, y1 + dy*j), z(0, 0);
            int count = -1;
            while ((++count < maxIters) && (norm(z) < 4.0))
                z = z*z + c;
            *image++ = count;
        }
}

void MandelbrotAVX(float x1, float y1, float x2, float y2,
    int width, int height, int maxIters, unsigned short * image)
{
    float dx = (x2 - x1) / width;
    float dy = (y2 - y1) / height;
    // round up width to next multiple of 8
    int roundedWidth = (width + 7) & ~7UL;

    float constants[] = { dx, dy, x1, y1, 1.0f, 4.0f };
    __m256 ymm0 = _mm256_broadcast_ss(constants);   // all dx
    __m256 ymm1 = _mm256_broadcast_ss(constants + 1); // all dy
    __m256 ymm2 = _mm256_broadcast_ss(constants + 2); // all x1
    __m256 ymm3 = _mm256_broadcast_ss(constants + 3); // all y1
    __m256 ymm4 = _mm256_broadcast_ss(constants + 4); // all 1's (iter increments)
    __m256 ymm5 = _mm256_broadcast_ss(constants + 5); // all 4's (comparisons)

    float incr[8] = { 0.0f,1.0f,2.0f,3.0f,4.0f,5.0f,6.0f,7.0f }; // used to reset the i position when j increases
    __m256 ymm6 = _mm256_xor_ps(ymm0, ymm0); // zero out j counter (ymm0 is just a dummy)

    for (int j = 0; j < height; j += 1)
    {
        __m256 ymm7 = _mm256_load_ps(incr);  // i counter set to 0,1,2,..,7
        for (int i = 0; i < roundedWidth; i += 8)
        {
            __m256 ymm8 = _mm256_mul_ps(ymm7, ymm0);  // x0 = (i+k)*dx 
            ymm8 = _mm256_add_ps(ymm8, ymm2);         // x0 = x1+(i+k)*dx
            __m256 ymm9 = _mm256_mul_ps(ymm6, ymm1);  // y0 = j*dy
            ymm9 = _mm256_add_ps(ymm9, ymm3);         // y0 = y1+j*dy
            __m256 ymm10 = _mm256_xor_ps(ymm0, ymm0);  // zero out iteration counter
            __m256 ymm11 = ymm10, ymm12 = ymm10;        // set initial xi=0, yi=0

            unsigned int test = 0;
            int iter = 0;
            do
            {
                __m256 ymm13 = _mm256_mul_ps(ymm11, ymm11); // xi*xi
                __m256 ymm14 = _mm256_mul_ps(ymm12, ymm12); // yi*yi
                __m256 ymm15 = _mm256_add_ps(ymm13, ymm14); // xi*xi+yi*yi

                                                            // xi*xi+yi*yi < 4 in each slot
                ymm15 = _mm256_cmp_ps(ymm15, ymm5, _CMP_LT_OQ);
                // now ymm15 has all 1s in the non overflowed locations
                test = _mm256_movemask_ps(ymm15) & 255;      // lower 8 bits are comparisons
                ymm15 = _mm256_and_ps(ymm15, ymm4);
                // get 1.0f or 0.0f in each field as counters
                // counters for each pixel iteration
                ymm10 = _mm256_add_ps(ymm10, ymm15);

                ymm15 = _mm256_mul_ps(ymm11, ymm12);        // xi*yi 
                ymm11 = _mm256_sub_ps(ymm13, ymm14);        // xi*xi-yi*yi
                ymm11 = _mm256_add_ps(ymm11, ymm8);         // xi <- xi*xi-yi*yi+x0 done!
                ymm12 = _mm256_add_ps(ymm15, ymm15);        // 2*xi*yi
                ymm12 = _mm256_add_ps(ymm12, ymm9);         // yi <- 2*xi*yi+y0 

                ++iter;
            } while ((test != 0) && (iter < maxIters));

            // convert iterations to output values
            __m256i ymm10i = _mm256_cvtps_epi32(ymm10);

            // write only where needed
            int top = (i + 7) < width ? 8 : width & 7;
            for (int k = 0; k < top; ++k)
                image[i + k + j*width] = ymm10i.m256i_i16[2 * k];

            // next i position - increment each slot by 8
            ymm7 = _mm256_add_ps(ymm7, ymm5);
            ymm7 = _mm256_add_ps(ymm7, ymm5);
        }
        ymm6 = _mm256_add_ps(ymm6, ymm4); // increment j counter
    }
}


int main()
{
    const int width = 512;
    const int height = 512;
    const int num_pixels = width*height;

    unsigned short image[num_pixels];

    //MandelbrotCPU(0.29768, 0.48364, 0.29778, 0.48354, width, height, 4096, image);
    MandelbrotAVX(0.29768, 0.48364, 0.29778, 0.48354, width, height, 4096, image);

    CImg<unsigned short> ci(width, height, 1, 3);
    unsigned short *p_ci = ci.data();
    unsigned short *p_image = image;
    for (int k = 0; k < 3; k++) {
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width; j++) {
                *p_ci = *p_image % (30 * (k + 1));
                p_ci++;
                p_image++;
            }
        }
        p_image = image;
    }
    ci.display();

    return 0;
}

参考

f:id:suzuzusu:20170719032606j:plain

software.intel.com

The CImg Library - C++ Template Image Processing Toolkit