% BoostLogPost20.m   2004/6/12
%
% DataGen20.mat を読んでテストデータの判別結果を図にする
%
% 1. LDF で全画素の事後確率を推定
% 2. 判別関数 F(x,y) による事後確率
%    p = exp(F(x,y))/sum(F(x,1:g))  と事後確率を定義
% 3. sum(log p_1 + ... + log p_n) を最大にする F を探索
% 4. 半径 r の円環状の近傍への重みを決定
% 5. これらを組み合わせて判別基準を作成
%
% 入力画像: DataGen20.mat , ... , DataGen23.mat
%
% 教師データの初期重み uniform = 1: 各点に 1/gg
%                      uniform = 0: カテゴリ g の点に 1/n_g
%
% 円環状の近傍は
% C:\MyDocuments\MatlabDat\NewToolBox\NishiiOriginal\NBDcircle.m
% が作成

format compact; clear; close all; diary off;
global bb gg cat lab n mis correct onen
global NBDi NBDj NBDsize ssi ssj len

delete BoostLogPost20.dia; diary BoostLogPost20.dia;

% 半径の2乗の動く範囲
nbdr2set = [1 2 4 5 8 9 10 13 16 17 18 20 25 26 29 32 34 36 37 40 45 49 50];
lennbd   = length(nbdr2set);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 教師データ, テストデータの読み込み
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load DataGen20.mat

bb = 4; gg = 3; cat = zeros(1, gg); pp = bb;
col  = [  1   0  0 ;  1 1 1;   0  1   0; 0 0 1; 1   1   0;   0 1   1 ; 1 0 1; ...
        1/2 0 0; 0 1/2 0; 0 0 1/2; 1/2 1/2 0; 0 1/2 1/2; 1/2 0 1/2];
mono = [  0 0 0 ; 1 1 1; 1/2 1/2 1/2];

test      = TestImage(:, :, :, 1);                % 画素(i,j)での多重分光画像
train     = TestImage(:, :, :, 2);                % 画素(i,j)での多重分光画像
datatrain = reshape(train, lenv*lenh, bb);        % 画素を縦に並べた train
datatest  = reshape(test , lenv*lenh, bb);        % 画素を縦に並べた test
clear TestImage

Ztest  = Ztest( 1:lenv, 1:lenh);          % 教師ラベル行列 ＝ テストラベル行列
Ztrain = Ztest;                           % ラベル行列

for g = 1:gg
 cat(g) = length(find(Ztrain(:) == g));
end % g
n = sum(cat);
fprintf('テストデータの各カテゴリの標本数 = ( %s) , 総標本数 = %d\n', ...
         num2str(cat), n)

onen = ones( n, 1);           % n x 1 の1ベクトル
lab  = zeros(n, gg);          % lab : 各グループへの属性を表す行列
tlab = zeros(n,  1);          % tlab: 各グループへの属性を表すベクトル
for g = 1:gg
 tmp  = find( Ztrain(:) == g ); lab(tmp, g) = 1;
 tlab = tlab + g * lab(:, g);
end % g


for uniform = 1:1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 教師データの初期重み: 一様 or 各グループの重みの和が 1/G
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if uniform == 1
 weight0 = ones(n, gg) / gg;
 tmp = find(Ztrain(:) > 0); weight0 = weight0(tmp, :); clear tmp lab
                              % データのない画素番号は飛ばす
 fprintf('\n\n*** 教師データの初期重み = 1/gg = 一様 ***\n')
else
% weight0 = zeros(n, gg);
% for y = 1:gg
%  for ye = [1:y-1 y+1:gg]
%   weight0(find(Ztrain(:) == y), ye) = 1/(gg*(gg-1)*cat(y));
%  end
% end
 fprintf('*** 教師データの初期重み = 各グループでの和が 1/G ***\n')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LDF を構成し，教師データ，テストデータを判別
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu = zeros(bb, gg); C = zeros(bb);
for g = 1:gg
 tmp  = find( Ztrain(:) == g );
 mu(:, g) = ( mean( datatrain(tmp, :) ) )';
 C = C + (length(tmp)-1) * cov( datatrain(tmp, :) );
end % g
C = C / n; Ci = inv(C);
fprintf('教師データにより推定された各群の平均ベクトル\n')
mu

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LDFtrainPost: n x gg, LDFtestPost: n x gg  事後確率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
LDFtrain = tlab; LDFtest = tlab; tmp1 = zeros(gg, 1); tmp2 = zeros(gg, 1);
LDFtrainPost = zeros(gg, n); LDFtestPost = zeros(gg, n);
for s = 1:n
 for g = 1:gg
  tmp1(g) = (datatrain(s, :) - mu(:, g)') * Ci * (datatrain(s, :) - mu(:, g)')';
  tmp2(g) = (datatest( s, :) - mu(:, g)') * Ci * (datatest( s, :) - mu(:, g)')';
 end
 [dmy, LDFtrain(s)] = min(tmp1); [dmy, LDFtest(s)] = min(tmp2);  % 画素 s を判別
 tmp1 = exp( -tmp1/2 ); tmp2 = exp( -tmp2/2 );
 LDFtrainPost(:, s) = log( tmp1 / sum( tmp1 ) );  % 事後確率の対数を計算
 LDFtestPost( :, s) = log( tmp2 / sum( tmp2 ) );  % 事後確率の対数を計算
end % s                       % LDF による教師データ，テストデータの誤判別率
LDFtrainPost = LDFtrainPost'; LDFtestPost = LDFtestPost';

betaNBD  = zeros(lennbd+1, 1);% 近傍の事後確率の信頼度を表す係数
                              % LDF および 近傍を用いたときの誤判別率を保存
errtrain = zeros(lennbd+1, 1); errtrain(1) = length(find(tlab ~= LDFtrain))/n;            errtest  = zeros(lennbd+1, 1); errtest( 1) = length(find(tlab ~= LDFtest ))/n;

FUSIONtrain = zeros(n, 1); FUSIONtest = zeros(n, 1);
tmp1 = LDFtrainPost; tmp2 = LDFtestPost;
for s = 1:n
 [dmy, FUSIONtrain(s)] = max(tmp1(s, :));
 [dmy, FUSIONtest(s) ] = max(tmp2(s, :));
end % s
fprintf('\nLDF    classification, radius^2 of NBD = 0, Training error = %f, Test error = %f\n', ...
        errtrain(1), errtest(1) )

%%%%%%%%%%%%%%%%%%%% test data の判別結果を画像にする %%%%%%%%%%%%%%%%%%%%
CC = reshape(FUSIONtest, [lenv, lenh]); nbdr2 = 0;
fname = ['testCOLORr2_', num2str(nbdr2), '.fig'];
figure, imshow(CC, colormap(col));
title('Test Data 20'); saveas(gcf, fname)
fname = ['testMONOr2_', num2str(nbdr2), '.fig'];
figure, imshow(CC, colormap(mono));
title('Test Data 20'); saveas(gcf, fname)
fname = ['testMONOepsr2_', num2str(nbdr2), '.eps'];
figure, imshow(CC, colormap(mono)); print(gcf, '-deps', fname)
pause(0.1); close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 重み weight0 で real adaboost による最適な beta と更新された重みを計算
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[betaLDF, weightLDF] = OptBetaLogPost(tlab, LDFtrainPost', weight0);

fprintf('LDF による事後確率の重み = %f\n\n', betaLDF)
betaNBD(1) = betaLDF;         % 先頭は LDF の信頼度

weightNBD = weightLDF;  % 教師データの LDFで判別した後の重みで逐次的に更新

LDFtrainPost2dim = reshape(LDFtrainPost', [gg, lenv, lenh]);% gg x lenv x lenh
LDFtestPost2dim  = reshape(LDFtestPost' , [gg, lenv, lenh]);
tmp1 = betaLDF*LDFtrainPost;  % 判別機，初期は線形判別，さらに
tmp2 = betaLDF*LDFtestPost;   %         近傍を大きくしていく度に更新

countnbd = 1;
for nbdr2 = nbdr2set;
 tic
 fprintf('円環状の近傍，radius^2 = %d, 教師データの初期重み uniform = %d\n', ...
          nbdr2, uniform)
 %%%%% 半径^2 = nbdr2-1 から nbdr2 までのリング状の近傍を生成 %%%%%
 NBDcircle(reshape(tlab, [lenv, lenh]), 'LDFtrainNBD', nbdr2);
 fprintf('\n円環状の近傍での事後確率への信頼度の計算\n')
 countnbd = countnbd + 1;

 % 近傍画素の事後確率により中心画素の事後確率を計算
 NBDtrainPost = zeros(gg, n); NBDtestPost = zeros(gg, n);
 for s = 1:len
  sumposttrain = zeros(gg, 1); sumposttest = zeros(gg, 1);
  for k = 1:NBDsize(s)
   sumposttrain = sumposttrain + LDFtrainPost2dim(:, NBDi(k, s), NBDj(k, s) );
   sumposttest  = sumposttest  + LDFtestPost2dim( :, NBDi(k, s), NBDj(k, s) );
  end % k
  NBDtrainPost(:, s) = sumposttrain/NBDsize(s);   % 平均した対数事後確率
  NBDtestPost( :, s) = sumposttest /NBDsize(s);   % 平均した対数事後確率
 end % s
 NBDtrainPost = NBDtrainPost'; NBDtestPost = NBDtestPost';
 [betaNBD(countnbd), weightNBD] = OptBetaLogPost(tlab, NBDtrainPost', weightLDF);

 tmp1 = tmp1 + betaNBD(countnbd) * NBDtrainPost;
 tmp2 = tmp2 + betaNBD(countnbd) * NBDtestPost;
 for s = 1:n
  [dmy, FUSIONtrain(s)] = max(tmp1(s, :));
  [dmy, FUSIONtest( s)] = max(tmp2(s, :));
 end % s
 errtrain(countnbd) = length(find(tlab ~= FUSIONtrain))/n;
 errtest( countnbd) = length(find(tlab ~= FUSIONtest ))/n;
 fprintf('Fusion classification, radius^2 of NBD = %d, Training error = %f, Test error = %f\n\n', ...
        nbdr2, errtrain(countnbd) , errtest(countnbd) )

 %%%%%%%%%%%%%%%%%%%% test data の判別結果を画像にする %%%%%%%%%%%%%%%%%%%%
 CC = reshape(FUSIONtest, [lenv, lenh]);
 fname = ['testCOLORr2_', num2str(nbdr2), '.fig'];
 figure, imshow(CC, colormap(col));
 title('Test Data 20'); saveas(gcf, fname)
 fname = ['testMONOr2_', num2str(nbdr2), '.fig'];
 figure, imshow(CC, colormap(mono));
 title('Test Data 20'); saveas(gcf, fname)
fname = ['testMONOepsr2_', num2str(nbdr2), '.eps'];
figure, imshow(CC, colormap(mono)); print(gcf, '-deps', fname)
 pause(0.1); close all
 toc
end % nbdr2
end % uniform


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 事後確率の平均の係数，トレイニングエラー，テストエラーのグラフ化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure, plot([0 nbdr2set], betaNBD')
title('coefficients for LDF and neighborhoods')
saveas(gcf, 'beta.fig')
pause(1)

figure, plot([0 nbdr2set], errtrain)
title('training error, uniform = 1')
saveas(gcf, 'errtrain.fig')
pause(1)

figure, plot([0 nbdr2set], errtest)
title('test error, uniform = 1')
saveas(gcf, 'errtest.fig')
pause(1)

%save RealboostNBD
diary off
