% 隠れ状態とマルコフランダム場演習用資料
% 科学技術振興事業団川人学習動態脳プロジェクト
% 岡田真人
% プログラム1用サンプル

clear
close all

length = 30	% データの長さ
                % 式(6)参照
stepUp = 11	% step-upポイント   bar_d_iが0から1に変化する点
                % 式(8)(9)参照
stepDown = 21	% step-downポイント bar_d_iが1から0に変化する点
                % 式(8)(9)参照
lower_val = 0.0 % i=1から10までとi=21から30までのbar_d_iの値
                % 式(8)参照
upper_val = 1.0	% i=11から20までのbar_d_iの値
                % 式(9)参照
sigma_d = 0.05	% d_iに関するガウスノイズの標準偏差
                % 式(10)参照
lambda = 6.0	% lambda
                % 式(1)参照
epsilon_f   = 0.0001
                % epsilon_f
                % 式(5)参照
loop = 100	% 繰り返し回数

% 配列 bar_d[] を作成する
% 式(8)(9)参照

% まず全てのbar_d_iにlower_val=0を代入する
bar_d = lower_val * ones(length,1);

% 次にi=11から20までのbar_d_iにupper_val=1を代入する
bar_d(stepUp:stepDown-1) = upper_val * ones(stepDown-stepUp,1);

% ガウス分布に従う乱数のタネをcputimeからひろう
randn('seed', cputime)

% d_iを設定する．式(10)参照
d = bar_d + ( sigma_d *randn( size(bar_d) ) );

% f_iの初期状態はd_iに設定する
f = d;

% 式(7)の周期的境界条件のために，i=0とi=31に対するバッファを用意する．
% 具体的にはd_iとf_iの各々に対して配列 dtmp[] ftmp[] を作成する。
% よって配列 dtmp[] ftmp[] は32個の要素を持つ．
% d_iとdtmp_(i+1) が対応する．ftmp[]も同様．
% 計算用の配列 dtmp[] ftmp[] を作成する。
% dtmpの最初の要素にd(30)を代入し，dtmpの最後の要素にd(1)を代入する．
dtmp = [d(length); d; d(1)];

% f_iの初期状態はd_iに設定するので，ftmpの初期状態を以下のように決める．
ftmp = dtmp;

% 実際のdtmpやftmpを使った計算範囲は2から31までになる．
% keisanは計算範囲を指定する．
keisan = 2:length+1;		

% まず，例えばloop=100なら100回オイラー法を繰り返す．
loop_count = loop;

% 繰り返し計算を続けるかどうかのフラグの設定
% flg=1で続ける それ以外でやめる
% まずは計算するから取りあえず flg=1
flg = 1;

% 図の指定
fig1 = figure;

% 繰り返し計算を続けるかどうかのフラグ flg=1で続ける それ以外でやめる
while flg == 1
  
% 式(5)を計算する．
  for cnt=1: loop
% 参考: keisanは計算範囲を指定する．keisan = 2:length+1;
% 以下の例だと，ftmp(2)がdf(1)に対応する．
    df = lambda*((ftmp(keisan-1)-ftmp(keisan)) ...
	         +(ftmp(keisan+1)-ftmp(keisan))) ...
	 - ( ftmp(keisan) - dtmp(keisan) );
    df = df*epsilon_f;
    % f[], ftmp[] を更新する。
    f = f + df;
    ftmp = [f(length); f; f(1)];
  end

  hold off
% d(fの初期値)を青で表示する
  plot(d, 'b');
  hold on
% fを赤で表示する
  plot(f, 'r');
% 横軸は画素の番号iである．
  xlabel('i')

% オイラー法を抜け出すか否かの処理
  str = sprintf('Exit? (loop=%d) y/n [y] : ', loop_count);
% yを入力すると抜け出す nを入力すると続ける処理
  key = input(str,'s');
  if isempty(key);
    key = 'y'	% 入力省略時の処理
  end
  
% yを入力すると抜け出す処理
  if key(1) == 'y'	% 終了
% flgを0に設定してwhile文から抜け出す
    flg = 0;

% nを入力すると続ける処理
  elseif key(1) == 'n'	% 続行
% 繰り返し回数(loop_count)をloop分だ増やす
    loop_count = loop_count + loop;
  end
  
end;