Implementation of Poisson Spike model using Rust
自分用のメモ
Poisson Spike model(ポアソンスパイクモデル)
ある時間内に発火する頻度が,ポアソン分布に従うスパイクモデルである. 実際のニューロンは不応期などがあり,連続で発火することはないが,このモデルはその辺を考慮していない. ただ,大脳新皮質のスパイクはポアソンスパイクモデルと似たようなスパイクが観測されるらしい.
ポアソン分布は以下のような式で表される.
は平均回数を表す.
以下は,ポアソンスパイクモデルの実装
fn generate_spike(fr: f64, dt: f64, t: f64, rng: &mut XorShiftRng) -> Vec<usize> { // fr: 発火頻度 Hz // dt: 微小量 // t: 生成する時系列の長さ sec let num = (t/dt) as usize; let v: Vec<usize> = rng.gen_iter().take(num).map(|x: f64| if x < fr*dt { 1 } else { 0 } ).collect(); // 1が発火を表す v }
例えばの発火頻度のポアソンスパイクは1秒間に平均で60回発火が起きる.
以下でのポアソン分布と1秒間に発火頻度のポアソンスパイクが発火した回数を比較した実験を行った.
extern crate gnuplot; extern crate rand; use rand::{Rng, XorShiftRng, SeedableRng}; use rand::distributions::{Poisson, Distribution}; use std::mem::transmute; pub fn generate_spike(fr: f64, dt: f64, t: f64, rng: &mut XorShiftRng) -> Vec<usize> { let num = (t/dt) as usize; let v: Vec<usize> = rng.gen_iter().take(num).map(|x: f64| if x < fr*dt { 1 } else { 0 } ).collect(); v } fn main() { let seed = 1; let x: u32 = 123456789; let y: u32 = 362436069; let z: u32 = 521288629; let w: u32 = seed; let seeds: [u32; 4] = [x, y, z, w]; const N: usize = 1000000; const LAMBDA: f64 = 60.0; let seeds = unsafe {transmute(seeds)}; let mut rng = XorShiftRng::from_seed(seeds); let mut fg = gnuplot::Figure::new(); { // poisson spike let mut spike_cnts = Vec::with_capacity(N); for _ in 0..N { let spikes = generate_spike(LAMBDA, 0.0001, 1.0, &mut rng); let spike_cnt = spikes.iter().sum::<usize>(); spike_cnts.push(spike_cnt); } let max_cnt = spike_cnts.iter().max().unwrap(); let xs = (0..*max_cnt).collect::<Vec<usize>>(); let mut spike_hist = vec![0;*max_cnt+1]; for spike_cnt in &spike_cnts { spike_hist[*spike_cnt] += 1; } let ax = fg.axes2d(); ax.lines(xs.iter(), spike_hist.iter(), &[gnuplot::Caption("poisson spike"), gnuplot::Color("blue")]); // poisson distributions let mut spike_cnts = Vec::with_capacity(N); for _ in 0..N { let poi = Poisson::new(LAMBDA); let v = poi.sample(&mut rng); spike_cnts.push(v); } let max_cnt = spike_cnts.iter().max().unwrap(); let mut spike_hist = vec![0;(*max_cnt+1) as usize]; for spike_cnt in &spike_cnts { spike_hist[*spike_cnt as usize] += 1; } ax.lines(xs.iter(), spike_hist.iter(), &[gnuplot::Caption("poisson distributions"), gnuplot::Color("red")]); } fg.echo_to_file("ps.plt"); }
良さそう.
参考
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をブロードキャストアドレスに飛ばしてなぜ返ってこないのか分からない