function Results = rapmusic_gui(StudySubject,GUI);
%RAPMUSIC_GUI - Execute RAPMUSIC using GUI inputs, BrainStorm MMII version
% function Results = rapmusic_gui(StudySubject,GUI);
% StudySubject is the conventional structure of strings,
%   StudySubject = find_brainstorm_structure(<studyname>).
% GUI is structure of parameters needed to run this function, 
%   see below.
% All strings represent .mat files that are selectively loaded from disk.
% The string names should represent fully qualified filenames.
% Parameters in GUI are:
%   .DataName, string of actual data file used
%   .Results, optional file to write results
%   .Segment, an index vector giving the column index numbers to use in 
%      the spatiotemporal data matrix.
%   .Order, ordered vector of multipolar orders to process, e.g. [-1 0 1]
%   .OrderHeadModel, one-one with .Order, the referential filename to the
%      corresponding HeadModel to use with each Order
%   .Rank, rank to use, e.g. 10
%   .Corr, correlation to use in processing, e.g. 0.98
%   .ChannelFlag, a vector of Flags in the data file to process, See Data.ChannelFlag
%   .DisplayGraphics, if non-zero, plots information as the algorithm runs
% The following three structure fields are optional for regularization purposes.
%   They may be null or missing to represent unused. If multiple fields are given, 
%   then precedence is given in the order given below.
%     .Condition, condition number to use in truncation, e.g. 100
%     .Energy, fractional energy to use in truncation, e.g. .95
%     .Tikhonov, Tikhonov condition number to use in regularization
%   If all are null, no regularization is performed in the RAP-MUSIC loops.
% Another flag that effects regularization is
%   .Column_norm, if non-zero then columns of gain matrix are normalized in the 
%      original model
%   If GUI.Results is non-null,then writes results to that file.

%<autobegin> ---------------------- 14-Jun-2004 17:12:32 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Inverse Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_color_scheme.m
%   toolbox\bst_imaging_display.m
%   toolbox\bst_layout.m
%   toolbox\bst_message_window.m
%   toolbox\bst_mriviewer.m
%   toolbox\bst_static_taskbar.m
%   toolbox\bst_wavedata_display.m
%   toolbox\colnorm.m
%   toolbox\extract_channels.m
%   toolbox\get_user_directory.m
%   toolbox\good_channel.m
%   toolbox\load_brainstorm_file.m
%   toolbox\norcol.m
%   toolbox\regcheck.m
%   toolbox\regsubspace.m
%   toolbox\results_visualization.m
%   toolbox\save_fieldnames.m
%   toolbox\str_repeater.m
%   toolbox\subcorr_gui.m
%
% Subfunctions in this file, in order of occurrence in file:
%   P = combs(v,m)
%   e = min_fgain_gui(L,Function,Channel,Param,Order,GUI,Uf,Ua,Whitener,imegsens,irefsens,ieegsens);
%
% Application data and their calls in this file:
%   'TileType'
%   'visfig'
%   
%   setappdata(CorticEnvFigure,'TileType','T')
%   setappdata(hTimeSeries,'TileType','T')
%   setappdata(hTimeSeries,'visfig',visfig);
%   
%   visfig = getappdata(hTimeSeries,'visfig');
%
% At Check-in: $Author: Mosher $  $Revision: 44 $  $Date: 6/14/04 3:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 14-Jun-2004 17:12:32 -----------------------

% John C. Mosher, Ph.D.
% Sylvain Baillet, Ph.D.

% ---------------------------------------------------------------------------------
%  6/03/99 JCM,   consolidating the regularization functions
%  6/09/99 JCM,   updating for .DataName and .Results
%  6/24/99 JCM,   updating for visualization of music metric
% 10/18/99 JCM,   adding the 2 synchronous dipole model
% 10/22/99 SB,    updating tessellation file management fo visualization
% JCM 10/26       updating 2 synchronous model handling
% JCM 11/9/99     updating referential filenames, dropping the SearchTess stuff.
% JCM 12/2/99     replacing fmins with fminsearch, fmins now outdated.
% SB Dec-Jan      CTF gradiometer stuff.
% JCM 1/31/00     synchronizing versions
% JCM 2/16/00     Segment is now an index vector, display graphics handled by
%                 results_update
% JCM 2/23/00     some cleanup
% JCM 3/19/00     Added in Baillet's fixes for CTF reference sensors
% JCM 3/20/00     Compatibility with results_update
% SB 4/27/01      Little clean-up and update on bad-channel referencing with 
%                 CTF channel structure
% JCM 2/25/2002   used function handles for subfunction for fminsearch, replaced "disp" with bst_message_window
% JCM 10-May-2002 major changes to headmodel, to conform to MMII approach
%                 i.e. no more "SearchGain", just "Gain", etc.
%                 Have to rethink the different gain matrices as well. The 
%                 magnetic dipole is on hiatus, we have just the current dipole 
%                 and current multipole, orders -1 and 1.
% JCM 13-May-2002 Added GUI.OrderHeadModel, which is a cell-array of
%                 referential strings to the HeadModel corresponding to GUI.Order
% JCM 14-May-2002 Moving HeadModel around to allow that each SourceOrder is called
%                 separately. SB and JCM also updated the HeadModel to have fields 
%                 SourceOrder,HeadModelType. RAP-MUSIC doesn't actually care which 
%                 HeadModelType search grid is called, but the SourceOrder is crucial.
% 24-Jun-2002 JCM Comments changes, upped Pairs to 5000, fixed bug with Data Projector
%                 and the displaying of results. Some of A is now DP (Data Projector)
%                 to distinguish those components of the new model from those of DP.
% 21-Oct-2002 JCM Major: Changed projection operation. Rather than project the signal
%                 subspace away from the found sources, we project away the *data*, then
%                 recompute the signal subspace. In ProjOption = 1 (hardwired just
%                 below these comments), we decrement the signal subspace rank by one for
%                 each recursion. In ProjOption = 2, we keep the rank the same. In
%                 theory, the full rank subspace will become garbage anyhow.
%                 In ProjOption = 0, now deprecated, we keep the old scheme, for 
%                 legacy testing.
% JCM 31-Oct-2002 Bug handling for missing UDP, (orthogonal data projector matrix).
% SB  24-Oct-2003 Changed call to bst_wave_display
% SB  05-Mar-2004 Changed call to bst_wave_display (again)
% JCM 03-May-2004 Cleaning comments. bst_mriviewer is used to display the
%                 final solution. bst_imaging_dispay and
%                 bst_wavedata_display are also called.
% ---------------------------------------------------------------------------------


% ------ HARDWIRED PARAMETERS ---------------
% maximum number of random pairs to search in the initial synchronous
%  grid. Set based on speed performance.
PAIRS = 5000;

% Projection Options, see comments above from 21-Oct-2002 JCM.
% Set ProjOption to 0, 1, or 2
%   0 deprecated, the old style of projecting just the signal subspace:
%   1 recommended, project the data, then find the signal subspace as the new
% dominant subspace vectors, where the rank of the new space has been 
% decremented each time:
%   2 experimental, project the data, then find the signal subspace, but
% keep the rank the same:
ProjOption = 1;

% DIMS is probably deprecated, now that each head model is separate.
DIMS = [3 3 12]; % the number of linear parameters per multipolar oders -1, 0, and 1.

% --------- END HARDWIRED ------------------


% keep the user from getting bored.
bst_message_window('Loading information to run RAP-MUSIC . . .');

Users = get_user_directory; % the default search paths

% load the head model information
if isempty(GUI.OrderHeadModel),
   % User must compute a head model prior to running RAP_MUSIC
   errordlg('Please compute a head model prior to running RAP_MUSIC')
   return
end

% which subject is this, informational
Subject = load_brainstorm_file(StudySubject.Subject);

if(0), % JCM 10-May-2002, won't display tesselated results anymore
   % Subject Information
   if isempty(Subject.Tesselation) % User must assign a Tesselation file to the subject -> a dummy one if necessary
      errordlg('Please assign first a Tesselation file to the Subject file')
      return
   end   
   
   if(GUI.DisplayGraphics)
      
      % Tesselation information, setup up FV for viewing the source solutions.
      %  Every subject must have a tesselation for use, even if surrogate.
      SubjectTess = load_brainstorm_file(Subject.Tesselation);
      %Find cortical tesselations
      % Look for tessellations beguinning with 'cortex' as a keyword -> update documentation
      icrtx =  strmatch('cortex',lower(SubjectTess.Comment)); 
      if(isempty(icrtx)),
         errordlg('Need a tesslation that begins with ''cortex''. Return to Head Modeler');
         return
      end
      
      % If multiple cortical files, Look for the smallest tessellation.
      if(length(icrtx) > 1)
         k = 0;
         for i = icrtx(:)',
            k = k+1;
            siz(k) = size(SubjectTess.Faces{i},1);
         end
         [ignore, imin] = min(siz);
      else
         imin = 1; % only one cortical tesselation
      end
      % assign the cortical tesselation to the structure
      FV = struct('faces',SubjectTess.Faces{icrtx(imin)},...
         'vertices',SubjectTess.Vertices{icrtx(imin)}');
      clear SubjectTess % don't need the rest of the Tesselations
      
      % Create a figure for viewing the source results
      hfdepthgauge = figure; % new figure, used to show depthgauges into cortex
      set(hfdepthgauge,'units','normalized');
      set(hfdepthgauge,'position',[.675 .55 .31 .31]); % upper right corner
      set(hfdepthgauge,'Name',GUI.Results);
      [hfdepthgauge,hpatch] = view_surface(GUI.Results,FV.faces,FV.vertices);
      set(hpatch,'FaceLighting','flat'); % much faster
      uimenu('Label','Toggle Face','callback','toggleface(get(gcbo,''UserData''))','UserData',hpatch);
      uimenu('Label','Toggle Lighting','callback','togglelight(get(gcbo,''UserData''))','UserData',hpatch);
      rotate3d on
      axis tight vis3d
      drawnow
   end
   
end % no more viewing cortical tesselation


% regularization check
GUI = regcheck(GUI); % returns GUI.REG set to the appropriate field

Study = load_brainstorm_file(StudySubject.Study); % and the study information
% don't need SubjectImage, SubjectTess, SubjectTriConn, SubjectVertConn

% now trim the data to just the desired channels
% Use the first head model to establish desired channels
[HeadModel, Channel, GoodChannel, DataFlag] = ...
   extract_channels(GUI.OrderHeadModel{1},StudySubject,GUI);
% GoodChannel represents that subset of channels desired for analysis
% Channel and HeadModel have already been trimmed to this set

% find indices of all MEG and EEG channels
imegsens = good_channel(Channel,[],'MEG');
irefsens = good_channel(Channel,[],'MEG REF'); % may return empty
ieegsens = good_channel(Channel,[],'EEG'); 

% Load the data set
Data = load_brainstorm_file(GUI.DataName);

if ~isempty(Data.Projector),
   if(size(Data.Projector,1) ~= length(GoodChannel)),
      bst_message_window('wrap',{...
            'Projector not same row dimension as GoodChannels.',...
            'Projector removed'});
      Data.Projector = [];
      cd(Users.STUDIES)
      save_fieldnames(Data,GUI.DataName); % save the empty projector
   end
end

% Trim the data to the requested time segments
F = Data.F(:,GUI.Segment);
F = F(GoodChannel,:); % good channels only
Time = Data.Time(GUI.Segment);


% done altering the data and parameters

% did the user provide an existing control. Decompose it
% JCM 24-Jun-2002, switched this A to DP (Data Projector),
% to keep clear what is DP and what is the Independent Topographies
% generated by this fit to the data.
if isfield(Data,'Projector'),
   if(isempty(Data.Projector)),
      DP = [];
   else
      DP = Data.Projector; % handy shorthand
   end
else % user did not give a projector
   DP = [];
end

% Is there a noise covariance for a whitener?

if ~isempty(Data.NoiseCov),
   Whitener = Data.NoiseCov(GoodChannel,:);
   Whitener = pinv(Whitener);
   % now whiten the data and the precomputed gain
   F = Whitener*F;
   HeadModel.Gain{1} = Whitener*double(HeadModel.Gain{1});
else
   Whitener = []; % 
end

% The data could have already been placed into a subspace
%  away from the projector. 
% . Noise: Should already be in a lower dimension
% . Other solution: Projector represents prior solution, so 
%     Data may not be projected.
% Let's ensure.

if(~isempty(DP)), % was there a Data Projector?
   UDP = orth(DP); % make sure that DP is orthogonalized
   F = F - UDP*(UDP'*F); % project data away from projector
else
   UDP = []; % orthog DP is also empty
end

% Now the data exist in a lower dimensional subspace away from the
%  projector. Proceed as normal, building a topographies matrix.
%

% Initialize the Independent Topographies matrix separate
%  from the DP matrix. In the code that follows, Ua effectively represents
%  orth([DP A]), but A starts off as null.

A = []; % A initialized for the loops
Ua = orth([DP A]); % DP and A orthogonalized, works for null case as well
Aorg = []; % the final topographies, no spatial filters applied

bst_message_window('overwrite','Head model loaded, loading data . . .');

% Preprocess the Data

% User has asked for an rank RANK subspace search.
% Form the signal subspace

if(size(F,1) >= size(F,2)), % skinny matrix
   [Uf,Sf,Vf] = svd(F,0); % economy
else % wide matrix
   [Vf,Sf,Uf] = svd(F.',0); % direct transpose
end

% form the reduced rank data
Fr = Uf(:,1:GUI.Rank)*(Sf(1:GUI.Rank,1:GUI.Rank)*Vf(:,1:GUI.Rank)');

% form UPF as the projected signal subspace
UPF = Uf(:,1:GUI.Rank);

% form the structure for the answers in the loops below. 
% pre-allocate some space. No problem if exceeded, just
% slows down some. Will remap these to Results when done.
[best(1:100)] = deal(struct('L',[],'Order',[],...
   'c',[],'O',[],'rank',[]));

% Entering into this loop, we have the predefined topology matrix A and Ua, possibly null
%  if the user gave no initial projector. We will build on this matrix.

nOrder = 1;  % order indexer in GUI.Order. Possible to search in any order, but
%  logically we only anticipate increasing modeling order, such as [-1 0 1];
% New MMII code, load particular headmodel here
%  The first headmodel was already loaded above

Topo_i = 1;  % topology indexer, which source we are processing

bst_message_window('append',{'','',...
      sprintf('************* Run at %s *************',datestr(now))});

while((nOrder <= length(GUI.Order)) & (Topo_i <= GUI.Rank)), % multipolar order maximum
   
   Order = GUI.Order(nOrder); % The source order
   
   % handle synchronous source flags.
   % Orders -1, 0, and 1 are current, magnetic, and 1st-order multipolar
   % Orders 2, 3, 4 are synchronous pairs for -1, 0, and 1.
   % Not elegant, but workable
   % iOrder is used for indexing the search grids and orders, 1:6 one-to-one with -1:4
   switch Order
   case {-1,0,1}
      % we're fine, but set Sync to 1
      Sync = 1; % single sources
      iOrder = Order + 2; % Order -1 is 1, Order 1 is 3, array indexer
   case 2 % synchronous current dipole
      Order = -1; % they are really dipoles
      Sync = 2; % want two synchronous sources
      % FIX, we really need a two dipole grid, not the denser single dipole grid
      iOrder = 1; % Use the synchronous current dipole grid
   case 3 % synchronous magnetic dipole
      Order = 0; % magnetic dipoles
      Sync = 2;
      iOrder = 2; % Use the synchronous magnetic dipole grid
   case 4 % synchronous multipolar
      Order = 1; % 1st-order multipoles
      Sync = 2;
      iOrder = 3; % Use the multipolar dipole grid
   otherwise % -1, 0, 1 orders
      errormsg(sprintf('Don''t know this order: %.0f',Order));
      return
   end
   
   bst_message_window('append',{'',sprintf(...
         'PROCESSING Source %.0f at Order %.0f (%.0f)',Topo_i,Order,Sync),...
         sprintf('with regularization of %s',GUI.REG),''});
   
   % these are not that big, and useful to remap
   if(0), % old calls
      Lgrid = double(HeadModel.GridLoc{iOrder});
      Ggrid = double(HeadModel.Gain{iOrder});
   else % new MMII calls
      Lgrid = double(HeadModel.GridLoc{1});
      Ggrid = double(HeadModel.Gain{1});
   end
   
   
   nLgrid = size(Lgrid,2); % how many are there
   
   % Has user requested column normalization of the grid
   if(GUI.Column_norm),
      % CHEAT, May 2002, to reimplement this, we need to store the "ignored" data
      [ignore,Ggrid] = colnorm(Ggrid);
   end
   
   % first scan the grid looking for a maximum
   % what is the rank of the signal subspace
   best(Topo_i).rank = size(UPF,2);
   
   % project the grid away from the existing model subspace
   if(~isempty(Ua)), % may be empty on the first iteration
      PGgrid = Ggrid - Ua*(Ua'*Ggrid);
   else % ua is empty
      PGgrid = Ggrid; % grid remains unprojected
   end
   
   if(Sync == 1), % searching for single sources
      % form grid of decomposed gain matrices  
      UGrid = zeros(size(Ggrid)); % decomposed grid
      % now must form and possibly truncate each grid point
      Dim = DIMS(iOrder); % how many columns per source
      for i = 1:nLgrid, % for each source location
         ndx = i*Dim+[(1-Dim):0];
         [Utemp,Stemp,Vtemp] = svd(PGgrid(:,ndx),0);
         Utemp = regsubspace(GUI,Utemp,Stemp); % trim to conditioned subspace
         % keep the null space columns to zero
         UGrid(:,ndx(1:size(Utemp,2))) = Utemp;
      end % for each source location 
      
      % We will bypass calling subcorr_gui here, for speed. As of 6/13/99, all
      %  forms of regularization use the unaltered subspaces formed above.
      % UPF is the signal subspace.
      % UPF is already orthogonal and possibly rank reduced.
      C = UPF'*UGrid; % matrix inner product against the grid
      c = zeros(1,nLgrid); % the correlations
      for i = 1:nLgrid, % for each source location
         cndx = i*Dim+[(1-Dim):0]; % which submatrix
         tmp = svd(C(:,cndx));  % decompose it
         c(i) = tmp(1); % maximum correlation between UPF and this source
      end
      [max_c,Di] = max(c); % maximum correlation over the grid
      best(Topo_i).L = Lgrid(:,Di); % location of best source
      best(Topo_i).c = max_c; % corresponding correlation
      best(Topo_i).Order = Order; % model order
      
      bst_message_window('append',sprintf(...
         'Initial solution on Grid: [%s], %.1f%%',...
         str_repeater('%.1f',',',best(Topo_i).L*1000),best(Topo_i).c*100));
      
      bst_message_window('Optimizing in a nonlinear directed search . . .');
      
      % Optimization call
      % optimize the solution found on the grid
      if(0), % outdated fmins call, as of Matlab 5.3
         best(Topo_i).L = fmins(@min_fgain_gui,...
            best(Topo_i).L,[],[],HeadModel.Function,Channel,HeadModel.Param,Order,...
            GUI,UPF,Ua,Whitener);
      else % the newer call
         %See OPTIMSET for details.  FMINSEARCH uses
         %these options: Display, TolX, TolFun, MaxFunEvals, and MaxIter.
         % We're dealing in meters, and 1e-4 is tenths of a millimeter.
         % Let's go for smaller. TolFun is correlation keep to six places
         Options = optimset('Display','off','TolX',1e-6,'TolFun',1e-6);
         % get the display off, and set the tolerance to sub millimeters
         best(Topo_i).L = fminsearch(@min_fgain_gui,...
            best(Topo_i).L,Options,HeadModel.Function,Channel,HeadModel.Param,Order,...
            GUI,UPF,Ua,Whitener,imegsens,irefsens, ieegsens);
      end
      
      
   elseif(Sync > 1), % multiple sources
      
      best(Topo_i).c = 0; %initialize the correlation
      % snipped from nchoosek, adjusted notation
      SingleDim = DIMS(iOrder); % dimension of a single source
      SetDim = SingleDim*Sync; % how many columns per source
      
      % FIX: area to research, better gridding.
      % we have two means of initializing the directed search, either exhaustively
      %  testing every combination of grid points, or seeking random pairs
      % Exhaustive can be too demanding if there are too many grid points.
      % For now, stick to random pairs.
      
      % pick a switch with a comment line
      % switch 'Exhaustive' % exhaustively test all combinations
      switch 'RandomPair' % randomly select pairs for testing
      case 'Exhaustive'
         
         bst_message_window('append',sprintf('Searching through %.0f exhaustive combinations. . .',...
            nchoosek(nLgrid,Sync)));
         
         vIndices = [1:nLgrid]; % index to all dipole numbers, row vec
         for Comb_idx = 1:nLgrid-Sync+1, % indexing through the combinations
            Q = combs(vIndices(Comb_idx+1:nLgrid),Sync-1); % local function call
            next_set = [vIndices(ones(size(Q,1),1),Comb_idx) Q];
            
            % code below virtually identical to RandomPair
            % bst_message_window('append',sprintf('Another %.0f sets',size(next_set,1)));
            % next_set has the next set of possible combinations in the grid
            % form grid of decomposed gain matrices  
            UGrid = zeros(size(Ggrid,1),SetDim*size(next_set,1)); % decomposed grid
            
            for i = 1:size(next_set,1), % for each source combination
               ndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
               sndx = zeros(1,SetDim); % set indexer
               for j = 1:length(next_set(i,:)), % each item in the set
                  sndx([(1-SingleDim):0]+j*SingleDim) = [(1-SingleDim):0] + next_set(i,j)*SingleDim;
               end
               Utemp = PGgrid(:,sndx); % extract the desired columns
               [Utemp,Stemp,Vtemp] = svd(Utemp,0); %decompose
               Utemp = regsubspace(GUI,Utemp,Stemp); % trim to conditioned subspace
               % keep the null space columns to zero
               UGrid(:,ndx(1:size(Utemp,2))) = Utemp;
            end % for each source combination
            
            C = UPF'*UGrid; % inner product
            c = zeros(1,size(next_set,1));
            for i = 1:size(next_set,1),
               cndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
               tmp = svd(C(:,cndx));
               c(i) = tmp(1);
            end
            [max_c,Di] = max(c); % maximum correlation in this set
            if(max_c > best(Topo_i).c), % better combination
               best(Topo_i).c = max_c;
               best(Topo_i).L = Lgrid(:,next_set(Di,:));
               best(Topo_i).Order = GUI.Order(nOrder);
            end
            
         end % next set of combination indices
         
      case 'RandomPair'
         bst_message_window('append',sprintf('Searching through at least %.0f random combinations. . .',PAIRS));
         iPair = 0; % indexer into the pairs
         
         while (iPair < PAIRS);
            % create random unique pairs in the grid
            [ignore,next_set] = sort(randn(1,nLgrid));
            nLgrid_trim = nLgrid - rem(nLgrid,Sync); % nearest multiple
            next_set = reshape(next_set(1:nLgrid_trim),nLgrid_trim/Sync,Sync);
            iPair = iPair + size(next_set,1); % next set of pairs.
            % each row of next_set is uniquely paired, no twins.
            
            % code below virtually identical to Exhaustive
            % bst_message_window('append',sprintf('Another %.0f sets',size(next_set,1)));
            % next_set has the next set of possible combinations in the grid
            % form grid of decomposed gain matrices  
            UGrid = zeros(size(Ggrid,1),SetDim*size(next_set,1)); % decomposed grid
            
            for i = 1:size(next_set,1), % for each source combination
               ndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
               sndx = zeros(1,SetDim); % set indexer
               for j = 1:length(next_set(i,:)), % each item in the set
                  sndx([(1-SingleDim):0]+j*SingleDim) = [(1-SingleDim):0] + next_set(i,j)*SingleDim;
               end
               Utemp = PGgrid(:,sndx); % extract the desired columns
               [Utemp,Stemp,Vtemp] = svd(Utemp,0); %decompose
               Utemp = regsubspace(GUI,Utemp,Stemp); % trim to conditioned subspace
               % keep the null space columns to zero
               UGrid(:,ndx(1:size(Utemp,2))) = Utemp;
            end % for each source combination
            
            C = UPF'*UGrid; % inner product
            c = zeros(1,size(next_set,1));
            for i = 1:size(next_set,1),
               cndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
               tmp = svd(C(:,cndx));
               c(i) = tmp(1);
            end
            [max_c,Di] = max(c); % maximum correlation in this set
            if(max_c > best(Topo_i).c), % better combination
               best(Topo_i).c = max_c;
               best(Topo_i).L = Lgrid(:,next_set(Di,:));
               best(Topo_i).Order = GUI.Order(nOrder);
            end
            
         end % next set of combination indices
         
      end % case switch of random vs. exhaustive searching
      
      bst_message_window('overwrite',sprintf(...
         'Initial solution on Grid: [%s], %.1f%%',...
         str_repeater('%.1f',',',best(Topo_i).L*1000),best(Topo_i).c*100));
      
      bst_message_window('Optimizing in a nonlinear directed search . . .');
      
      % Optimization call
      % optimize the solution found on the grid
      if(0), % outdated fmins call, as of Matlab 5.3
         best(Topo_i).L = fmins(@min_fgain_gui,...
            best(Topo_i).L,[],[],HeadModel.Function,Channel,HeadModel.Param,Order,...
            GUI,UPF,Ua,Whitener);
      else % the newer call
         %See OPTIMSET for details.  FMINSEARCH uses
         %these options: Display, TolX, TolFun, MaxFunEvals, and MaxIter.
         % We're dealing in meters, and 1e-4 is tenths of a millimeter.
         % Let's go for smaller.
         Options = optimset('Display','off','TolX',1e-6,'TolFun',1e-6);
         % get the display off, and set the tolerance to sub millimeters
         best(Topo_i).L = fminsearch(@min_fgain_gui,...
            best(Topo_i).L,Options,HeadModel.Function,Channel,HeadModel.Param,Order,...
            GUI,UPF,Ua,Whitener,imegsens,irefsens,ieegsens);
      end
      
   end %  handling the synchronization cases
   
   % we have completed the localization, ready to see if it's any good
   % we have an improved localization
   new_corr = -min_fgain_gui(best(Topo_i).L,HeadModel.Function,...
      Channel,HeadModel.Param,Order,GUI,UPF,Ua,Whitener,imegsens,irefsens, ieegsens);
   bst_message_window('overwrite',sprintf(...
      'Improved solution: [%s], %.1f%%',...
      str_repeater('%.1f',',',best(Topo_i).L*1000),new_corr*100));
   
   if(new_corr < best(Topo_i).c-100*eps),
      % somehown fmins came back more poorly than the grid. probably an error somewhere
      bst_message_window('append',' ')
      bst_message_window('append','Why is the new corr less than the gridded?')
      bst_message_window('append',sprintf('New corr: %.2f%%, gridded: %.2f%%',...
         100*new_corr, 100*best(Topo_i).c));
      bst_message_window('append','Pausing at keyboard for investigation, ''dbcont'' to continue.')
      keyboard
   end
   % assign the new correlation to the best of the results
   best(Topo_i).c = new_corr;
   
   if(best(Topo_i).c >= GUI.Corr) % we have an acceptable solution!
      
      % CHEAT: use the first name only
      % recompute to get the corresponding orientation
      G1 = feval(HeadModel.Function{1},best(Topo_i).L,Channel,HeadModel.Param,Order,imegsens,irefsens, ieegsens); % single source
      G1org = G1; % original topography, no whitener, no projector
      
      if(GUI.Column_norm), % column normalization
         [ignore,G1] = colnorm(G1);
      end
      
      if ~isempty(Whitener),
         G1 = Whitener*G1;
      end
      
      if ~isempty(Ua), % if there already existing projectors
         PpG1 = G1 - Ua*(Ua'*G1); % project away from existing subspace
      else
         PpG1 = G1;
      end  
      
      % Note that here we get the orientation in the reduced subspace.
      % get the associated orientation, regularized if need be.
      [Co,X,Y] = subcorr_gui(PpG1,UPF,GUI);
      best(Topo_i).O = X(:,1)/norm(X(:,1)); % best unit norm orientation
      
      bst_message_window('append',{sprintf(...
            'Orientation [%s]',...
            str_repeater('%.2f',',',best(Topo_i).O)),...
            sprintf('Rank %.0f Source Order %.0f',...
            best(Topo_i).rank,Order)});
      
      
      % need to form the current topography, append to existing
      %  and project the signal subspace away from it
      A = [A G1*best(Topo_i).O];
      Aorg = [Aorg G1org*best(Topo_i).O]; % unfiltered for visualization
      
      Ua = orth([DP A]); % orthogonalize the existing space with the projector
      
      
      % re-form UPF as the projected signal subspace      
      % 21 October 2002 Projection Options
      % See Comments at top, plus hardwired section at top
      
      switch ProjOption
      case 0 
         % legacy, now a deprecated method
         UPF = Uf(:,1:GUI.Rank); % get the original signal subspace
         UPF = UPF - Ua*(Ua'*UPF); % project away existing subspace
         % orthogonalize the reduced subspace
         [UPF,SPF,VPF] = svd(UPF,0);
         
         
         % %%%%%%%%%% NOT A CHEAT BUT LOOK AT THIS %%%%
         % we've had a lot of discussions about the next line.  Do we trim off the weakest space
         %  in the next recursion. Earlier concerns were about regularization
         % JCM thoughts 2/23/00, why not run regsubspace?
         switch 'reduce'
         case 'reduce'
            UPF = UPF(:,1:(end-Topo_i)); % remove the weakest vector spaces
         case 'retain'
            % do nothing, don't truncate
         end
         
      case {1,2}
         % recommended 21 Oct 2002
         % User has asked for an rank RANK subspace search.
         % Form the signal subspace
         
         ReducedF = F - Ua*(Ua'*F); % The data in the new subspace
         
         if(size(ReducedF,1) >= size(ReducedF,2)), % skinny matrix
            [UPF,ignore,ignore] = svd(ReducedF,0); % economy
         else % wide matrix
            [ignore,ignore,UPF] = svd(ReducedF.',0); % direct transpose
         end
         
         % form UPF from the signal subspace
         UPF = UPF(:,1:GUI.Rank);
         
         if(ProjOption == 1), % remove the weakest vector spaces from the subspace
            UPF = UPF(:,1:(end-Topo_i)); % remove the weakest vector spaces
         end
         
         
      otherwise
         error('Unknown Projection Method hardwired into code.');
         bst_message_window('Contact brainstorm@sipi.usc.edu about error');
      end
      
      
      
      if(GUI.DisplayGraphics),
         
         % setup the GUIs for the display of time series
         global data % Get display parameters from Viewer
         if Topo_i == 1
            hTimeSeries = figure;
            setappdata(hTimeSeries,'TileType','T')
            bst_color_scheme(hTimeSeries);
            bst_layout('align',hTimeSeries,2,2,4)
            [ignore,ResultsFileName] = fileparts(GUI.Results);
            set(hTimeSeries,'Name',ResultsFileName);
            data.Results.Display.Handles.Figures.TimeSeries = hTimeSeries;
         end
         
         % Graphics option, could be a lot better:
         % let user see a solution
         % need to regularize
         if(~isempty(UDP)),
            Ap = A - UDP*(UDP'*A); % project model away
         else
            Ap = A; % no projector
         end
         % Data already projected away
         S = [pinv(Ap)*Fr]'; % time series associated with the reduced rank data
         Frsyn = Ap*S'; % the data synthesized in the whitened, projected space
         Frem = Fr-Frsyn; % remainder
         
         SourceLoc = cell(1,Topo_i);
         SourceOrientation = cell(1,Topo_i);
         for i = 1:Topo_i,
            SourceLoc{i} = best(i).L; % map each solution into a cell.
            SourceOrientation{i} = best(i).O;
         end
         
         % temporary results structure
         if(0) % no longer displaying this way
            tResults = struct('IndepTopo',Ap,...
               'TimeSeries',S,'Time',GUI.Segment,'SourceOrder',[best(1:Topo_i).Order],...
               'Comment',sprintf('RAP-MUSIC, rank %.0f, corr %.3f',GUI.Rank,GUI.Corr),...
               'Subject',Subject,'Study',Study,'SourceCorrelation',[best(1:Topo_i).c],...
               'SourceLoc',[],'SourceOrientation',[],'GUI',GUI,'Function',[]);
            tResults.SourceLoc = SourceLoc; % would otherwise DEAL
            tResults.SourceOrientation = SourceOrientation;
            tResults.Function = mfilename; % name of this calling routine
            % send the reduced rank time series. User threw away the rest            
            results_time_series('update',hTimeSeries,struct('F',Fr,'Time',Data.Time),tResults);
         end
         
         
         % -> Display source time series
         
         D = num2cell(S',2);
         for src = 1:length(D)
            SrcName{src} = sprintf('Source %d',src);
         end
         GoodTypes = {Channel(1).Type};
         chans = good_channel(Channel,[],GoodTypes{1});
         
         DisplayStruct = data.Measures.Display;
         %bst_layout('align',data.Results.Display.Figure(1),2,2,4);
         
         OPTIONS = struct(...
            'ModalityLabel',{GoodTypes},...
            'WaveLabel',{SrcName},...
            'Channel', data.Measures.Channel(data.Measures.Display.SelectedChannels{DisplayStruct.Modality}),...
            'IndepTopo',Aorg,...
            'Time', Time, ...
            'MoveCursor','on'...
            );
        OPTIONS.SelectedChannels = data.Measures.Display.SelectedChannels(data.Measures.Display.Modality);
         OPTIONS.DataType = 'source';
         
         figure(hTimeSeries), clf
         figs = data.Results.Display.Handles.Figures;
         
         [data.Results.Display.Handles,OPTIONS] = ...
            bst_wavedata_display(D,OPTIONS);
         data.Results.Display.Handles.Figures = figs; clear figs
         % <- DONE
         
          
         % Now display original and modeled data + residuals
         DisplayStruct = data.Results.Display;
         dOPTIONS = struct(...
            'ModalityLabel',{GoodTypes},...
            'WaveLabel',{{'Original data','Modeled',sprintf('Residuals: %3.2f%%',...
                  100*mean(norcol(Frem))/mean(norcol(Fr)))}},...
            'Channel', data.Measures.Channel(data.Measures.Display.SelectedChannels{data.Measures.Display.Modality}),...
            'Time', Time, ...
            'MoveCursor','on'...
            );
         dOPTIONS.SelectedChannels =  cell(1,3);
         dOPTIONS.SelectedChannels{data.Measures.Display.Modality} = chans;
         % -> Fill out parameters before call to specific display subfunction
         DefaultLabels = {'MEG','EEG','OTHER'}; % Default modality labels
         
         figure(data.Measures.Display.Handles.Figures.TimeSeries), clf
         figs = data.Measures.Display.Handles.Figures;
         % was Time(GUI.Segment)`
         dOPTIONS.DataType = 'surface';
         [data.Measures.Display.Handles,dOPTIONS]  = bst_wavedata_display({Fr,Frsyn,Frem},dOPTIONS);         
         
         data.Measures.Display.Flags = struct('SensorMarkers',OPTIONS.SensorMarkers,...
            'SensorLabels',OPTIONS.SensorLabels, 'ShowContours', OPTIONS.ShowContours,...
             'AbsoluteColormap',OPTIONS.AbsoluteColormap, 'Colorbar',OPTIONS.Colorbar);
         data.Measures.Display.Handles.Figures =  figs;
         
         drawnow
         
         if(0), % 10-May-2002 no more spatial display
            if(hfdepthgauge), %is there a figure to draw to?
               depth_point = best(Topo_i).L;
               figure(hfdepthgauge)
               hlinedepth = zeros(size(depth_point,2),1); % for each source location
               for ihline = 1:length(hlinedepth),
                  [ignore,ignore,ignore,ignore,hlinedepth(ihline)] = ...
                     depthgauge(depth_point(:,ihline),...
                     depth_point(:,ihline)/norm(depth_point(:,ihline)),...
                     100/1000,5/1000,1000,gca);
               end
               
               % set the color of the depthgauge to that of the default line colors in the axes
               tempcolor = get(gca,'ColorOrder');
               mod_i = mod(Topo_i-1,size(tempcolor,1))+1; % modulo the length
               set(hlinedepth,'Color',tempcolor(mod_i,:),'MarkerFaceColor',tempcolor(mod_i,:))
            end
            
            drawnow
            pause(.25) % see if this helps force a screen image to appear
         end % no more spatial display
         
         % Here is the new spatial display (Nov-11 2003 - SB)
         % JCM mods 13 Nov to account for missing tesselations
         TessFile = bst_static_taskbar('GET','TESS');
         existsTessFile = 0; % assume no
         if ~isempty(TessFile) & exist(TessFile,'file')
            existsTessFile = 1; % does exist
         end           
         
         if existsTessFile, 
            if Topo_i == 1 % Fisrt source found %data.Results.Display.CorticalEnvelope
               CorticEnvFigure = figure;
               set(CorticEnvFigure,'visible','off' ); drawnow
               setappdata(CorticEnvFigure,'TileType','T')
               bst_layout('align',CorticEnvFigure,1,2,1);
               bst_layout('align',data.Results.Display.Handles.Figures.TimeSeries,2,2,4);
               bst_layout('align',data.Measures.Display.Handles.Figures.TimeSeries,2,2,2)
               
               ImageGrid.GridLoc = {TessFile};
               bst_message_window('wrap',{...
                     ' ',...
                     sprintf('Loading %s tessellation file',TessFile)});
               
               Users = get_user_directory;
               load(fullfile(Users.SUBJECTS,TessFile),'Comment')
               
               ImageGrid.iGrid = 1;
               if length(Comment)>1 % Several surfaces in this tessellation file
                  [ImageGrid.iGrid,OK] = listdlg('PromptString','Select an envelope from list:',...
                     'SelectionMode','single',...
                     'ListString',Comment);
                  if OK == 0
                     return
                  end
               end
               envOPTIONS.Dipoles.SourceLoc = {best(1).L};
               envOPTIONS.Dipoles.Orientation = {best(1).O};
               envOPTIONS.Dipoles.TimeSeries = [D{:}]';
               envOPTIONS.Time = 1;
               envOPTIONS.currDensityTransparency = 1;
               envOPTIONS.Dipoles.Color = data.Results.Display.Handles.LineColor;
               
               [PlotHandles,envOPTIONS] = bst_imaging_display(...
                  CorticEnvFigure,ImageGrid,...
                  0,envOPTIONS);
               envOPTIONS.hPatch = PlotHandles.hPatch; 
               
            else % Update display with next sources
               
               envOPTIONS.Dipoles.SourceLoc = {best(1:Topo_i).L};
               envOPTIONS.Dipoles.Orientation = {best(1:Topo_i).O};
               envOPTIONS.Dipoles.TimeSeries = reshape([D{:}],length(D{1}),Topo_i);
               envOPTIONS.Dipoles.Color = data.Results.Display.Handles.LineColor;
               
               [PlotHandles,envOPTIONS] = bst_imaging_display(...
                  CorticEnvFigure,ImageGrid,...
                  0,envOPTIONS);
               %envOPTIONS.hPatch = envPlotHandles.hPatch; 
               
            end % if first topo or subsequent topographies
         else
            
            CorticEnvFigure = []; % no tesselation figure
            
         end % if a tesselation file exists
         
      end % if user wanted graphics
      
      % Note that the search grid will always be projected at the beginning of the loop
      
      Topo_i = Topo_i + 1; % increment the counter to the next source topography
      
   else % correlation too low, increment order
      
      bst_message_window('append',{'',...
            sprintf('Correlation too low at Order %.0f.  Next Order',...
            GUI.Order(nOrder))})
      nOrder = nOrder + 1; % order indexer
      if(nOrder <= length(GUI.Order)), % more orders left to search
         % then load the next head model associated with the new order
         % The Channel, GoodChannel, DataFlag don't change, just the HeadModel call
         HeadModel = extract_channels(...
            GUI.OrderHeadModel{nOrder},StudySubject,GUI); % next HeadModel
      end
      
   end % checking cross correlation value
   
end % while valid orders

if GUI.DisplayGraphics & existsTessFile,
   data.Results.Display.Figure(2) = CorticEnvFigure;
   data.Results.Display.OPTIONS = envOPTIONS;
end


% done with processing, now mop up.

best(Topo_i:end) = []; % remove last failed source and pre-allocated zeros

if(isempty(A)),% nothing, honey
   S = [];
   msgbox('No Sources Found','Notice','modal');
   return
end


if(~isempty(UDP)),
   Ap = A - UDP*(UDP'*A); % project model away
else
   Ap = A; % no projector
end

% The data have already been projected away from the projector
S = [pinv(Ap)*Fr]'; % reduced rank time series, Hermitian transpose

% Now synthesize the data without the projector applied
Fsynth = Aorg * S'; % the synthesized data, in the unwhitened unprojected space


% now map to results

Comment = sprintf('RAP-MUSIC, rank %.0f, corr %.3f',GUI.Rank,GUI.Corr);
Date = datestr(datenum(now),0);
Subject = Subject;
Study = Study;
Time = [GUI.Segment];
ChannelFlag = GUI.ChannelFlag;
NoiseCov = [];
SourceCov = [];
Projector = Data.Projector;
SourceOrder = [best.Order];
SourceCorrelation = [best.c]; % the source correlations
SourceLoc = cell(1,length(best));
SourceOrientation = cell(1,length(best));
Aall = zeros(size(Ap));
for i = 1:length(best),
   SourceLoc{i} = best(i).L; % map each solution into a cell.
   SourceOrientation{i} = best(i).O;
end
ModelGain = [];
IndepTopo = Aorg; % Fixed orientation gain matrix
TimeSeries = S; % time series normalized by gain columns
PatchNdx = [];
PatchAmp = [];
ImageGridAmp = [];
ImageGridTime =[];
Function = mfilename; % name of this calling routine

Results = struct('Comment',Comment,'Function',mfilename,'StudySubject',StudySubject,...
   'DataFlag',DataFlag,'GUI',GUI,'Date',Date,'Subject',Subject,'Study',Study,'Time',Time,...
   'ChannelFlag',ChannelFlag,'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
   'Projector',Projector,'SourceOrder',SourceOrder,'SourceCorrelation',SourceCorrelation,...
   'SourceOrientation',[],'ModelGain',ModelGain,...
   'IndepTopo',IndepTopo,'TimeSeries',TimeSeries,'PatchNdx',PatchNdx,...
   'PatchAmp',PatchAmp,'ImageGridAmp',ImageGridAmp,'ImageGridTime',...
   ImageGridTime,'Fsynth',Fsynth);
Results.SourceOrientation = SourceOrientation; % struct command would deal out this variable
Results.SourceLoc = SourceLoc; % ditto

if(isfield(GUI,'Results')), %user provided such a string
   if(~isempty(GUI.Results)), % and it's not empty    
      save(fullfile(Users.STUDIES,GUI.Results),...
         'Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
         'DataFlag','ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder','SourceCorrelation',...
         'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
         'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
      
      bst_message_window('append',{'',sprintf('Wrote results to %s',GUI.Results)}); % to the command window
   end
end

bst_message_window('append',{'',sprintf('RAP-MUSIC ALGORITHM DONE')}); % console screen

if GUI.DisplayGraphics,

   figure(hTimeSeries); % bring to the front
   drawnow
   visfig = getappdata(hTimeSeries,'visfig');
   if isempty(visfig) | ~ishandle(visfig), % not created or gone
      visfig = results_visualization;
      setappdata(hTimeSeries,'visfig',visfig); % save it to prevent duplication
   end
   vishandles = guidata(visfig);
   results_visualization('loadparametric_Callback',...
      vishandles.loadparametric, [], vishandles, Results.GUI.Results);
   
end

if GUI.FinalDisplayGraphics,
   
    bst_mriviewer('LoadSourceMap_Callback','','','',Results.GUI.Results)
    
    %     visfig = results_visualization;
    %     vishandles = guidata(visfig);
    %     results_visualization('loadparametric_Callback',...
    %         vishandles.loadparametric, [], vishandles, Results.GUI.Results);
   
end




%---------------------------------------------------------------------
% snipped from nchoosek function in Matlab
function P = combs(v,m)
%COMBS  All possible combinations.
%   COMBS(1:N,M) or COMBS(V,M) where V is a row vector of length N,
%   creates a matrix with N!/((N-M)! M!) rows and M columns containing
%   all possible combinations of N elements taken M at a time.
%
%   This function is only practical for situations where M is less
%   than about 15.

if nargin~=2, error('Requires 2 input arguments.'); end

v = v(:).'; % Make sure v is a row vector.
n = length(v);
if n == m
   P = v;
elseif m == 1
   P = v.';
else
   P = [];
   if m < n & m > 1
      for k = 1:n-m+1
         Q = combs(v(k+1:n),m-1); % recursively calling
         P = [P; [v(ones(size(Q,1),1),k) Q]];
      end
   end
end





% --------------- THE FMINS FUNCTION, USE WITH FUNCTION HANDLE --------------------
function e = min_fgain_gui(L,Function,Channel,Param,Order,GUI,Uf,Ua,Whitener,imegsens,irefsens,ieegsens);
%MIN_FGAIN_GUI minimize the error using the a specified gain model
% LIMITED: ONLY LOOKS AT THE FIRST FUNCTION CALL AND APPLIES EQUALL TO ALL CHANNELS!
% function e = min_fgain_gui(L,Function,Channel,Param,Order,GUI,Uf,Ua);
% call as fmins('min_fgain_gui',L,[],[],Function,Channel,Param,Order,GUI,Uf,Ua);
% Uf and Ua must already be orthogonal
% Uf is the subspace with which to correlate
% Ua is an existing subspace to project away from
% Function,Channel,Param,Order
% Function is the name of the gain matrix call, e.g. "os_meg"

% John C. Mosher, Ph.D.
% See Copyright.m file for information
% $Date: 6/14/04 3:38p $ $Revision: 44 $

% generate the gain matrix

global GoodChannelMEG GoodChannelEEG 

%%%%% CHEAT  %%%%%%%%
% Fix is to look for all calls to the same function and call in clusters. 
if(0) % too expensive too call one at a time.
   DIMS = [3 3 12];
   Dims= DIMS(Order+2);
   G = zeros(length(Channel),Dims);
   for i =1:length(Channel),
      G(i,:) = feval(Function{i},L,Channel(i),Param(i),Order);
   end
else
   % CHEAT: CALL ONLY THE FIRST FUNCTION
   global ChannelFlag
   ChannelFlag = GUI.ChannelFlag; % Need info on bad channels for CTF gradient computation.
   G= feval(Function{1},L,Channel,Param,Order,imegsens,irefsens,ieegsens);
end

if ~isempty(Whitener),
   G = Whitener*G; % data were already whitened
end


if(0), % disabled by JCM is irrelevant
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % The following is relevant for BEM models only
   if ~isempty(findstr(lower(Function{1}),'meg'))
      if length(GoodChannelMEG) < size(G,1)  % There are bad MEG Channels
         G = G(GoodChannelMEG-GoodChannelMEG(1)+1,:);
      end
   else
      if ~isempty(findstr(lower(Function{1}),'eeg'))
         if length(GoodChannelEEG) < size(G,1)  % There are bad EEG Channels
            G = G(GoodChannelEEG-GoodChannelEEG(1)+1,:);
         end
      end   
   end
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end


if(GUI.Column_norm), % scale the model
   [ignore,G] = colnorm(G);
end

% if the user provided previous recursions
%  project away from it
if(exist('Ua','var')),
   if(~isempty(Ua))
      G = G - Ua*(Ua'*G);
   end
end

% since this function is repeatedly called in directed minimization, reluctant
%  to call subcorr_gui here. More efficient here because we know that Uf is
%  already orthogonal.

% Orthogonalize
[Ug,Sg,Vg] = svd(G,0);

% now set Ug to the reduced rank called for by the regularization
Ug = regsubspace(GUI,Ug,Sg);

% subspace correlations
c = svd(Ug'*Uf);

if(isempty(c)), % nothing valid
   e = 0;
else
   e = -c(1); % negative maximum for minimizing
end
