[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

    brainstorm                         C:\BrainStorm_2001\Toolbox\brainstorm.m
    bst_color_scheme                   C:\BrainStorm_2001\Toolbox\bst_color_scheme.m
    bst_imaging_display                C:\BrainStorm_2001\Toolbox\bst_imaging_display.m
    bst_layout                         C:\BrainStorm_2001\Toolbox\bst_layout.m
    bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m
    bst_mriviewer                      C:\BrainStorm_2001\Toolbox\bst_mriviewer.m
    bst_static_taskbar                 C:\BrainStorm_2001\Toolbox\bst_static_taskbar.m
    bst_wavedata_display               C:\BrainStorm_2001\Toolbox\bst_wavedata_display.m
    colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m
    depthgauge                         C:\BrainStorm_2001\Toolbox\depthgauge.m
    extract_channels                   C:\BrainStorm_2001\Toolbox\extract_channels.m
    get_user_directory                 C:\BrainStorm_2001\Toolbox\get_user_directory.m
    good_channel                       C:\BrainStorm_2001\Toolbox\good_channel.m
    load_brainstorm_file               C:\BrainStorm_2001\Toolbox\load_brainstorm_file.m
    norcol                             C:\BrainStorm_2001\Toolbox\norcol.m
    regcheck                           C:\BrainStorm_2001\Toolbox\regcheck.m
    regsubspace                        C:\BrainStorm_2001\Toolbox\regsubspace.m
    results_time_series                C:\BrainStorm_2001\Toolbox\results_time_series.m
    results_visualization              C:\BrainStorm_2001\Toolbox\results_visualization.m
    save_fieldnames                    C:\BrainStorm_2001\Toolbox\save_fieldnames.m
    str_repeater                       C:\BrainStorm_2001\Toolbox\str_repeater.m
    subcorr_gui                        C:\BrainStorm_2001\Toolbox\subcorr_gui.m
    toggleface                         C:\BrainStorm_2001\Toolbox\toggleface.m
    togglelight                        C:\BrainStorm_2001\Toolbox\togglelight.m
    view_surface                       C:\BrainStorm_2001\Toolbox\view_surface.m

This function is called by

    rap_music_parameters               C:\BrainStorm_2001\Toolbox\rap_music_parameters.m
    results_time_series                C:\BrainStorm_2001\Toolbox\results_time_series.m
    results_update                     C:\BrainStorm_2001\Toolbox\results_update.m

### 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