MMA Advent Calendar 2016/12/6
そろそろブログまとめようと思ったので移転記事です.
移転元 https://wiki.mma.club.uec.ac.jp/su_zu/MMA%20AdventCalendar%202016
MMA Advent Calendar 2016/12/6
これはMMA Advent Calendar 2016の記事です
PowerShellでぐいぐいっとGUI動かす方法を紹介する。Widnows使っている人は気軽にPowerShellで実行してみて頂ければと思います。
先に言っておくが最近MSが頑張っていてLinuxでもPowerShellが動くがこれから行うものがバリバリ.NETに依存しているのでWindowsを用意してほしい。ちなみにWPFはWineでも対応していない。Windowsを宗教的理由で使えない方はまあできるんだ分かる分かるみたいな気持ちで見ていただければ幸いです。
サクッと.NET使いたいときはPowerShellが便利だと私は思っている。それならGUIも意外といけるのではと思いPowerShellでGUIを動かしてみた 最初に結果を言ってしまうがまあできるが高機能なGUIはやはりしんどいし書いてて面倒くさくなったので当初Twitterクライアント作ろうかなと思ったがTwitter検索ツールに変更した。
WPF
WPFとはXAMLというXMLベースで書くことができるユーザーインターフェイス部分とロジック部分を分離してGUIを作成する仕組みです。WPF以前にはWindows Formsなどがあったが時代はXMLベースのGUI作成に流れが傾いている気がする(?)
まずXAMLを書く。そのまま書いてUIを組むのはしんどいのでVisual StudioとかBlendとかを作って頑張って欲しい。そしてそのXAMLを読み込んでwindowを作る。すると以下のようなものを作ることができる。
Add-Type -AssemblyName presentationframework [xml]$xaml = @' <Window xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation" xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:d="http://schemas.microsoft.com/expression/blend/2008" xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006" Title="mmaAdventCalendar" Height="350" Width="600"> <StackPanel> <TextBlock FontSize="30" Text="炎炎炎!!!✌('ω'✌ )三✌('ω')✌三( ✌'ω')✌"></TextBlock> <Image Source="http://d250g2.com/d250g2.jpg"></Image> </StackPanel> </Window> '@ $reader=(New-Object System.Xml.XmlNodeReader $xaml) $window=[Windows.Markup.XamlReader]::Load( $reader ) $window.ShowDialog() | Out-Null
Twitter検索ツール
途中で飽きて雑なコードになっているのでまあこんなことできるよぐらいに思ってみてくれれば幸いです。
なぜか分からないが以下のスクリプトを一つのファイルにしてPowerShellで実行するとエラーがでる。最初にAdd-Type -AssemblyName presentationframeworkでDLL読み込んでから動かすと普通に動くので誰か教えて欲しい。
Add-Type -AssemblyName presentationframework class Twi { static [xml]$xaml = @' <Window xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation" xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:d="http://schemas.microsoft.com/expression/blend/2008" xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006" Title="mmaAdventCalendar" Height="350" Width="600"> <Grid> <Grid.RowDefinitions> <RowDefinition Height="50"></RowDefinition> <RowDefinition></RowDefinition> </Grid.RowDefinitions> <Grid Grid.Row="0"> <Grid.ColumnDefinitions> <ColumnDefinition></ColumnDefinition> <ColumnDefinition Width="50"></ColumnDefinition> </Grid.ColumnDefinitions> <TextBox Grid.Column="0" FontSize="30" Name="keyword"></TextBox> <Button Grid.Column="1" Name="search">search</Button> </Grid> <ScrollViewer Grid.Row="1"> <Grid> <Grid.RowDefinitions> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> <RowDefinition Height="100"></RowDefinition> </Grid.RowDefinitions> <Grid.ColumnDefinitions> <ColumnDefinition Width="50"></ColumnDefinition> <ColumnDefinition></ColumnDefinition> </Grid.ColumnDefinitions> <Canvas Grid.Row="0" Grid.Column="0"> <Image Width="50" Height="50" Name="img0"></Image> </Canvas> <Canvas Grid.Row="1" Grid.Column="0"> <Image Width="50" Height="50" Name="img1"></Image> </Canvas> <Canvas Grid.Row="2" Grid.Column="0"> <Image Width="50" Height="50" Name="img2"></Image> </Canvas> <Canvas Grid.Row="3" Grid.Column="0"> <Image Width="50" Height="50" Name="img3"></Image> </Canvas> <Canvas Grid.Row="4" Grid.Column="0"> <Image Width="50" Height="50" Name="img4"></Image> </Canvas> <Canvas Grid.Row="5" Grid.Column="0"> <Image Width="50" Height="50" Name="img5"></Image> </Canvas> <Canvas Grid.Row="6" Grid.Column="0"> <Image Width="50" Height="50" Name="img6"></Image> </Canvas> <Canvas Grid.Row="7" Grid.Column="0"> <Image Width="50" Height="50" Name="img7"></Image> </Canvas> <Canvas Grid.Row="8" Grid.Column="0"> <Image Width="50" Height="50" Name="img8"></Image> </Canvas> <Canvas Grid.Row="9" Grid.Column="0"> <Image Width="50" Height="50" Name="img9"></Image> </Canvas> <Canvas Grid.Row="10" Grid.Column="0"> <Image Width="50" Height="50" Name="img10"></Image> </Canvas> <Canvas Grid.Row="11" Grid.Column="0"> <Image Width="50" Height="50" Name="img11"></Image> </Canvas> <Canvas Grid.Row="12" Grid.Column="0"> <Image Width="50" Height="50" Name="img12"></Image> </Canvas> <Canvas Grid.Row="13" Grid.Column="0"> <Image Width="50" Height="50" Name="img13"></Image> </Canvas> <Canvas Grid.Row="14" Grid.Column="0"> <Image Width="50" Height="50" Name="img14"></Image> </Canvas> <Canvas Grid.Row="15" Grid.Column="0"> <Image Width="50" Height="50" Name="img15"></Image> </Canvas> <Canvas Grid.Row="16" Grid.Column="0"> <Image Width="50" Height="50" Name="img16"></Image> </Canvas> <Canvas Grid.Row="17" Grid.Column="0"> <Image Width="50" Height="50" Name="img17"></Image> </Canvas> <Canvas Grid.Row="18" Grid.Column="0"> <Image Width="50" Height="50" Name="img18"></Image> </Canvas> <Canvas Grid.Row="19" Grid.Column="0"> <Image Width="50" Height="50" Name="img19"></Image> </Canvas> <TextBox Grid.Row="0" Grid.Column="1" Name="text0" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="1" Grid.Column="1" Name="text1" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="2" Grid.Column="1" Name="text2" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="3" Grid.Column="1" Name="text3" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="4" Grid.Column="1" Name="text4" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="5" Grid.Column="1" Name="text5" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="6" Grid.Column="1" Name="text6" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="7" Grid.Column="1" Name="text7" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="8" Grid.Column="1" Name="text8" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="9" Grid.Column="1" Name="text9" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="10" Grid.Column="1" Name="text10" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="11" Grid.Column="1" Name="text11" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="12" Grid.Column="1" Name="text12" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="13" Grid.Column="1" Name="text13" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="14" Grid.Column="1" Name="text14" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="15" Grid.Column="1" Name="text15" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="16" Grid.Column="1" Name="text16" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="17" Grid.Column="1" Name="text17" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="18" Grid.Column="1" Name="text18" IsReadOnly="True" TextWrapping="Wrap"></TextBox> <TextBox Grid.Row="19" Grid.Column="1" Name="text19" IsReadOnly="True" TextWrapping="Wrap"></TextBox> </Grid> </ScrollViewer> </Grid> </Window> '@ [System.Windows.Controls.ContentControl]$window [string]$consumer_key = "OlCmnCziFLdEe308ZvEJckX6U" [string]$consumer_secret = "knoLolvTxQ1Had8jCcPbvWyD4lYt6f6anw27uAXmWc5GapNkId" [Hashtable]$header [System.Windows.Controls.Primitives.TextBoxBase]$keyword [System.Windows.Controls.Primitives.ButtonBase]$search [System.Windows.FrameworkElement[]]$imgs [System.Windows.FrameworkElement[]]$texts [System.Management.Automation.PSMethodInfo]$searched Twi(){ $this.init() } Twi($consumer_key, $consumer_secret){ $this.consumer_key = $consumer_key $this.consumer_secret = $consumer_secret $this.init() } init(){ $this.imgs = New-Object System.Windows.FrameworkElement[] 20 $this.texts = New-Object System.Windows.FrameworkElement[] 20 $bearer = [Twi]::get_bearer($this.consumer_key, $this.consumer_secret) $this.header = @{ "Authorization" = "Bearer " + $bearer; } $reader=(New-Object System.Xml.XmlNodeReader ([Twi]::xaml)) $this.window=[Windows.Markup.XamlReader]::Load( $reader ) for ($i=0; $i -lt 20; $i++){ $this.imgs[$i] = $this.window.FindName("img" + $i.ToString()) $this.texts[$i] = $this.window.FindName("text" + $i.ToString()) } $this.keyword = $this.window.FindName("keyword") $this.search = $this.window.FindName("search") $this.keyword.Text = "UEC lang:ja" $this.searched = $this.search.add_Click $t = $this $this.searched.Invoke({ $t.update() }) $this.update() $this.window.ShowDialog()| Out-Null } [Object[]]getStatuses(){ [string]$q = $this.keyword.Text $q = $q.Replace(" ", "+") $uri = "https://api.twitter.com/1.1/search/tweets.json?count=20&q=" + $q $response = Invoke-RestMethod -Method Get -Headers $this.header -Uri $uri return $response.statuses } update(){ $statuses = $this.getStatuses() foreach ($status in $statuses){ $index=$statuses.IndexOf($status) $this.imgs[$index].Source = $status.user.profile_image_url_https $this.texts[$index].Text = $status.text } } static [string]base64($str){ $bytes = [System.Text.Encoding]::UTF8.GetBytes($str) return [System.Convert]::ToBase64String($bytes) } static [string]get_bearer($consumer_key, $consumer_secret){ $basic = [Twi]::base64($consumer_key + ":" + $consumer_secret) $auth_header = @{ "Authorization" = "Basic " + $basic; "Content-Type" = "application/x-www-form-urlencoded;charset=UTF-8" } $auth_body = "grant_type=client_credentials" $response = Invoke-RestMethod -Method Post -Headers $auth_header -Body $auth_body -Uri "https://api.twitter.com/oauth2/token" return $response.access_token; } } [Twi]::new()
MMA入って思うこと
みんなMMAに入って思うこと書いてる気がするので僕も書く。
最初にMMAの部室に行ったときは誰か覚えていないがSSHとVPNの説明してもらって「苦しんで覚えるc言語」はいい本だよって教えてもらってすぐその時に生協に行って買って部室で読みながらコード書いていた気がする。そのあとAOJとか説いていってナップサック問題を自力で解こうとするが動的計画法を学んでこんなの自力で解けるわけなかったんだと時間を無駄にした感を味わってこれはしんどいなと思って一時期AOJやめてGUIとかLinuxの環境構築楽しんでいた気がする。当初MMAでのIDで呼び合う文化に抵抗がありあんまり気が進まなかったが環境に適応していっているのかわからないが今ではIDは分かるけど本名わからないが当たり前になってる。あとMMAに入りはじめのころは部員同士の「わかる!!」とか永遠手を手を叩ハイコンテクストな会話にもついていけなかったが今では分かるし分かりすぎるので人間なんとかなるもんだなぁと思う。
明日はbenevolent0505さんが記事を書いてくれるみたいですね。
STDP学習ベースのSNNでMNISTを分類する
雑に論文みて勉強したので何も保証しません.SNNの情報全然無いような気がするのでしんどい.
Spiking Neural Network(SNN)
ニューロンの活動電位(スパイク)に注目したニューラルネットワーク.現状,計算量が多いのでチップを作ったりなどのハードウェアの研究が盛んな傾向にあるような気がする.IBMのTruneNorthとかIntelのLoihiなどはSNNにあたる.
Spike-Timing-Dependent Plasticity(STDP)
STDPはスパイクタイミングに依存してシナプスの重みが変わる仕組みである. ざっくり言うと細胞Aが発火して,細胞Bが発火するとA->B間のシナプス結合が強まるLTP(長期増強)が起きて,逆に細胞Bが発火して,細胞Aが発火するとA->B間のシナプス結合が弱まるLTD(長期抑圧)が起きる.
シナプス前細胞の発火時間を,シナプス後細胞の発火時間をとする.を以下のように定義する.
シナプスの結合の重みの更新は以下の式で表される.STDPの式は色々提案されているが今回は基本的なものを使用する
ニューロンモデル
ネットワーク
入力はMNISTのそれぞれのピクセルデータを最大60Hzのポアソンスパイクモデルにして入力させる. よってneural codingはrate codingを使用する. 入力を受け取る興奮性の細胞が発火した場合,他の細胞を抑制させる(winner-take-all).
学習後,最終的により多くの反応(発火回数)を示すニューロンをその教師ラベルと紐づける. テストでは,入力に対して最も多くの反応を示したニューロンを抽出し,そのニューロンと紐付けられているラベルと比較する.
結果
以下が学習させた900個のニューロンの重みを表した図である.
accuracy: 0.775
ソースコード
参考
- Frontiers | Unsupervised learning of digit recognition using spike-timing-dependent plasticity | Frontiers in Computational Neuroscience
- Simplified spiking neural network architecture and STDP learning algorithm applied to image classification
- github.com
- Spike-timing dependent plasticity - Scholarpedia
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
- 脳の計算論