Implementation of FitzHugh-Nagumo model using Rust
自分用のまとめ
FitzHugh-Nagumo(フィッツヒュー・南雲モデル)
HH型モデルを簡略化したモデルである.
それぞれのパラメータは以下を参考にして次のように設定した.
http://www.chem.scphys.kyoto-u.ac.jp/nonnonWWW/b8/07f/majima/majima.pdf
初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード
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"); }
参考
Implementation of Integrate-and-fire model using Rust
自分用のまとめ
Integrate-and-fire(積分発火モデル)
入力の時間積分によって細胞の膜電位が増加し,ある閾値を超えるとスパイクが生成されるという現象を再現したモデルである. HH型モデルのような生物的なニューロンモデルとは異なった挙動をするので注意する必要がある.
それぞれの定数は以下のように設定した.
初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード
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"); }
Implementation of Izhikevich model using Rust
自分用のまとめ
Izhikevich(イジケヴィッチモデル)
Hodgkin–Huxley modelなどのコンダクタンス依存型モデルは複雑な式で計算量もそれなりにいるので,簡易な形に落とし込んだモデルがIzhikevich modelである. 以下のような2変数の常微分方程式で表すことができる.
ここで変数とは,Izニューロンが発火するたびに以下のような規則で状態を更新する.
の値が30になった場合,の値をcにし,にを加算する.
if v = 30 mV then v <- c, u <- u+d
それぞれの定数は以下のように設定した.
初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード
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"); }
参考
- Which Model to Use for Cortical Spiking Neurons?
- Fundamentals of Computational Neuroscience
- 脳の計算論
Implementation of Hodgkin–Huxley model using Rust
自分用のまとめ
Hodgkin–Huxley model(ホジキン・ハクスレーモデル)
神経細胞の細胞膜をコンデンサ,イオンチャネルを変化する抵抗素子として考えた電気回路モデル(回路図は割愛)である. ナトリウム,カリウムを通すイオンチャネルが存在し,それぞれに平衡電位を持つ.
キルヒホッフの法則の法則から以下の式が成り立つ.
,はナトリウムチャンネル,カリウムチャンネルが通り流れ出す電流で,は漏れ電流である.
それぞれの電流は以下の式で表される.
変数は以下のxを読み替えた同じ式に従う.
それぞれの定数は以下のように設定した.
初期値問題にはルンゲ=クッタ法を使用した. 以下実装したコード
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"); }
参考
- 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()
すごい無理矢理感があるので,誰かもっと良い方法を知っていたら教えてほしい.
ARP tableを監視して部屋にいるか簡易通知する方法
行きたい先に目的の人がいるかどうか簡易的に知るにはどうすべきか考えたときに同じサブネットにラップトップやスマートフォンを接続しているなら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; }
参考