[Master Index] [Index for Toolbox]

rapmusic_gui

(Toolbox/rapmusic_gui.m in BrainStorm 2.0 (Alpha))


Function Synopsis

Results = rapmusic_gui(StudySubject,GUI);

Help Text

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.

Cross-Reference Information

This function calls
This function is called by

Listing of function C:\BrainStorm_2001\Toolbox\rapmusic_gui.m

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

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON