%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function B = nlem(X)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 非線形EM アルゴリズム(？)で混合データ x を分離する行列 %%%%%
%%%%%% B を求めるプログラム %%%%%
%%%%%% 前処理の部分はCardosoのプログラムから一部拝借 %%%%%
%%%%%% %%%%%
%%%%%% 潜在変数の prior を Gaussianから変更．非線形性が入ることで分 %%%%%
%%%%%% 離が可能になる(はず)．Olshausen & Field のアルゴリズム． %%%%%
%%%%%% 樺島祥介 %%%%%
%%%%%% checked & modified by 池田思朗 %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% データ行列のサイズを調べる
[N,P]=size(X);
% 定行列の設定
Id=eye(N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% 前処理 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% 各行の平均値を引く %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MX=mean(X');
X=X-mean(X')'*ones(1,P);
%%%%%%%%%%%%%%%%%% ホワイトニング %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[U,D] = eig((X*X')/P) ; % 共分散行列の固有値と固有ベクトル
% U: 固有ベクトル D: 固有値
[puiss,k] = sort(diag(D)) ; % 固有値の大きさ順に並べ替え
rangeW = 1:N ;
scales = sqrt(puiss(rangeW)); % 固有値の平方根
W = diag(1./scales) * U(1:N,k(rangeW))';
% W=固有値の平方根で規格化＋固有値基底への変換を行うホワイトニング行列
X = W*X; % ホワイトニング実行
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% 学習 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global C D Lambda Sigma; % 潜在変数(活動度 a )を求める内ループ(in_loop)に
%%%%%%%%%%%%%%%%%%%% 受け渡すグローバル変数の宣言
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 分離行列の初期条件
B=eye(N); %単位行列に選ぶ
L=0.01; %学習率
Sigma=0.01; %活動度の特徴的範囲．適当に調整．
Lambda=0.5; % Prior の強さ． 適当に調整．
% 活動度ベクトルのダイナミクスを解く必要があるので学習例を 1000 に
% 減らしている．
for t=1:1000
u=X(:,t);
% u に対する潜在変数(活動度ベクトル)を求める．方程式はin_loop.m 参照．
C=B'*u;
D=B'*B;
a=fmins('in_loop',(D\C));
% ヘッブ学習．ただし，直交行列に制限して．modified by shiro
B=(Id-L*((u-B*a)*a'-a*(u-B*a)'))*B;
end;
% ホワイトニングを含む転置行列に変換
B = B'*W ;
return;