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

clear
close all

length = 30	% データの長さ
                % 式(6)(26)参照
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)(11)(14)参照
amp_phi = 0.2   % 位相の初期値の広がり 一様分布 +- amp_phi*pi
                % 式(28)参照
average = 0.0   % 位相の初期値の平均 0 に設定しても一般性を失わない．
J = 1.0	        % 位相の滑らか拘束
                % 式(11)のJ^Rに対応
epsilon_f   = 0.0001
                % epsilon_f
                % 式(23)参照
epsilon_phi = 0.1
                % epsilon_phi
                % 式(25)参照
loop = 1000	% 繰り返し回数

% 配列 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) ) );

% phi_iの初期値を格納するinit_phi[]を作成する．
% 式(28)参照
% まず全ての要素を average = 0.0 に設定する．
init_phi = average * ones(length,1);

% 一様乱数のタネをcputimeからひろう
rand('seed', cputime)

% [-amp_phi*pi, amp_phi*pi]の一様乱数を加えinit_phi[]を作成する．
init_phi = init_phi + ( amp_phi * 2.0 * pi * (rand(size(init_phi)) - 0.5 ));

% 位相を格納するphi[]の初期状態はinit_phi[]に設定する
phi      = init_phi;

% 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;

% 式(27)の周期的境界条件のために，i=0とi=31に対するバッファを用意する．
% 具体的にはphi_i対して配列 phi_tmp[] を作成する。
% よって配列 phi_tmp[] は32個の要素を持つ．
% phi_iとphi_tmp_(i+1) が対応する．
% 計算用の配列 phi_tmp[] を作成する。
% phi_tmpの最初の要素にphi(30)を代入し，phi_tmpの最後の要素にphi(1)を代入する．
phi_tmp = [phi(length); phi; phi(1)];

% 実際の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
  
% 式(23)(24)(25)を計算する．
  for cnt=1: loop
    % dphi[] を更新する。
    dphi = J * (sin(phi_tmp(keisan-1) - phi_tmp(keisan)) + ...
	sin(phi_tmp(keisan+1) - phi_tmp(keisan)))/2.0 ...
	- lambda * ((ftmp(keisan-1)-ftmp(keisan)).^2.* ...
	sin(phi_tmp(keisan-1) - phi_tmp(keisan)) + ...
	(ftmp(keisan+1)-ftmp(keisan)).^2.* ...
	sin(phi_tmp(keisan+1) - phi_tmp(keisan))) / 4.0 ;
    dphi = dphi*epsilon_phi;
    % df[] を更新する。
    df = lambda*(   (1.0+cos(phi_tmp(keisan-1)-phi_tmp(keisan))).* ...
	           (ftmp(keisan-1)-ftmp(keisan)) ...
	         + (1.0+cos(phi_tmp(keisan+1)-phi_tmp(keisan))).* ...
		   (ftmp(keisan+1)-ftmp(keisan))) / 2.0  ...
	 - ( ftmp(keisan) - dtmp(keisan) );
    df = df*epsilon_f;
    % phi[], phi_tmp[] を更新する。
    phi = phi + dphi;
    % f[], ftmp[] を更新する。
    f = f + df;
    ftmp = [f(length); f; f(1)];
    phi_tmp = [phi(length); phi; phi(1)];
  end

% 位相をpiで規格化する
  normal_init_phi = init_phi / pi;
  normal_phi      = phi / pi;

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

% phiの表示
  hold off 
  fig1 = figure(2);
% init_phi(phiの初期値)を青で表示する ただしpiで規格化済み
  plot(normal_init_phi, 'b');
  hold on 
% phiを赤で表示する ただしpiで規格化済み
  plot(normal_phi,'r');

% オイラー法を抜け出すか否かの処理
  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;
% flgを0に設定してwhile文から抜け出す
  elseif key(1) == 'n'	% 続行
    loop_count = loop_count + loop;
  end
  
end;