%% 補助点を格子状に配置した補助変数法 (KISS-GP) のデモ
% 以下のコードは GPML のデモと解説を行う文書
% http://www.gaussianprocess.org/gpml/code/matlab/doc/
% の Section 4c を改変したものです.
% GPML-MATLAB は MATLAB によるガウス過程回帰関連ツールボックスです。
% バージョンは gpml-matlab-v4.1-2017-10-19 に基づいています。

% gpml をインストールしたディレクトリにて
% startup.m を実行して path を通したうえで以下を実行してください。

% 最初に、訓練用データを生成します。
% 訓練用データは2次元空間上の  [-2,2]x[-3,3] の領域に
% 格子状に配置されている 15600 pixel のスカラー値です。
% 15600点ともなると通常の GP でまともに取り扱うことは計算量的に
% 困難です。

 x1 = linspace(-2,2,120); x2 = linspace(-3,3,130);
  x = apxGrid('expand',{x1',x2'});
  y = sin(x(:,2)) + x(:,1) + 0.1*gpml_randn(1,[size(x,1),1]);

cc = [min(y),max(y)]*1.5;
figure(1)
imagesc( x1, x2, reshape(y,120,130), cc )
title('Observations on 2D grid (120x130=15600)')
% title('2次元グリッド(120×130=15600点)上の観測データ')
colormap hsv
colorbar

% 次に、テスト出力を与える点を定義します。
% テスト出力は [-4,4]x[-6,6] の領域上の格子点とします。
% 訓練用データ点の外側における補外結果を見るために、ひろい領域を設定してあります。

xs1 = linspace(-4,4,100); xs2 = linspace(-6,6,110);
 xs = apxGrid('expand',{xs1',xs2'});
 
% いよいよ、KISS GP 用に補助点の格子状配置を定義します。
% 補助点は2次元空間上で、訓練用データとテスト用データの分布する
% 範囲をじゅうぶんにカバーするように配置します。

 xg = apxGrid('create',[x;xs],true,[50,70]);

% 訓練データ点、テスト入力点、補助点の配置を可視化する図です。
% 色温度画像は訓練用データ、
% 赤点はテスト用データ点、
% 青色グリッド (50×70) は補助点の配置を示しています。

figure(2)
imagesc( x1, x2, reshape(y,120,130), cc )
title('Observations on 2D grid (120x130=15600)')
% title('2次元グリッド(120×130=15600点)上の観測データ')
colormap hsv
hold on
for i=1:length(xg{1})
    x000 = xg{1}(i);
    plot( x000+[0,0], [max(xg{2}),min(xg{2})], 'b-')
end
for i=1:length(xg{2})
    y000 = xg{2}(i);
    plot( [max(xg{1}),min(xg{1})], y000+[0,0],  'b-')
end
plot( xs(:,1), xs(:,2), 'r.')
axis([-5,5,-8,8])

% ここまで、データセットの定義を行いました。
% ここから、GPモデルとハイパーパラメタの初期設定を行います。
% まず、共分散関数 cov を定義します。 
% 次に格子状配置 xg のもとで共分散を計算する関数 apxGrid
% を使って共分散関数 cov をラッピングします。
% covg は apxGrid に cov と xg を適用したものであり、
% cov の代わりに使うと GP回帰のなかで効率的な共分散行列計算をしてくれます。

% ここで covSEiso は等方的な距離を用いた単純がガウスカーネルを意味します。
  cov = {{@covSEiso},{@covSEiso}}; covg = {@apxGrid,cov,xg};
  mean = {@meanZero}; lik = {@likGauss};
  hyp.cov = zeros(4,1); hyp.mean = []; hyp.lik = log(0.1);

% ハイパーパラメタの最適化を行います。
% 通常の GP モデルの場合と同様に infGaussLik という関数を用います。
% 繰り返し回数の最大値や収束判定用のしきい値などを事前に決めておきます。
% 共役勾配法を適用しています。

  opt.cg_maxit = 200; opt.cg_tol = 5e-3;
  infg = @(varargin) infGaussLik(varargin{:},opt);
  hyp = minimize(hyp,@gp,-50,infg,[],covg,[],x,y);

% 最後に、グリッド上の補間によって予測値を出力します。
% infGrid クラスの predict メソッドでは
% 格子状配置を用いた高速な計算をします。

  [post,nlZ,dnlZ] = infGrid(hyp,{@meanZero},covg,{@likGauss},x,y,opt);
  [fm,fs2,ym,ys2] = post.predict(xs);

% 結果を図示します。
figure(3)
imagesc( xs1, xs2, reshape(fm, 100,110), cc )
title('Estimated output on test locations')
% title( 'テスト位置における出力推定値' )
hold on
colormap hsv
plot(2*[1,-1,-1,1,1],3*[1,1,-1,-1,1],'k-')
  
% 図のファイル出力
for i=1:3
    figure(i)
    fg = gcf;
    fg.Position = [50, 70, 600,720];
    fg.PaperPositionMode = 'auto';
    saveas(gcf, sprintf('kissgp-%d.png',i),'png')
end
