runicatest.m

含有多种ICA算法的eeglab工具箱

       fprintf('Using starting weights named on commandline ...\n');
       fprintf('Returning the identity matrix in variable "sphere" ...\n');
      end
  end
  sphere = eye(chans,chans);
  if verbose,
      fprintf('Returned variable "sphere" will be the identity matrix.\n');
  end
end
%
%%%%%%%%%%%%%%%%%%%%%%%% Initialize ICA training %%%%%%%%%%%%%%%%%%%%%%%%%
%
  lastt=fix((datalength/block-1)*block+1);
  BI=block*eye(ncomps,ncomps);
  delta=zeros(1,chans*ncomps);
  changes = [];
  degconst = 180./pi;
  startweights = weights;
  prevweights = startweights;
  oldweights = startweights;
  prevwtchange = zeros(chans,ncomps);
  oldwtchange = zeros(chans,ncomps);
  lrates = zeros(1,maxsteps);
  onesrow = ones(1,block);
  bias = zeros(ncomps,1);
  signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
  for k=1:nsub
      signs(k) = -1;
  end
  if extended & extblocks < 0 & verbose,
    fprintf('Fixed extended-ICA sign assignments:  ');
    for k=1:ncomps
       fprintf('%d ',signs(k));	
    end; fprintf('\n');
  end
  signs = diag(signs); % make a diagonal matrix
  oldsigns = zeros(size(signs));;
  signcount = 0;              % counter for same-signs
  signcounts = [];
  urextblocks = extblocks;    % original value, for resets
  old_kk = zeros(1,ncomps);   % for kurtosis momemtum
%
%%%%%%%% ICA training loop using the logistic sigmoid %%%%%%%%%%%%%%%%%%%
%
  if verbose,
    fprintf('Beginning ICA training ...');
    if extended,
        fprintf(' first training step may be slow ...\n');
    else
        fprintf('\n');
    end
  end
  step=0;
  laststep=0; 
  blockno = 1;  % running block counter for kurtosis interrupts

  while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % permute=randperm(datalength); % shuffle data order at each step
      bootstrap = round(datalength*rand(1,datalength)+0.5); % draw a new bootstrap dataset at each step

    for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%%
      pause(0);
      if ~isempty(get(0, 'currentfigure')) & strcmp(get(gcf, 'tag'), 'stop')
          close; error('USER ABORT');
      end;
      if biasflag                                                   
       u=weights*data(:,bootstrap(t:t+block-1)) + bias*onesrow;      
      else                                                             
       u=weights*data(:,bootstrap(t:t+block-1));                      
      end                                                              
      if ~extended
       %%%%%%%%%%%%%%%%%%% Logistic ICA weight update %%%%%%%%%%%%%%%%%%%
       y=1./(1+exp(-u));                                                %
       weights=weights+lrate*(BI+(1-2*y)*u')*weights;                   %
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      else % extended-ICA
       %%%%%%%%%%%%%%%%%%% Extended-ICA weight update %%%%%%%%%%%%%%%%%%%
       y=tanh(u);                                                       %
       weights = weights + lrate*(BI-signs*y*u'-u*u')*weights;          %
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      end
      if biasflag 
         if ~extended
       %%%%%%%%%%%%%%%%%%%%%%%% Logistic ICA bias %%%%%%%%%%%%%%%%%%%%%%%
             bias = bias + lrate*sum((1-2*y)')'; % for logistic nonlin. %
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         else % extended
       %%%%%%%%%%%%%%%%%%% Extended-ICA bias %%%%%%%%%%%%%%%%%%%%%%%%%%%%
             bias = bias + lrate*sum((-2*y)')';  % for tanh() nonlin.   %
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         end                                    
      end

      if momentum > 0 %%%%%%%%% Add momentum %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        weights = weights + momentum*prevwtchange;                
        prevwtchange = weights-prevweights;                      
        prevweights = weights;                                  
      end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      if max(max(abs(weights))) > MAX_WEIGHT
          wts_blowup = 1;
          change = nochange;
      end
      if extended & ~wts_blowup
       %
       %%%%%%%%%%% Extended-ICA kurtosis estimation %%%%%%%%%%%%%%%%%%%%%
       %
       if extblocks > 0 & rem(blockno,extblocks) == 0, 
                                  % recompute signs vector using kurtosis
        if kurtsize < frames % 12-22-99 rand() size suggestion by M. Spratling
           rp = fix(rand(1,kurtsize)*datalength);  % pick random subset
           % Accout for the possibility of a 0 generation by rand
           ou = find(rp == 0);
           while ~isempty(ou) % 1-11-00 suggestion by J. Foucher
              rp(ou) = fix(rand(1,length(ou))*datalength);
              ou = find(rp == 0);
           end
           partact=weights*data(:,rp(1:kurtsize));
        else                                        % for small data sets,
            partact=weights*data;                   % use whole data
        end
        m2=mean(partact'.^2).^2; 
        m4= mean(partact'.^4);
        kk= (m4./m2)-3.0;                           % kurtosis estimates
        if extmomentum
            kk = extmomentum*old_kk + (1.0-extmomentum)*kk; % use momentum
            old_kk = kk;
        end
        signs=diag(sign(kk+signsbias));             % pick component signs
        if signs == oldsigns,
           signcount = signcount+1;
        else
           signcount = 0;
        end
        oldsigns = signs;
        signcounts = [signcounts signcount];
        if signcount >= SIGNCOUNT_THRESHOLD,
         extblocks = fix(extblocks * SIGNCOUNT_STEP);% make kurt() estimation
         signcount = 0;                             % less frequent if sign
        end                                         % is not changing
       end % extblocks > 0 & . . .
      end % if extended %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      blockno = blockno + 1;
      if wts_blowup
         break
      end
    end % training block %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if ~wts_blowup
      oldwtchange = weights-oldweights;
      step=step+1; 
      %
      %%%%%%% Compute and print weight and update angle changes %%%%%%%%%
      %
      lrates(1,step) = lrate;
      angledelta=0.;
      delta=reshape(oldwtchange,1,chans*ncomps);
      change=delta*delta'; 
    end
    %
    %%%%%%%%%%%%%%%%%%%%%% Restart if weights blow up %%%%%%%%%%%%%%%%%%%%
    %
    if wts_blowup | isnan(change)|isinf(change),  % if weights blow up,
        fprintf('');
        step = 0;                          % start again
        change = nochange;
        wts_blowup = 0;                    % re-initialize variables
        blockno = 1;
        lrate = lrate*DEFAULT_RESTART_FAC; % with lower learning rate
        weights = startweights;            % and original weight matrix
        oldweights = startweights;            
        change = nochange;
        oldwtchange = zeros(chans,ncomps);
        delta=zeros(1,chans*ncomps);
        olddelta = delta;
        extblocks = urextblocks;
        prevweights = startweights;
        prevwtchange = zeros(chans,ncomps);
        lrates = zeros(1,maxsteps);
        bias = zeros(ncomps,1);
        if extended
           signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
           for k=1:nsub
               signs(k) = -1;
           end
           signs = diag(signs); % make a diagonal matrix
           oldsigns = zeros(size(signs));;
        end
        if lrate> MIN_LRATE
          r = rank(data);
          if r<ncomps
           fprintf('Data has rank %d. Cannot compute %d components.\n',...
                                  r,ncomps);
           return
          else
           fprintf(...
            'Lowering learning rate to %g and starting again.\n',lrate);
          end
        else
          fprintf( ...
          'runica(): QUITTING - weight matrix may not be invertible!\n');
          return;
        end
    else % if weights in bounds 
     %
     %%%%%%%%%%%%% Print weight update information %%%%%%%%%%%%%%%%%%%%%%
     %
     if step> 2 
         angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
     end
     if verbose,
      if step > 2, 
          if ~extended,
              fprintf(...
         'step %d - lrate %5f, wchange %7.6f, angledelta %4.1f deg\n', ...
          step,lrate,change,degconst*angledelta);
          else
              fprintf(...
'step %d - lrate %5f, wchange %7.6f, angledelta %4.1f deg, %d subgauss\n',...
          step,lrate,change,degconst*angledelta,(ncomps-sum(diag(signs)))/2);
          end
      elseif ~extended
         fprintf(...
          'step %d - lrate %5f, wchange %7.6f\n',step,lrate,change);
      else
         fprintf(...
          'step %d - lrate %5f, wchange %7.6f, %d subgauss\n',...
           step,lrate,change,(ncomps-sum(diag(signs)))/2);
      end % step > 2
     end; % if verbose
%
%%%%%%%%%%%%%%%%%%%% Save current values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
      changes = [changes change];
      oldweights = weights;
%
%%%%%%%%%%%%%%%%%%%% Anneal learning rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
      if degconst*angledelta > annealdeg,  
        lrate = lrate*annealstep;          % anneal learning rate
        olddelta   = delta;                % accumulate angledelta until
        oldchange  = change;               %  annealdeg is reached
      elseif step == 1                     % on first step only
        olddelta   = delta;                % initialize 
        oldchange  = change;               
      end
%
%%%%%%%%%%%%%%%%%%%% Apply stopping rule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
      if step >2 & change < nochange,      % apply stopping rule
        laststep=step;            
           step=maxsteps;                  % stop when weights stabilize
      elseif change > DEFAULT_BLOWUP,      % if weights blow up,
        lrate=lrate*DEFAULT_BLOWUP_FAC;    % keep trying 
      end;                                 % with a smaller learning rate
    end; % end if weights in bounds

  end; % end training %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if ~laststep
    laststep = step;
  end;
  lrates = lrates(1,1:laststep);           % truncate lrate history vector
  %
  %%%%%%%%%%%%%% Orient components towards max positive activation %%%%%%
  %
  if strcmp(posactflag,'on')
      [activations,winvout,weights] = posact(data,weights);
       % changes signs of activations and weights to make activations
       % net rms-positive
  else
       activations = weights*data;
  end
  %
  %%%%%%%%%%%%%% If pcaflag, compose PCA and ICA matrices %%%%%%%%%%%%%%%
  %
  if strcmp(pcaflag,'on')
    fprintf('Composing the eigenvector, weights, and sphere matrices\n');
    fprintf('  into a single rectangular weights matrix; sphere=eye(%d)\n'...
                                                                  ,chans);
    weights= weights*sphere*eigenvectors(:,1:ncomps)'; 
    sphere = eye(urchans);
  end
  %
  %%%%%% Sort components in descending order of max projected variance %%%%
  %
  if verbose,
    fprintf(...
   'Sorting components in descending order of mean projected variance ...\n');
  end
  if wts_passed == 0
    %
    %%%%%%%%%%%%%%%%%%%% Find mean variances %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    meanvar  = zeros(ncomps,1);      % size of the projections
    if ncomps == urchans % if weights are square . . .
      winv = inv(weights*sphere);
    else
      fprintf('Using pseudo-inverse of weight matrix to rank order component projections.\n');
      winv = pinv(weights*sphere);
  end
  for s=1:ncomps
    if verbose,
      fprintf('%d ',s);         % construct single-component data matrix
    end
   					% project to scalp, then add row means 
    compproj = winv(:,s)*activations(s,:);
    meanvar(s) = mean(sum(compproj.*compproj)/(size(compproj,1)-1));
                                            % compute mean variance 
  end                                         % at all scalp channels
  if verbose,
   fprintf('\n');
  end
  %
  %%%%%%%%%%%%%% Sort components by mean variance %%%%%%%%%%%%%%%%%%%%%%%%
  %
  [sortvar, windex] = sort(meanvar);
  windex = windex(ncomps:-1:1); % order large to small 
  meanvar = meanvar(windex);
  % 
  %%%%%%%%%%%%%%%%%%%%% Filter data using final weights %%%%%%%%%%%%%%%%%%
  %
  if nargout>2, % if activations are to be returned
   if verbose,
     fprintf('Permuting the activation wave forms ...\n');
   end
   activations = activations(windex,:);
  else
   clear activations
  end
  weights = weights(windex,:);% reorder the weight matrix
  bias  = bias(windex);		% reorder them
  signs = diag(signs);        % vectorize the signs matrix
  signs = signs(windex);      % reorder them
else
  fprintf('Components not ordered by variance.\n');
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    
return

if nargout > 6
  u=weights*data + bias*ones(1,frames);      
  y = zeros(size(u));
  for c=1:chans
    for f=1:frames
       y(c,f) = 1/(1+exp(-u(c,f)));
    end
  end
end