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