% runica() - Perform Independent Component Analysis (ICA) decomposition % of input data using the logistic infomax ICA algorithm of % Bell & Sejnowski (1995) with the natural gradient feature % of Amari, Cichocki & Yang, or optionally the extended-ICA % algorithm of Lee, Girolami & Sejnowski, with optional PCA % dimension reduction. Annealing based on weight changes is % used to automate the separation process. % Usage: % >> [weights,sphere] = runica(data); % train using defaults % else % >> [weights,sphere,compvars,bias,signs,lrates,activations] ... % = runica(data,'Key1',Value1',...); % Input: % data = input data (chans,frames*epochs). % Note that if data consists of multiple discontinuous epochs, % each epoch should be separately baseline-zero'd using % >> data = rmbase(data,frames,basevector); % % Optional keywords [argument]: % 'extended' = [N] perform tanh() "extended-ICA" with sign estimation % N training blocks. If N > 0, automatically estimate the % number of sub-Gaussian sources. If N < 0, fix number of % sub-Gaussian comps to -N [faster than N>0] (default|0 -> off) % 'pca' = [N] decompose a principal component (default -> 0=off) % subspace of the data. Value is the number of PCs to retain. % 'sphering' = ['on'/'off'] flag sphering of data (default -> 'on') % 'weights' = [W] initial weight matrix (default -> eye()) % (Note: if 'sphering' 'off', default -> spher()) % 'lrate' = [rate] initial ICA learning rate (<< 1) (default -> heuristic) % 'block' = [N] ICA block size (<< datalength) (default -> heuristic) % 'anneal' = annealing constant (0,1] (defaults -> 0.90, or 0.98, extended) % controls speed of convergence % 'annealdeg' = [N] degrees weight change for annealing (default -> 70) % 'stop' = [f] stop training when weight-change < this (default -> 1e-6 % if less than 33 channel and 1E-7 otherwise) % 'maxsteps' = [N] max number of ICA training steps (default -> 512) % 'bias' = ['on'/'off'] perform bias adjustment (default -> 'on') % 'momentum' = [0 0) % 'specgram' = [srate loHz hiHz frames winframes] decompose a complex time/frequency % transform of the data (Note: winframes must divide frames) % (defaults [srate 0 srate/2 size(data,2) size(data,2)]) % 'posact' = make all component activations net-positive(default 'off'} % Requires time and memory; posact() may be applied separately. % 'ncomps' = [N] number of ICA components to compute (default -> chans or 'pca' arg) % using rectangular ICA decomposition. This parameter may return % strange results. This is because the weight matrix is rectangular % instead of being square. Do not use except to try to fix the problem. % 'verbose' = give ascii messages ('on'/'off') (default -> 'on') % % Outputs: [Note: RO means output in reverse order of projected mean variance % unless starting weight matrix passed ('weights' above)] % weights = ICA weight matrix (comps,chans) [RO] % sphere = data sphering matrix (chans,chans) = spher(data) % Note that unmixing_matrix = weights*sphere {if sphering off -> eye(chans)} % compvars = back-projected component variances [RO] % bias = vector of final (ncomps) online bias [RO] (default = zeros()) % signs = extended-ICA signs for components [RO] (default = ones()) % [ -1 = sub-Gaussian; 1 = super-Gaussian] % lrates = vector of learning rates used at each training step [RO] % activations = activation time courses of the output components (ncomps,frames*epochs) % % Authors: Scott Makeig with contributions from Tony Bell, Te-Won Lee, % Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky, Delorme Arnaud, % CNL/The Salk Institute, La Jolla, 1996- % Uses: posact() % Reference (please cite): % % Makeig, S., Bell, A.J., Jung, T-P and Sejnowski, T.J., % "Independent component analysis of electroencephalographic data," % In: D. Touretzky, M. Mozer and M. Hasselmo (Eds). Advances in Neural % Information Processing Systems 8:145-151, MIT Press, Cambridge, MA (1996). % % Toolbox Citation: % % Makeig, Scott et al. "EEGLAB: ICA Toolbox for Psychophysiological Research". % WWW Site, Swartz Center for Computational Neuroscience, Institute of Neural % Computation, University of San Diego California % , 2000. [World Wide Web Publication]. % % For more information: % http://www.sccn.ucsd.edu/eeglab/icafaq.html - FAQ on ICA/EEG % http://www.sccn.ucsd.edu/eeglab/icabib.html - mss. on ICA & biosignals % http://www.cnl.salk.edu/~tony/ica.html - math. mss. on ICA % Copyright (C) 1996 Scott Makeig et al, SCCN/INC/UCSD, scott@sccn.ucsd.edu % % This program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA % $Log: runica.m,v $ % Revision 1.29 2007/02/20 14:18:56 arno % nothing % % Revision 1.28 2006/06/20 02:00:57 toby % Modified for speed, about 15% faster % % Revision 1.27 2006/01/17 16:55:21 scott % re-coded 'posact' processing to be efficient -sm % % Revision 1.26 2006/01/17 16:29:05 scott % re-wrote 'posact' processing for efficiency -sm % % Revision 1.25 2005/10/11 16:17:07 arno % ncomps warning % % Revision 1.24 2005/07/17 15:03:38 scott % shortened computation of cmponent variances % overwrite data with activations to save memory % add back rowmeans when computing activations % change default 'posact' to 'off' to save memory -sm % % Revision 1.23 2004/10/06 20:41:10 scott % rm debug -sm % % Revision 1.22 2004/10/06 20:40:32 scott % same -sm % % Revision 1.21 2004/10/06 20:39:20 scott % fixed block size to integer -sm % % Revision 1.20 2004/06/29 16:44:09 scott % randomize data shuffling by clock state % % Revision 1.19 2004/05/16 01:13:19 scott % typo % % Revision 1.18 2004/05/16 01:09:30 scott % typo % % Revision 1.17 2004/05/16 01:07:44 scott % disable new buggy fprintf (eval(ps)) feature - for emergency use -sm % % Revision 1.16 2004/05/09 22:49:45 scott % made wchange printout hold enough decimal places when 'stop' is small % % Revision 1.15 2004/05/09 18:18:36 scott % added printout of frames per weight -sm % % Revision 1.14 2003/12/15 23:28:34 arno % debug nochangeupdate % % Revision 1.13 2003/12/11 17:51:11 arno % stoping rule debug if more than 32 channels % % Revision 1.12 2003/10/23 15:48:45 arno % indents % % Revision 1.11 2003/10/19 17:14:15 scott % cosmetics, help msg % % Revision 1.10 2003/10/03 18:21:25 arno % releasing constraint pca 1 & rem(nargin,2) == 0) fprintf('runica(): Even number of input arguments???') return end for i = 3:2:nargin % for each Keyword Keyword = eval(['p',int2str((i-3)/2 +1)]); Value = eval(['v',int2str((i-3)/2 +1)]); if ~isstr(Keyword) fprintf('runica(): keywords must be strings') return end Keyword = lower(Keyword); % convert upper or mixed case to lower if strcmp(Keyword,'weights') | strcmp(Keyword,'weight') if isstr(Value) fprintf(... 'runica(): weights value must be a weight matrix or sphere') return else weights = Value; wts_passed =1; end elseif strcmp(Keyword,'ncomps') if isstr(Value) fprintf('runica(): ncomps value must be an integer') return end if ncomps < urchans & ncomps ~= Value fprintf('runica(): Use either PCA or ICA dimension reduction'); return end fprintf('*****************************************************************************************'); fprintf('************** WARNING: NCOMPS OPTION OFTEN DOES NOT RETURN ACCURATE RESULTS ************'); fprintf('************** WARNING: IF YOU FIND THE PROBLEM, PLEASE LET US KNOW ************'); fprintf('*****************************************************************************************'); ncomps = Value; if ~ncomps, ncomps = chans; end elseif strcmp(Keyword,'pca') if ncomps < urchans & ncomps ~= Value fprintf('runica(): Use either PCA or ICA dimension reduction'); return end if isstr(Value) fprintf(... 'runica(): pca value should be the number of principal components to retain') return end pcaflag = 'on'; ncomps = Value; if ncomps > chans | ncomps < 1, fprintf('runica(): pca value must be in range [1,%d]\n',chans) return end chans = ncomps; elseif strcmp(Keyword,'posact') if ~isstr(Value) fprintf('runica(): posact value must be on or off') return else Value = lower(Value); if ~strcmp(Value,'on') & ~strcmp(Value,'off'), fprintf('runica(): posact value must be on or off') return end posactflag = Value; end elseif strcmp(Keyword,'lrate') if isstr(Value) fprintf('runica(): lrate value must be a number') return end lrate = Value; if lrate>MAX_LRATE | lrate <0, fprintf('runica(): lrate value is out of bounds'); return end if ~lrate, lrate = DEFAULT_LRATE; end elseif strcmp(Keyword,'block') | strcmp(Keyword,'blocksize') if isstr(Value) fprintf('runica(): block size value must be a number') return end block = floor(Value); if ~block, block = DEFAULT_BLOCK; end elseif strcmp(Keyword,'stop') | strcmp(Keyword,'nochange') ... | strcmp(Keyword,'stopping') if isstr(Value) fprintf('runica(): stop wchange value must be a number') return end nochange = Value; elseif strcmp(Keyword,'maxsteps') | strcmp(Keyword,'steps') if isstr(Value) fprintf('runica(): maxsteps value must be an integer') return end maxsteps = Value; if ~maxsteps, maxsteps = DEFAULT_MAXSTEPS; end if maxsteps < 0 fprintf('runica(): maxsteps value (%d) must be a positive integer',maxsteps) return end elseif strcmp(Keyword,'anneal') | strcmp(Keyword,'annealstep') if isstr(Value) fprintf('runica(): anneal step value (%2.4f) must be a number (0,1)',Value) return end annealstep = Value; if annealstep <=0 | annealstep > 1, fprintf('runica(): anneal step value (%2.4f) must be (0,1]',annealstep) return end elseif strcmp(Keyword,'annealdeg') | strcmp(Keyword,'degrees') if isstr(Value) fprintf('runica(): annealdeg value must be a number') return end annealdeg = Value; if ~annealdeg, annealdeg = DEFAULT_ANNEALDEG; elseif annealdeg > 180 | annealdeg < 0 fprintf('runica(): annealdeg (%3.1f) is out of bounds [0,180]',... annealdeg); return end elseif strcmp(Keyword,'momentum') if isstr(Value) fprintf('runica(): momentum value must be a number') return end momentum = Value; if momentum > 1.0 | momentum < 0 fprintf('runica(): momentum value is out of bounds [0,1]') return end elseif strcmp(Keyword,'sphering') | strcmp(Keyword,'sphereing') ... | strcmp(Keyword,'sphere') if ~isstr(Value) fprintf('runica(): sphering value must be on, off, or none') return else Value = lower(Value); if ~strcmp(Value,'on') & ~strcmp(Value,'off') & ~strcmp(Value,'none'), fprintf('runica(): sphering value must be on or off') return end sphering = Value; end elseif strcmp(Keyword,'bias') if ~isstr(Value) fprintf('runica(): bias value must be on or off') return else Value = lower(Value); if strcmp(Value,'on') biasflag = 1; elseif strcmp(Value,'off'), biasflag = 0; else fprintf('runica(): bias value must be on or off') return end end elseif strcmp(Keyword,'specgram') | strcmp(Keyword,'spec') if ~exist('specgram') < 2 % if ~exist or defined workspace variable fprintf(... 'runica(): MATLAB Sig. Proc. Toolbox function "specgram" not found.\n') return end if isstr(Value) fprintf('runica(): specgram argument must be a vector') return end srate = Value(1); if (srate < 0) fprintf('runica(): specgram srate (%4.1f) must be >=0',srate) return end if length(Value)>1 loHz = Value(2); if (loHz < 0 | loHz > srate/2) fprintf('runica(): specgram loHz must be >=0 and <= srate/2 (%4.1f)',srate/2) return end else loHz = 0; % default end if length(Value)>2 hiHz = Value(3); if (hiHz < loHz | hiHz > srate/2) fprintf('runica(): specgram hiHz must be >=loHz (%4.1f) and <= srate/2 (%4.1f)',loHz,srate/2) return end else hiHz = srate/2; % default end if length(Value)>3 Hzframes = Value(5); if (Hzframes<0 | Hzframes > size(data,2)) fprintf('runica(): specgram frames must be >=0 and <= data length (%d)',size(data,2)) return end else Hzframes = size(data,2); % default end if length(Value)>4 Hzwinlen = Value(4); if rem(Hzframes,Hzwinlen) % if winlen doesn't divide frames fprintf('runica(): specgram Hzinc must divide frames (%d)',Hzframes) return end else Hzwinlen = Hzframes; % default end Specgramflag = 1; % set flag to perform specgram() elseif strcmp(Keyword,'extended') | strcmp(Keyword,'extend') if isstr(Value) fprintf('runica(): extended value must be an integer (+/-)') return else extended = 1; % turn on extended-ICA extblocks = fix(Value); % number of blocks per kurt() compute if extblocks < 0 nsub = -1*fix(extblocks); % fix this many sub-Gauss comps elseif ~extblocks, extended = 0; % turn extended-ICA off elseif kurtsize>frames, % length of kurtosis calculation kurtsize = frames; if kurtsize < MIN_KURTSIZE fprintf(... 'runica() warning: kurtosis values inexact for << %d points.\n',... MIN_KURTSIZE); end end end elseif strcmp(Keyword,'verbose') if ~isstr(Value) fprintf('runica(): verbose flag value must be on or off') return elseif strcmp(Value,'on'), verbose = 1; elseif strcmp(Value,'off'), verbose = 0; else fprintf('runica(): verbose flag value must be on or off') return end else fprintf('runica(): unknown flag') return end end % %%%%%%%%%%%%%%%%%%%%%%%% Initialize weights, etc. %%%%%%%%%%%%%%%%%%%%%%%% % if ~annealstep, if ~extended, annealstep = DEFAULT_ANNEALSTEP; % defaults defined above else annealstep = DEFAULT_EXTANNEAL; % defaults defined above end end % else use annealstep from commandline if ~annealdeg, annealdeg = DEFAULT_ANNEALDEG - momentum*90; % heuristic if annealdeg < 0, annealdeg = 0; end end if ncomps > chans | ncomps < 1 fprintf('runica(): number of components must be 1 to %d.\n',chans); return end if weights ~= 0, % initialize weights % starting weights are being passed to runica() from the commandline if verbose, fprintf('Using starting weight matrix named in argument list ...\n') end if chans>ncomps & weights ~=0, [r,c]=size(weights); if r~=ncomps | c~=chans, fprintf(... 'runica(): weight matrix must have %d rows, %d columns.\n', ... chans,ncomps); return; end end end; % %%%%%%%%%%%%%%%%%%%%% Check keyword values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % if frames frames, fprintf('runica(): block size exceeds data length!\n'); return elseif floor(epochs) ~= epochs, fprintf('runica(): data length is not a multiple of the epoch length!\n'); return elseif nsub > ncomps fprintf('runica(): there can be at most %d sub-Gaussian components!\n',ncomps); return end; % % adjust nochange if necessary % if isnan(nochange) if ncomps > 32 nochange = 1E-7; nochangeupdated = 1; % for fprinting purposes else nochangeupdated = 1; % for fprinting purposes nochange = DEFAULT_STOP; end; else nochangeupdated = 0; end; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Process the data %%%%%%%%%%%%%%%%%%%%%%%%%% % if verbose, fprintf( ... '\nInput data size [%d,%d] = %d channels, %d frames/n', ... chans,frames,chans,frames); if strcmp(pcaflag,'on') fprintf('After PCA dimension reduction,\n finding '); else fprintf('Finding '); end if ~extended fprintf('%d ICA components using logistic ICA.\n',ncomps); else % if extended fprintf('%d ICA components using extended ICA.\n',ncomps); if extblocks > 0 fprintf(... 'Kurtosis will be calculated initially every %d blocks using %d data points.\n',... extblocks, kurtsize); else fprintf(... 'Kurtosis will not be calculated. Exactly %d sub-Gaussian components assumed.\n',... nsub); end end fprintf('Decomposing %d frames per ICA weight ((%d)^2 = %d weights, %d frames)\n',... floor(frames/ncomps.^2),ncomps.^2,frames); fprintf('Initial learning rate will be %g, block size %d.\n',... lrate,block); if momentum>0, fprintf('Momentum will be %g.\n',momentum); end fprintf( ... 'Learning rate will be multiplied by %g whenever angledelta >= %g deg.\n', ... annealstep,annealdeg); if nochangeupdated fprintf('More than 32 channels: default stopping weight change 1E-7\n'); end; fprintf('Training will end when wchange < %g or after %d steps.\n', ... nochange,maxsteps); if biasflag, fprintf('Online bias adjustment will be used.\n'); else fprintf('Online bias adjustment will not be used.\n'); end end % %%%%%%%%%%%%%%%%% Remove overall row means of data %%%%%%%%%%%%%%%%%%%%%%% % if verbose, fprintf('Removing mean of each channel ...\n'); end rowmeans = mean(data'); data = data - rowmeans'*ones(1,frames); % subtract row means if verbose, fprintf('Final training data range: %g to %g\n', ... min(min(data)),max(max(data))); end % %%%%%%%%%%%%%%%%%%% Perform PCA reduction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % if strcmp(pcaflag,'on') fprintf('Reducing the data to %d principal dimensions...\n',ncomps); [eigenvectors,eigenvalues,data] = pcsquash(data,ncomps); % make data its projection onto the ncomps-dim principal subspace end % %%%%%%%%%%%%%%%%%%% Perform specgram transformation %%%%%%%%%%%%%%%%%%%%%%% % if exist('Specgramflag') == 1 % [P F T] = SPECGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) % MATLAB Sig Proc Toolbox % Hzwinlen = fix(srate/Hzinc); % CHANGED FROM THIS 12/18/00 -sm Hzfftlen = 2^(ceil(log(Hzwinlen)/log(2))); % make FFT length next higher 2^k Hzoverlap = 0; % use sequential windows % % Get freqs and times from 1st channel analysis % [tmp,freqs,tms] = specgram(data(1,:),Hzfftlen,srate,Hzwinlen,Hzoverlap); fs = find(freqs>=loHz & freqs <= hiHz); if isempty(fs) fprintf('runica(): specified frequency range too narrow!\n'); return end; specdata = reshape(tmp(fs,:),1,length(fs)*size(tmp,2)); specdata = [real(specdata) imag(specdata)]; % fprintf(' size(fs) = %d,%d\n',size(fs,1),size(fs,2)); % fprintf(' size(tmp) = %d,%d\n',size(tmp,1),size(tmp,2)); % % Loop through remaining channels % for ch=2:chans [tmp] = specgram(data(ch,:),Hzwinlen,srate,Hzwinlen,Hzoverlap); tmp = reshape((tmp(fs,:)),1,length(fs)*size(tmp,2)); specdata = [specdata;[real(tmp) imag(tmp)]]; % channels are rows end % % Print specgram confirmation and details % fprintf(... 'Converted data to %d channels by %d=2*%dx%d points spectrogram data.\n',... chans,2*length(fs)*length(tms),length(fs),length(tms)); if length(fs) > 1 fprintf(... ' Low Hz %g, high Hz %g, Hz incr %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),freqs(fs(2))-freqs(fs(1)),Hzwinlen); else fprintf(... ' Low Hz %g, high Hz %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),Hzwinlen); end % % Replace data with specdata % data = specdata; datalength=size(data,2); end % %%%%%%%%%%%%%%%%%%% Perform sphering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % if strcmp(sphering,'on'), %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if verbose, fprintf('Computing the sphering matrix...\n'); end sphere = 2.0*inv(sqrtm(cov(data'))); % find the "sphering" matrix = spher() if ~weights, if verbose, fprintf('Starting weights are the identity matrix ...\n'); end weights = eye(ncomps,chans); % begin with the identity matrix else % weights given on commandline if verbose, fprintf('Using starting weights named on commandline ...\n'); end end if verbose, fprintf('Sphering the data ...\n'); end data = sphere*data; % decorrelate the electrode signals by 'sphereing' them elseif strcmp(sphering,'off') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if ~weights % is starting weights not given if verbose, fprintf('Using the sphering matrix as the starting weight matrix ...\n'); fprintf('Returning the identity matrix in variable "sphere" ...\n'); end sphere = 2.0*inv(sqrtm(cov(data'))); % find the "sphering" matrix = spher() weights = eye(ncomps,chans)*sphere; % begin with the identity matrix sphere = eye(chans); % return the identity matrix else % weights ~= 0 if verbose, fprintf('Using starting weights from commandline ...\n'); fprintf('Returning the identity matrix in variable "sphere" ...\n'); end sphere = eye(chans); % return the identity matrix end elseif strcmp(sphering,'none') sphere = eye(chans); % return the identity matrix if ~weights if verbose, fprintf('Starting weights are the identity matrix ...\n'); fprintf('Returning the identity matrix in variable "sphere" ...\n'); end weights = eye(ncomps,chans); % begin with the identity matrix else % weights ~= 0 if verbose, 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 rand('state',sum(100*clock)); % set the random number generator state to % a position dependent on the system clock %% Compute ICA Weights if biasflag & extended while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% timeperm=randperm(datalength); % shuffle data order at each step for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%% pause(0); if ~isempty(get(0, 'currentfigure')) if strcmp(get(gcf, 'tag'), 'stop') close; error('USER ABORT'); end; end u=weights*data(:,timeperm(t:t+block-1)) + bias*onesrow; y=tanh(u); weights = weights + lrate*(BI-signs*y*u'-u*u')*weights; bias = bias + lrate*sum((-2*y)')'; % for tanh() nonlin. 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 ~wts_blowup % %%%%%%%%%%% Extended-ICA kurtosis estimation %%%%%%%%%%%%%%%%%%%%% %while step < maxsteps 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 & ~wts_blowup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blockno = blockno + 1; if wts_blowup break end end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); 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));; if lrate> MIN_LRATE r = rank(data); if r 2 angledelta=acos((delta*olddelta')/sqrt(change*oldchange)); end if verbose, places = -floor(log10(nochange)); if step > 2 ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, angledelta %4.1f deg, %d subgauss\n',... step, lrate, degconst*angledelta,... places+1,places, (ncomps-sum(diag(signs)))/2); else ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, %d subgauss\n',... step, lrate, places+1,places, (ncomps-sum(diag(signs)))/2); end % step > 2 fprintf('step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ... step, lrate, change, degconst*angledelta); % fprintf(ps,change); % <---- BUG !!!! 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 while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %% Compute ICA Weights if biasflag & ~extended while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% timeperm=randperm(datalength); % shuffle data order 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; u=weights*data(:,timeperm(t:t+block-1)) + bias*onesrow; y=1./(1+exp(-u)); weights = weights + lrate*(BI+(1-2*y)*u')*weights; bias = bias + lrate*sum((1-2*y)')'; % for logistic nonlin. % 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 blockno = blockno + 1; if wts_blowup break end end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 lrate> MIN_LRATE r = rank(data); if r 2 angledelta=acos((delta*olddelta')/sqrt(change*oldchange)); end if verbose, places = -floor(log10(nochange)); if step > 2 ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, angledelta %4.1f deg\n', ... step, lrate, places+1,places, degconst*angledelta); else ps = sprintf('step %d - lrate %5f, wchange %%%d.%df\n',... step, lrate, places+1,places ); end % step > 2 fprintf('step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ... step, lrate, change, degconst*angledelta); % fprintf(ps,change); % <---- BUG !!!! 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 while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %% Compute ICA Weights if ~biasflag & extended while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% timeperm=randperm(datalength); % shuffle data order 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 u=weights*data(:,timeperm(t:t+block-1)); y=tanh(u); % weights = weights + lrate*(BI-signs*y*u'-u*u')*weights; 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 ~wts_blowup % %%%%%%%%%%% Extended-ICA kurtosis estimation %%%%%%%%%%%%%%%%%%%%% %while step < maxsteps 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 ~wts_blowup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blockno = blockno + 1; if wts_blowup break end end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); 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)); if lrate> MIN_LRATE r = rank(data); if r 2 angledelta=acos((delta*olddelta')/sqrt(change*oldchange)); end if verbose, places = -floor(log10(nochange)); if step > 2 ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, angledelta %4.1f deg, %d subgauss\n',... step, lrate, degconst*angledelta,... places+1,places, (ncomps-sum(diag(signs)))/2); else ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, %d subgauss\n',... step, lrate, places+1,places, (ncomps-sum(diag(signs)))/2); end % step > 2 fprintf('step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ... step, lrate, change, degconst*angledelta); % fprintf(ps,change); % <---- BUG !!!! 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 while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Compute ICA Weights if ~biasflag & ~extended while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% timeperm=randperm(datalength); % shuffle data order 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; u=weights*data(:,timeperm(t:t+block-1)); y=1./(1+exp(-u)); % weights = weights + lrate*(BI+(1-2*y)*u')*weights; 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 blockno = blockno + 1; if wts_blowup break end end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 lrate> MIN_LRATE r = rank(data); if r 2 angledelta=acos((delta*olddelta')/sqrt(change*oldchange)); end if verbose, places = -floor(log10(nochange)); if step > 2 ps = sprintf('step %d - lrate %5f, wchange %%%d.%df, angledelta %4.1f deg\n', ... step, lrate, places+1,places, degconst*angledelta); else ps = sprintf('step %d - lrate %5f, wchange %%%d.%df\n',... step, lrate, places+1,places ); end % step > 2 fprintf('step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ... step, lrate, change, degconst*angledelta); % fprintf(ps,change); % <---- BUG !!!! 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 while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %% Finalize Computed Data for Output if ~laststep laststep = step; end; lrates = lrates(1,1:laststep); % truncate lrate history vector % %%%%%%%%%%%%%% Orient components towards max positive activation %%%%%% % if strcmp(pcaflag,'off') sr = sphere * rowmeans'; for r = 1:ncomps data(r,:) = data(r,:)+sr(r); % add back row means end data = weights*data; % make activations from data; -sm 7/05 % add back the row means removed from data before sphering else ser = sphere*eigenvectors(:,1:ncomps)'*rowmeans'; for r = 1:ncomps data(r,:) = data(r,:)+ser(r); % add back row means end data = weights*data; % make activations from sphered and pca'd data; -sm 7/05 % add back the row means removed from data before sphering end % % NOTE: Now 'data' are the component activations = weights*sphere*raw_data % % %%%%%%%%%%%%%% 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 % %%%%%%%%%%%%%%%%%%%% 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 % % compute variances without backprojecting to save time and memory -sm 7/05 % meanvar = sum(winv.^2).*sum((data').^2)/(ncomps-1)^2; % NB: data is now activations -sm 7/05 % %%%%%%%%%%%%%% Sort components by mean variance %%%%%%%%%%%%%%%%%%%%%%%% % [sortvar, windex] = sort(meanvar); windex = windex(ncomps:-1:1); % order large to small meanvar = meanvar(windex); % %%%%%%%%%%%% re-orient max(abs(activations)) to >=0 ('posact') %%%%%%%% % if strcmp(posactflag,'on') % default is now off to save processing and memory fprintf('Making the max(abs(activations)) positive ...\n'); [tmp ix] = max(abs(data')); % = max abs activations signsflipped = 0; for r=1:ncomps if sign(data(r,ix(r))) < 0 if nargout>6 % if activations are to be returned (only) data(r,:) = -1*data(r,:); % flip activations so max(abs()) is >= 0 end winv(:,r) = -1*winv(:,r); % flip component maps signsflipped = 1; end end if signsflipped == 1 weights = pinv(winv)*inv(sphere); % re-invert the component maps end % [data,winvout,weights] = posact(data,weights); % overwrite data with activations % changes signs of activations (now = data) and weights % to make activations (data) net rms-positive % can call this outside of runica() - though it is inefficient! end % %%%%%%%%%%%%%%%%%%%%% Filter data using final weights %%%%%%%%%%%%%%%%%% % if nargout>6, % if activations are to be returned if verbose, fprintf('Permuting the activation wave forms ...\n'); end data = data(windex,:); % data is now activations -sm 7/05 else clear data 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 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % return % %%%%%%%%%%%%%%%%%% DEPRECATED return nonlinearly-transformed data %%%%%%%%%%%%%%%% % if nargout > 7 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