Matlab
お役立ちリンク
「メモ代わり」さんの「個人的MATLAB tips」：かゆいところに手がとどくことが書いてある！神経科学者っぽいです。
こまごましたこと
判別用
isnan NaN があるかどうか
isinf Inf or -inf があるかどうか
isinteger 整数値があるかどうか
デバッグ用
error('test') エラーメッセージが出て処理が止まる
warning('test') warning メッセージが出る。処理は止まらない
pause(n) n 秒処理が止まる
fprintf('model %d 終了 \n', model) コマンドにメッセージを出す
構造体
parameter.kf
parameter.kb
という感じで書くと、parameterと1つの変数に二種類の値を入れられる。便利。
並列化のメモ
といっても parfor を用いた並列化。Matlab Distributed Computing Serverは使ったことがありません。それでも数倍シミュレーションが早くなると感動します。
- Hyper-threading対応のCPUなら、4 core/cpuが8 core/cpu になりたくさんコアが使える。BIOSをいじろう。
- 全コアをparforでまわしちゃうと、ほかの仕事ができなくなる・・・
- 並列化計算の意味をよく考えて、並列計算にすべきところをParforにする。Parforで回るときは、完全にバラバラに回る。したがって、そのparforループ内は互いに独立に完結しないといけない。
- 例えば、初期条件を数個ふってみて、それをparfor内で並列に計算させることは可能。なぜならそれらの初期条件の間で情報のやり取りは無いから。
- あるループの情報が次のループに必要な場合とかは並列化無理。例えばオイラー法のように、逐次結果を更新していくような計算の場合は並列計算できない。
- Parforは、以下のように宣言する
- parfor n = 1:1:10 って感じ。もちろん、n が整数であること。
- Parforのなかの変数、データの取り扱いに注意する。
- parfor内の計算結果を出力するためには、sliced variableにする。
- 例えば、
- data = cell(10,1);
- parfor n = 1:10
- y = ode23(...)
- data{n} = y;
- end
- な感じ。 n のところが sliced variable っぽい？
- セル配列にするとsliced variableにしやすい気がする。
if matlabpool('size') ~= 0
matlabpool close;
end
非線形Fitting（式が分かっている場合）
例として、解離定数の式。
- [Kd, r, J, COVB, mse] = nlinfit(xdata, ydata, @(Kd, x) x./(Kd(1)+x), Kd0)
- % Kd パラメーター; r 残差; J ヤコビアン; COVB 共分散行列, mse 誤差項の推定分散
- ci = nlparci(Kd, r, 'Jacobian', J)
- % ci 95%信頼区間（デフォルト値）;
非線形Fitting（式が分かっていない場合。例えば微分方程式の解析解がないときなど。）
局所的な最適化でよければ fminsearch を使う。
このとき、CellDesignerで出力した.mファイルをode functionファイルとして用いる場合は、functionのなかにfunctionを使う入れ子関数にして、パラメーターをふり、自分で最適化関数を作って最小化する。
なかなかめんどくさい・・・
一例
function fitpara = optimization_test()
% 最適化のプログラム例
% 入れ子関数を使って、解析解の無い微分方程式の最適化をする。
% 構成としては、入れ子関数を使っている。
% optimization_test 実行ファイル。functionファイル
% -model.m ODEを含むモデル部分。functionファイル。
% -costfcn.m コスト関数。function。
% 実験例
% In vitro の ERK リン酸化のパラメーターを最適化により求める。
% 逐次反応モデル
% [npERK]' = -k1*[npERK]
% [pYERK]' = k1*[npERK]-k2*[pYERK]
% [pTpYERK]' = k2*[pYERK]
% ちなみに、このモデルは解析解があるのでfminsearchしなくても実は最適化できる・・・
% 実験結果（column 1, npERK; column 2, pYERK; column 3, pTpYERK）
init = ...
[2.20096 0 0
2.117027895 0.067465619 0.016466485
2.058372739 0.142102058 0.000485203
1.542959495 0.627330246 0.030670259
1.050852802 1.045774858 0.10433234
0.490494103 1.344476356 0.365989541
0.17618356 1.30477363 0.72000281
0.055793023 0.839737568 1.305429409];
% 実験結果の反応時間の値 (単位は秒）
tspan = [0;50;100;300;600;1200;1800;3600];
% 初期値パラメーター
x0 = [init(1,1) init(1,2) init(1,3)];
% 最適化処理のための初期値設定
para = [0.01 0.01]; % この初期値近傍から最適化パラメーターを探索する。
% 最適化処理のためのオプション。自分で調べて。
options = optimset('Display','OFF','PlotFcns',@optimplotfval,'TolX',1e-3,'TolFun',1e-3,'MaxIter',250);
% 最適化処理。
[fitpara, norm] = fminsearch(@costfcn, para, options);
% コスト関数。errorの値を最小にするようにサーチする。
function error = costfcn(para)
k1 = para(1);
k2 = para(2);
% ここで ode 関数を使う。
[t, fittedData] = ode23(@model,tspan,x0);
% 入れ子関数。関数の中に関数を入れる。これだと、ode関数の中の変数を最適化できる。
function xdot = model(time,fittedData)
xdot = zeros(3,1);
% 初期値の代入
npERK = fittedData(1);
pYERK = fittedData(2);
pTpYERK = fittedData(3);
%npERK
xdot(1) = -k1*npERK;
%pYERK
xdot(2) = k1*npERK-k2*pYERK;
%pTpYERK
xdot(3) = k2*pYERK;
end
% 目的関数。この関数を最小化するようにパラメーターをフィッティングする。
% つまり最小二乗法。
% sum一回だけでは3行それぞれの和が出るので、
% 二回かけることで、その三行分のまとめた和が出る。
error = sum(sum((fittedData(:,1:3)-init(:,1:3)).^2));
end
% グラフ化
figure(1)
plot(tspan, init(:,1),'bo')
hold on
plot(tspan, init(:,2),'go')
plot(tspan, init(:,3),'ro')
plot(tspan, fittedData(:,1),'b-')
plot(tspan, fittedData(:,2),'g-')
plot(tspan, fittedData(:,3),'r-')
hold off
xlabel('Time (sec)');
ylabel('ERK concentration [/\muM]');
end
グラフ用のoptions （自分用）
% Figureウィンドウの大きさ指定。結構便利。
% gcf は current figure のプロパティを指定するときに使う
set(gcf,'Position',[50,300,800,300]);
%================================================
% plot
%================================================
% グラフ化
h = plot(X_axis, Y_axis);
% 一般
h = plot(x, y, 'linestyle_marker_color','linewidhth',2)
% ラインスタイル ('-', '--', ':'); など
% マーカタイプ ('x', '*', 'o','+'); など
% カラー識別子('c', 'm', 'y', 'k', 'r', 'g', 'b', 'w');
h = plot(x, y, 'go--'); % 緑の破線、丸マーカつき
% 線幅、色
set(h,'linewidth',2,'color','blue');
% X, Y 軸
xlabel('Time(min)');
ylabel('kcat-cRaf [/\muM/sec]');
xlim([0 1]);
ylim([0 2]);
axis([0 12 -0.5 1]);
% gca は current axis のプロパティを指定するときに使う
set(gca,'XScale','Log'); % X 軸片ログ
set(gca,'YScale','Log'); % Y 軸片ログ
% タイトル
title(['ERK-' char(C(c)) ']);
% レジェンド
legend(h,'First','Second','Third');
%================================================
% hist
%================================================
histではなくBar関数を使うと色や重なりなどを変更しやすくなる。
%================================================
% pcolor (Heat mapを作るときに使う。)
%================================================
% pcolor 表示
% X_axis (１行X列ベクトル), Y_axis (Y行１列ベクトル), C(X,Y) 行列。
% C の次元をちゃんとX,Yに合うようにすること。
h = pcolor(X_axis, Y_axis, C);
% ラインスタイルをなしにする。ほぼ必須。
set(h,'linestyle','none');
% グラフの横にカラーバーつける。必須
colorbar;
% カラーバーの上限、下限を決める
caxis([0, 0.74]);
% 新しいMatlabでは色が変なのでjetにする
colormap('jet')
% 通常は左下が原点になるが、axis ijをいれると左上が原点になる。
axis ij;
% pcolorで行列データを格子状に表現する。こちらを参考にさせていただきました。
A = [1,0,0;0,1,0;0,0.5,0.7]; % 視覚化対象の行列
[x,y] = size(A); % 行列の大きさを測る
pcolor([A,zeros(x,1);zeros(1,y+1)]); % 行列サイズを１行１列を大きくしてプロット
Saveするときのデータに日付を入れる
データをsaveするときに日付を入れてoverwriteを防ぐ
save(['allData' datestr(now, 'yyyymmdd_HHMMSS')], '-v7.3');
で、allData20121223_210517.matというファイル名で保存される。
2012年12月23日21時5分17秒って感じ。
データ型
Matlabとか以前の話ですが・・・
- Zeros: デフォルトで double 。
- int8: -128 ~ 127。 めちゃ軽くなる。
CellDesignerから出力された.mファイルを使ってシミュレーションする
- 下の方にあるいらない行を削る
- function z = pow ....以降削る。
- ファイル名と関数名（mファイルの一行目）を一致させる。
- ode関数でまわす
- 関数名のみタイプすると初期値が得られる（引数nargin）
- [y, t] = ode23(@function_name, [0 100], x0)
- % x0は初期値
- ode23やode45で遅い場合は、微分方程式がstiffである可能性が高い（陽解法）。ode15s （陰解法）とかにすると爆速になったりする。計算誤差が大きくなるみたいだけど、今のところ誤差はほとんど感じない。