function [] = RUN_TerashimaAndKatoMethod()
%RUN_TERASHIMAANDKATOMETHOD
%
% 作成者：宮澤幸希（m-kouki@moegi.waseda.jp）
% 作成日時：2011/07/29

IWswitch = 1;       % 1 なら入力重みを考慮したヒストグラム頂点の推定を行う、0 なら入力数のみを考慮する
calcWeight = 0;     % 【IWswitch = 1 のときのみ】入力重みを計算する際、両端の仮想隣接ノードの重みを同じにするなら 1、ゼロにするなら 0

%正規分布する入力値を作る（引数の詳細は makeDist.m を参照）
[P, Plabel] = makeDist(5, [-1 0 2 5 8], [-4 0 8 -1 2], [3000 5000 1000 1000 1000], 'randn', 0);

% 入力層の要素の範囲を生成
iprange = [min(P') ; max(P')]';   %[最小値 最大値]を特徴量の数だけ行方向に連結
 
%SOMのパラメータを決める
outputSizeX = 30;       %SOMのX軸（横軸、列成分）のセル数
outputSizeY = 1;        %SOMのY軸（縦軸、行成分）のセル数
studyTotalNum = 100;    %順序付けフェイズのステップ数（デフォルトは1000）
reinforceRadius = 1;    %調整フェイズの近傍距離（デフォルトは1）
maxRewardO = 0.9;       %順序付けフェイズの学習比（デフォルトは0.9）
maxRewardT = 0.5;       %調整フェイズの学習比（デフォルトは0.02）
totalEpochs = 1;        %エポック数

%SOM作成（四角形の (outputSizeX)行 (outputSizeY)列ネットワーク）、初期重み値設定
net = newsom(iprange,[outputSizeX outputSizeY],'gridtop','boxdist',maxRewardO, studyTotalNum, maxRewardT, reinforceRadius);
%学習回数の設定
net.trainParam.epochs = totalEpochs;

%初期値をプロット
figure(1);
hold on                         %上のグラフをホールドしたまま、下のグラフも同じ画面にプロット
plotsom(net.IW{1,1},net.layers{1}.distances);    %SOMのベクトルをプロット（赤の点が出力層ノードのベクトル）
text(0.6,0.5,'全ノードの初期位置')
hold off

%学習を開始する
net.trainParam.epochs = totalEpochs;
net = train(net,P);

%評価用 入力ベクトルセットで、各入力の勝利ニューロンを確認
% 学習結果のシミュレーション
Pclass = vec2ind( sim(net,P) )
% 分類結果とネットワークの値をプロット
figure(2);
hold on;
graphlabel = [ '.b' ; '.g' ; '.c' ; '.m' ; '.y' ; '.k' ; '.b' ; '.g' ; '.c' ; '.m' ; '.y' ; '.k' ];
for count = 1 : 1 : outputSizeX
    thisgraphlabel = rem(count,length(graphlabel)) + 1;
    wNsells = find( Pclass == count );
    plot( P(1,wNsells), P(2,wNsells), graphlabel(thisgraphlabel, :) );
end
plotsom(net.IW{1,1},net.layers{1}.distances);    %SOMのベクトルをプロット（赤の点が出力層ノードのベクトル）
title('分類結果 カテゴリ推定前');
hold off;

%密度ヒストグラムによるカテゴリ数推定
% 各ノードに対応付けられた入力ベクトルの数を計算してプロット（ヒストグラム作成）
figure(3);
data = [zeros(length(P), 2) P'];    % Pの前に2列のゼロデータを挿入して unificateSom のデータフォーマットに合わせる
[ net, categoryNum ] = unificateSom2( data, net, 0, IWswitch, calcWeight );
%UNIFICATESOM 学習結果のカテゴリーを統合する（カテゴリー中心以外のノードの重みを 0 にする）
% 第一引数 : 学習用データ、形式は【フレーム番号, 時間, 特徴量（任意の次元数）】
% 第ニ引数 : 学習が終わった SOM のネットワーク構造体（ train の出力）
% 第三引数 : ヒストグラムの平滑化をするなら 1、しないなら 0
% 第四引数 :【任意指定】1 なら入力重みを考慮したヒストグラム頂点の推定を行う、0 なら入力数のみを考慮する
% 第五引数 :【inputWeight を指定したときのみ】入力重みを計算する際、両端の仮想隣接ノードの重みを同じにするなら 1、ゼロにするなら 0

%AICによる統合後クラスターの妥当性の判定
PclassUni = vec2ind( sim(net,P) );
unique(PclassUni)
categoryNode = unique(PclassUni);

% 寺島法で獲得したカテゴリーの隣あう各ペアについて順番にみていく
for count = 1 : 1 : length(categoryNode) - 1;
    category1 = categoryNode(count);    % 今回のペアのひとつめ
    category2 = categoryNode(count+1);  % 今回のペアのふたつめ
    disp( strcat( int2str(category1), ' と ', int2str(category2), 'を一つとみなすか、別々とみなすか' ) );
    % 「統合した一混合GMM」と「別々のニ混合GMM」のどちらが良いか、AICを使って評価する
    %  まずは category1 と category2 に属するデータだけを取り出す
    wNsells1 = find( PclassUni == category1 );
    wNsells2 = find( PclassUni == category2 );
    
    figure(4);
    subplot(2,3,[1 2 4 5]);
    plot(P(1, wNsells1), P(2, wNsells1), '.', P(1, wNsells2), P(2, wNsells2), '.');
    %disp( ' * Enter キーを押して下さい * ' );
    %pause();
    X = [P(:, wNsells1) P(:, wNsells2)]';
    %plot(X(:,1),X(:,2),'.');
    
    %  GMM
    options = statset('Display','final');			% STATS構造体の表示設定
    %  最適なカテゴリ数を AIC or BIC で推定
    AIC = zeros(1,2);                               % AICを格納するセル配列
    BIC = zeros(1,2);                               % BICを格納するセル配列
    for k = 1:1:2
        obj = gmdistribution.fit(X,k,'Options',options);	% カテゴリ（コンポーネント）数kで推定
        %GMMをプロット
        subplot(2,3,(k*3));
        scatter(X(:,1),X(:,2),'.');
        hold on
        h = ezcontour(@(x,y)pdf(obj,[x y]),[-8 6],[-8 6]);	% 確率密度関数を視覚化
        hold off
        AIC(k)= obj.AIC;                            % AICを求める
        BIC(k)= obj.BIC;                            % BICを求める
    end
    AIC
    categoryNumAIC = find(AIC == min(AIC));         % AICが最小になるのは何番目かを調べる
    disp( strcat( 'AIC : カテゴリ数は ', int2str(categoryNumAIC), 'が妥当' ) );
    BIC
    categoryNumBIC = find(BIC == min(BIC));         % BICが最小になるのは何番目かを調べる
    disp( strcat( 'BIC : カテゴリ数は ', int2str(categoryNumBIC), 'が妥当' ) );
    %注：上の推定があまりうまくいかない。完全な正規分布ならカテゴリ数 1 が妥当と正しく判断されるのだが…。
    %    詳細な検討は testGMMandAIC.m を参照。

    disp( ' * Enter キーを押して下さい * ' );
    pause();
end

% 結果をプロット
%figure(4);
%hold on;
%graphlabel = [ '.b' ; '.g' ; '.r' ; '.c' ; '.m' ; '.y' ; '.k' ; '.b' ; '.g' ; '.r' ; '.c' ; '.m' ; '.y' ; '.k' ];
%for count = 1 : 1 : length(categoryNode)
%    thisCategoryNode = categoryNode(count);     % 頂点ノードの番号
%    thisgraphlabel = rem(count,length(graphlabel)) + 1;
%    wNsells = find( PclassUni == thisCategoryNode );
%    plot( P(1,wNsells), P(2,wNsells), graphlabel(thisgraphlabel, :) );
%end
%title('分類結果 カテゴリ推定後');
%hold off;

disp( ' 画面を閉じます... * Enter キーを押して下さい * ' );
pause();
close all;

end