function [S,A,loglikelihood,Sigma,chi,exitflag]=ica_adatap(X,prior,par,draw);% ICA_ADATAP  Mean field independent component analysis (ICA)%    [S,A,LL,SIGMA,CHI,EXITFLAG]=ICA_ADATAP(X) performs linear %    instanteneous mixing ICA by solving the adaptive TAP or mean field %    or variational mean field equations [1-4].%   %    Input and output arguments: %      %    X        : Mixed signal%    A        : Estimated mixing matrix%    S        : Estimated source mean values                                   %    SIGMA    : Estimated noise covariance%    LL       : Log likelihood, log P(X|A,SIGMA), divided by  %               number of samples - mean field estimate%    CHI      : Estimated source covariances%    EXITFLAG : 0, no convergence; 1,  maximum number of iterations  %               reached and 2, convergence criteria met.  %%    Outputs are sorted in descending order according to 'energy'%    sum(A.^2,1))'.*sum(S.^2,2).%%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PRIOR) allows for specification %    of priors on A, S, SIGMA and the method to be used. PRIOR.method %    overrides other choices of priors. PRIOR.A, PRIOR.S, Prior.Sigma %    and PRIOR.method are character strings that can take the following %    values%    %    PRIOR.A :            constant      constant mixing matrix %                         free          free likelihood optimization%                         positive      constrained likelihood %                                       optimization (uses QUADPROG)%                         MacKay        corresponding to MacKay's %                                       hyperparameter update%%    PRIOR.S :            bigauss       sum of two Gaussians, default%                                       variance 1 centered at +/-1.  %                         binary        binary +/-1 %                         binary_01     binary 0/1             %                         combi         combinations of other priors%                                       E.g. M1 binary_01's and M2%                                       exponential's are specified by%                                       PRIOR.S1='binary_01';%                                       PRIOR.M1=M1; %                                       PRIOR.S2='exponential'; %                                       PRIOR.M2=M2;%                         exponential   exponential (positive)%                         heavy_tail    heavy_tailed (non-analytic %  				        power law tail)%                         Laplace       Laplace (double exponential)%                         Gauss         Gaussian for factor analysis%                                       PRIOR.Sigma='diagonal' and%                                       probabilistic PCA %                                       PRIOR.Sigma='isotropic'%%    PRIOR.Sigma :        constant      constant noise covariance%                         free          free likelihood optimization %                         isotropic     isotropic noise covariance%                                       (scalar noise variance)%                         diagonal      different noise variance for %                                       each sensor.%%    PRIOR.method :       constant      constant A and Sigma %                                       for test sets. A and Sigma%                                       should initialized, see %                                       below. Sets prior.A='constant'%                                       and prior.Sigma='constant'.%                         fa            factor analysis (FA) sets%                                       prior.A='free', %                                       prior.S='Gauss' and%                                       prior.Sigma='diagonal'.%                         neg_kurtosis  negative kurtosis set%                                       prior.S='bigauss'.%                         pos_kurtosis  positive kurtosis sets%                                       prior.S='heavy_tail'. %                         positive      positive ICA sets%                                       prior.S='exponential' and%                                       prior.A='positive'.%                         ppca          probabilistic PCA (pPCA)%                                       sets prior.A='free', %                                       prior.S='Gauss' and %                                       prior.Sigma='isotropic'.%%    Default values are PRIOR.A='free', PRIOR.S='Laplace',%    PRIOR.Sigma='isotropic' and PRIOR.method not set.%%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PRIOR,PAR) is used for giving %    additional arguments. PAR is a structure array with the %    following fields %   %       sources           Number of sources (default is quadratic%                         size(X,1)). This parameter is overridden %                         by size(PAR.A_init,2) if PAR.A_init is%                         defined.%       solver            method for solving mean field equations%                         can take the following character string%                         values, default is sequential: %%                         beliefprop2   Belief propagation -      %                                       Sequential second order %                                       method by T. Minka. %                         sequential    normal sequential iterative%                                       update (first order method) %                         %       A_init            Initial mixing matrix (default is Toeplitz%                         type with small randmom values).        %       S_init            Initial source values (default is zero).%       Sigma_init        Initial noise covariance (scalar for %                         isotropic). Default is Sigma_rel times%                         empirical covariance.%       Sigma_rel         See above, default 1 for PRIOR.A=%                         'A_positive' and 100 otherwise. %       max_ite           Maximum of A and SIGMA updates (default%                         is 20)%       S_max_ite         Maximum number of S update steps in each%                         E-step (default is 100)%       tol               Termination tolerance for program in %                         terms of relative change in 1-norm of %                         A and Sigma (default is 10^-3).%       S_tol             Termination tolerance for S updates in %                         terms of the squared error of the S%                         mean field equations (default is 10^-10).%%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PRIOR,PAR,DRAW) with DRAW %    different from zero will output runtime information. DRAW%    equal to 2 will output the runtime value of the log %    likelihood. Default value is 1.%%    The memory requirements are O(M^2*N+D*N), where D is the %    number of sensors, N the number of samples (X is D*N) and M is %    the number of sources. The computational complexity is %    O(M^3*N) plus for positive mixing matrix, D times the %    computational complexity of a M-dimensional quadratic %    programming problem.   %%    It is recommended as a preprocessing step to normalize the %    data, e.g. by dividing by the largest element X=X/max(max(X)).%% References [1-4]:%%    Tractable Approximations for Probabilistic Models: The Adaptive %    Thouless-Anderson-Palmer Mean Field Approach%    M. Opper and O. Winther%    Phys. Rev. Lett. 86, 3695-3699 (2001).%%    Adaptive and Self-averaging Thouless-Anderson-Palmer Mean Field %    Theory for Probabilistic Modeling%    M. Opper and O. Winther%    Physical Review E 64, 056131 (2001). % %    Mean Field Approaches to Independent Component Analysis%    Pedro A.d.F.R. H鴍en-S鴕ensen, Ole Winther and Lars Kai Hansen%    Neural Computation 14, 889-918 (2002).%  %    TAP Gibbs Free Energy, Belief Propagation and Sparsity%    L. Csato, M. Opper and O. Winther%    In Advances in Neural Information Processing Systems 14 (NIPS'2001), %    MIT Press (2002). % - by Ole Winther 2002 - IMM, Technical University of Denmark% - http://isp.imm.dtu.dk/staff/winther/% - version 2.2% Uses adaptive TAP convention refs. [1,2] besides in the definition of the mean % function where it uses the ICA-convention of ref. [3]. [D,N]=size(X);                                                % number of sensor, samplestry prior=prior; catch prior=[]; end                          % make sure prior is defined%try par=par; catch par=[]; end                               % make sure par is defined[prior]=method_init(prior);                                   % initialize method[A,A_cmd,A_prior,M]=A_init(prior,par,D);                      % initialize AA_norm_old=norm(abs(A),1);                                    % quantify chance in A[S,S_mean_cmd_all,S_mean_cmd_seq,likelihood_cmd,likelihood_arg,... S_arg_seq,cS_prior,cM,S_prior]=S_init(prior,par,M,N);        % initialize S[Sigma,Sigma_cmd,Sigma_eta,Sigma_prior,Sigma_rel]=...Sigma_init(prior,par,A_prior,S_prior,X,A,S);                  % initialize Sigmadim_Sigma=size(Sigma,1)*size(Sigma,2);                        % number of elements in Sigma Sigma_norm_old=norm(abs(Sigma),1);                            % quantify chance in Sigma.try max_ite=par.max_ite; catch max_ite=20; end                % default maximum number of EM-steps try S_max_ite=par.S_max_ite; catch S_max_ite=100; end         % default maximum number of S updates in each E-step try S_min_ite=par.S_min_ite; catch S_min_ite=1; end           % default minimum number of S updates in each E-step try tol=par.tol; catch tol=10^-3; end                         % termination tolerance relative chance in norm(Sigma)+norm(A) try S_tol=par.S_tol; catch S_tol=10^-10; end; S_tolN=S_tol/N; % termination tolerance chance in dS.*dS try                                                           % default solver is sequential  switch par.solver    case {'beliefprop2','convergent'}      solver=par.solver;        otherwise      solver='sequential';  endcatch solver='sequential'; endtry draw=draw; catch draw=1; end                              % set draw % here starts the algorithm! if strcmp(solver,'beliefprop2')  fprintf(' ICA Adaptive TAP - %s solver\n',solver); else    fprintf(' ICA Linear Response - %s solver\n',solver); endfprintf(' %s source prior\n',S_prior); %s A and %s noise optimization\n',S_prior,A_prior,Sigma_prior); if strcmp(S_prior,'combi')  for i=1:length(cM)     if cM(i)==1 fprintf('   1 %s source prior\n',cS_prior{i});      else fprintf('   %i %s source priors\n',cM(i),cS_prior{i}); end  endendfprintf(' %s A and %s noise optimization\n',A_prior,Sigma_prior); fprintf(' sensors: D = %i\n sources: M = %i\n samples: N = %i\n\n',D,M,N); dS=zeros(M,N); S_old=zeros(M,N);chi=zeros(M,M,N); diagchi=zeros(M,N);G=zeros(M,M,N);diagG=zeros(M,N);hpred=zeros(M,1); V=zeros(M,N); dV=zeros(M,N); traceSS=zeros(M,M); tracechi=zeros(M,M);dA_rel=Inf; dSigma_rel=Inf;if draw   if dim_Sigma==D    fprintf(' noise variance = %g \n\n',sum(Sigma)/D);  else           fprintf(' noise variance = %g \n\n',trace(Sigma)/size(Sigma,1));  end  endEM_step=0;while EM_step<max_ite & (dA_rel>tol | dSigma_rel>tol)   EM_step=EM_step+1;  if dim_Sigma==D    InvSigma=1./Sigma;      for i=1:M  % set up matrices for M-step.      for j=i:M J(i,j)=-(A(:,i))'*(InvSigma.*A(:,j)); J(j,i)=J(i,j); end      for k=1:N h(i,k)=(A(:,i))'*(InvSigma.*X(:,k)); end    end        else    InvSigma=inv(Sigma);     J=-A'*InvSigma*A;    h=A'*InvSigma*X;  end  if EM_step>1 % heuristic for initializing V so that V<V_0    Vrel=repmat(diag(J)./diagJ,1,N); V=Vrel.*V;  else     V=zeros(M,N);  end  diagJ=diag(J);    V_0=-repmat(diag(J),1,N);  J=J-diag(diag(J));   S_ite=0;    dSdS=Inf;  dSdSN=Inf*ones(N,1);  I=1:N;  switch solver     case 'beliefprop2'      alpha=S;           Omega=zeros(M,N);      for i=1:N        G(:,:,i)=J;      end