[Master Index] [Index for Toolbox]

# least_squares_fit

## (Toolbox/least_squares_fit.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    Results = least_squares_fit(StudySubject,GUI);

### Help Text

    LEAST_SQUARES_FIT - multiple dipole fit using fminsearch
     function Results = least_squares_fit(StudySubject,GUI);
     StudySubject is the conventional structure of strings.
     GUI is structure of parameters needed to run this function
     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.
     .Temporal, string of '(S)tationary' or '(M)oving'.
     .Order, cell array of source types -1, 0, 1
     .StartingPoints, cell array, one to one with Order, each cell comprising
        a 3 x 1 vector of starting points. A column with NaNs indicates
        a random grid search.
     .Rank, rank to use, e.g. 10, for the data space
     .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.
     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

    colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m
    depthgauge                         C:\BrainStorm_2001\Toolbox\depthgauge.m
    fronorm                            C:\BrainStorm_2001\Toolbox\fronorm.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
    regcheck                           C:\BrainStorm_2001\Toolbox\regcheck.m
    regsubspace                        C:\BrainStorm_2001\Toolbox\regsubspace.m
    results_update                     C:\BrainStorm_2001\Toolbox\results_update.m
    str_repeater                       C:\BrainStorm_2001\Toolbox\str_repeater.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

    least_squares_parameters_cb        C:\BrainStorm_2001\Toolbox\least_squares_parameters_cb.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\least_squares_fit.m

    function Results = least_squares_fit(StudySubject,GUI);
    %LEAST_SQUARES_FIT - multiple dipole fit using fminsearch
    % function Results = least_squares_fit(StudySubject,GUI);
    % StudySubject is the conventional structure of strings.
    % GUI is structure of parameters needed to run this function
    % 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.
    % .Temporal, string of '(S)tationary' or '(M)oving'.
    % .Order, cell array of source types -1, 0, 1
    % .StartingPoints, cell array, one to one with Order, each cell comprising
    %    a 3 x 1 vector of starting points. A column with NaNs indicates
    %    a random grid search.
    % .Rank, rank to use, e.g. 10, for the data space
    % .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.
    % 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> ---------------------- 26-May-2004 11:30:48 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Inverse Modeling
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\colnorm.m
    %   toolbox\depthgauge.m
    %   toolbox\fronorm.m
    %   toolbox\get_user_directory.m
    %   toolbox\good_channel.m
    %   toolbox\load_brainstorm_file.m
    %   toolbox\regcheck.m
    %   toolbox\regsubspace.m
    %   toolbox\results_update.m
    %   toolbox\str_repeater.m
    %   toolbox\view_surface.m
    %
    % Subfunctions in this file, in order of occurrence in file:
    %   e = min_least_squares_fit(L,Function,Channel,Param,Order,GUI,UfSf,Ua);
    %
    % At Check-in: $Author: Mosher $  $Revision: 15 $  $Date: 5/26/04 9:59a $
    %
    % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-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> ------------------------ 26-May-2004 11:30:48 -----------------------

    % Author: John C. Mosher, Ph.D.

    % ----------------------------- Script History ---------------------------------
    %  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
    % Sylvain Baillet, Ph.D.
    % 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.
    % Silvin 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 Created least-squares from the rapmusic_gui function
    % JCM 3/19/00 Heavily modified to achieve basic least-squares
    % JCM 3/20/00 Works with current dipoles only
    % JCM 2/25/2002 used function handle to fminsearch for subfunction
    % JCM 19-May-2004 comments cleaning
    % ----------------------------- Script History ---------------------------------

    % Hardwired parameters:
    % maximum number of random pairs to search in the initial synchronous
    %  grid. Set based on speed performance.
    PAIRS = 4000;

    DIMS = [3 3 12]; % the number of sources per multipolar oders -1, 0, and 1.

    NumSources = length(GUI.Order);
    if(any(diff([GUI.Order{:}]))),
       errordlg('Sorry, can not handle mixed models for now.')
       % the issue is the design of search grids for mixed types
       return
    end

    % keep the user from getting bored.
    hwait = waitbar(.25,'Loading information to run LEAST-SQUARES');
    set(hwait,'units','normalized');
    Vhwait = get(hwait,'position'); % get position in normalized units
    set(hwait,'position',[.015 .85 Vhwait(3:4)]); % set to upper corner
    drawnow

    Users = get_user_directory; % the default search paths

    % load the head model information. Do a partial load to bypass ImageGain loads.

    if isempty(StudySubject.HeadModel),
       % User must compute a head model prior to running RAP_MUSIC
       errordlg('Please compute a head model prior to running RAP_MUSIC')
       close(hwait);
       return
    end

    HeadModel = load(fullfile(Users.STUDIES,StudySubject.HeadModel),'Function','Param',...
       'SearchGridLoc','SearchGain'); % expensive load, if SearchGain is not stored as binary

    % Subject Information
    Subject = load_brainstorm_file(StudySubject.Subject);
    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')
       close(hwait);
       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');
          close(hwait);
          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

    waitbar(0.5,hwait);
    drawnow

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

    Study = load_brainstorm_file(StudySubject.Study); % and the study information

    % Load the data set
    Data = load_brainstorm_file(GUI.DataName);
    % Trim the data to the requested time segments
    F = Data.F(:,GUI.Segment);
    Time = Data.Time(GUI.Segment);

    if(0), % not using symmetric axes yet (for stereo locations)
       % what is the axis of symmetry
       switch deblank(lower(Data.Device))
       case 'ctf_axial'
          % left/right symmetry is in the y axis
          AxisOfSymmetry = [1 -1 1]';
       case 'neuromag_planar_122'
          AxisOfSymmetry = [-1 1 1]';
       otherwise
          errordlg(sprintf(...
             'Unknown Data Device: %s\nDefaulting to x-axis symmetry',Data.Device))
          AxisOfSymmetry = [-1 1 1]';
       end
       % Multiplying a location by the AxisOfSymmetry yields the symmetric point
    end

    % load up the information needed for the gain matrix
    Channel = load_brainstorm_file(StudySubject.Channel);
    Channel = Channel.Channel; % reassign structure

    % now alter the data according to the bad channels
    % Note that Data.ChannelFlag has the original channel flag information on
    %  Bad and not Bad. GUI.ChannelFlag set to one says analyze that channel
    GoodChannel = find(GUI.ChannelFlag == 1); %only those selected for analysis

    % Are there special reference channel considerations?
    %  See Channel.mat structure description in the ParameterDescriptions document.
    % Assign the reference information of Ctf channels to the first MEG channel
    if(0),
    %if(isfield(Channel,'RefChannel')),
       % there might possibly be CTF reference channel information
       %good channels indexer to ALL MEG channels, good and bad
       MEGChannel = good_channel(Channel,ones(size(GUI.ChannelFlag)),'MEG');

       % Critical case where the reference information should be stored
       %  in the first good MEG channel
       Channel(GoodChannel(1)).RefChannel = Channel(MEGChannel(1)).RefChannel;

       tmp = Channel(MEGChannel(1)).grad3;
       Channel(GoodChannel(1)).grad3 = tmp(:,GoodChannel-min(MEGChannel)+1);
       clear tmp
    end

    % Determine which modality has been chosen for source localization here
    GoodChannelMEG = good_channel(Channel,GUI.ChannelFlag,'MEG'); %good channels indexer
    GoodChannelEEG = good_channel(Channel,GUI.ChannelFlag,'EEG');

    if isempty(GoodChannelMEG)
       %EEG
       GoodChannel = GoodChannelEEG;
       DataFlag = 'EEG';

       % Reference Channel
       ref = strmatch('AVERAGE REF',{Channel.Comment});
       if length(ref) > 2 % Average reference
          ref = 0;
          HeadModel.Param = HeadModel.Param(GoodChannel);
          HeadModel.Function = HeadModel.Function(GoodChannel);
          Channel = Channel(GoodChannel);
       else
          ref = strmatch('EEG REF',{Channel.Comment});
          HeadModel.Param = HeadModel.Param(sort([GoodChannel,ref]));
          HeadModel.Function = HeadModel.Function(sort([GoodChannel,ref]));
          Channel = Channel(sort([GoodChannel,ref]));
       end

    else
       %MEG
       GoodChannel = GoodChannelMEG;
       if(isfield(Channel,'RefChannel')),
          % there might possibly be CTF reference channel information
          MEGChannel = good_channel(Channel,ones(size(GUI.ChannelFlag)),'MEG'); %good channels indexer

          if MEGChannel(1) ~= GoodChannel(1) % Critical case where the reference information should be stored in the first good MEG channel
             if(~isempty(Channel(MEGChannel(1)).RefChannel)), % there are reference sensors in Channel(1)
                Channel(GoodChannel(1)).RefChannel = Channel(MEGChannel(1)).RefChannel;
             end
             tmp = Channel(MEGChannel(1)).grad3;
             Channel(GoodChannel(1)).grad3 = tmp(:,GoodChannel-min(MEGChannel)+1);
             clear tmp
          end
       end
       DataFlag = 'MEG';
       HeadModel.Param = HeadModel.Param(GoodChannel);
       HeadModel.Function = HeadModel.Function(GoodChannel);
       Channel = Channel(GoodChannel);
    end

    % result is that GoodChannel(1) now has a loaded RefChannel and grad3 fields
    %  that are stripped down to just the number of good channels

    F = F(GoodChannel,:); % good channels only
    if(~isempty(Data.Projector)),
       Data.Projector = Data.Projector(GoodChannel,:);
    end
    for i = 1:length(HeadModel.SearchGain),
       if(~isstr(HeadModel.SearchGain{i})), % might be a filename of binary
          HeadModel.SearchGain{i} = HeadModel.SearchGain{i}(GoodChannel,:);
       end
    end
    % done altering the data and parameters
    % The search grid will be loaded and stripped below, depending on which
    %  model order is called.

    % Does the data contain an existing projector. Decompose it
    if(isfield(Data,'Projector')),
       if(isempty(Data.Projector)),
          A = [];
       else
          A = Data.Projector; % initialize
       end
    else % user did not give a projector
       A = [];
    end
    Ua = orth(A); % includes empty A

    % Data are assumed to have already passed through. Do so anyway
    if(~isempty(Ua)),
       F = F - Ua*(Ua'*F);
    end

    waitbar(0.75,hwait);
    drawnow

    % 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);
    end

    % form the decomposed form of the data
    UfSf = Uf(:,1:GUI.Rank)*Sf(1:GUI.Rank,1:GUI.Rank);

    % form the reduced rank data
    Fr = UfSf*Vf(:,1:GUI.Rank)';
    Fr2 = sum(Fr(:).^2); % total variance of data

    % finish off the waitbar, entering into the main run
    waitbar(1,hwait);
    drawnow
    close(hwait); % done loading, now the runs
    drawnow

    % pre-allocate some space. Will remap these to Results when done.
    [best(1:(NumSources*12))] = deal(struct('L',[],'Order',[],...
       'c',[],'O',[],'rank',[]));

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

    nOrder = 1;  % order indexer in GUI.Order, for NumSources.
    % We use nOrder as an indexer into StartingPoints
    Topo_i = 1;  % topology indexer, which source we are processing
    % Topo_i is only 1, but again we are allowing for future expansion
    %  for multiple runs of the topographies, such as in boot strapping
    %  or multiple restarts. Thus Topo_i represents the number of times
    %  we run this topography. In RAP-MUSIC terms, least-squares represents a single
    %  topography model.

    iOrder = GUI.Order{nOrder} + 2; % Order -1 has index of 1, etc.

    disp(sprintf('\n\n ************* Least-Squares Run at %s *************',datestr(now)));

    % First we will search over the searchgrid, to get a good localization point.
    % The complexity of the search depends on the number of NaN columns in the
    % GUI.StartingPoints. If only one, then we will exhaustively search.
    % If greater than one, then we will search up to PAIRS sets, before
    % launching the directed search

    disp(sprintf(...
       '\nPROCESSING Source %.0f at Order %.0f (%.0f) with regularization of %s\n',...
       Topo_i,GUI.Order{nOrder},length(GUI.StartingPoints),GUI.REG));

    % these are not that big, and useful to remap
    Lgrid = HeadModel.SearchGridLoc{iOrder};
    Ggrid = HeadModel.SearchGain{iOrder};
    if(isstr(Ggrid)), % it's a filename, load it
       fid = fopen(fullfile(Users.STUDIES,Ggrid),'rb','ieee-be');
       if(~fid),
          errordlg(sprintf('Unable to open %s',Ggrid));
          return
       end
       tempRows = fread(fid,1,'uint32');
       Ggrid = fread(fid,[tempRows, inf],'float32');

       fclose(fid);
       Ggrid = Ggrid(GoodChannel,:); % strip out the bad channels
    end

    nLgrid = size(Lgrid,2); % how many are there

    % quick double check for sanity
    if(nLgrid*DIMS(iOrder) ~= size(Ggrid,2)),
       errordlg(sprintf(...
          'Size of grid %s does not match model order (%s x %s)',...
          size(Ggrid,2),nLgrid,DIMS(iOrder)));
       return
    end

    % Scan the grid looking for a maximum
    % what is the rank of the signal subspace
    best(Topo_i).rank = GUI.Rank;

    % How many random sets do we need?
    LStart = [GUI.StartingPoints{:}]; % collect them all
    RandomSets = sum(isnan(LStart(1,:))); % number of cols that are NaN

    if(RandomSets == 0),
       % user has prespecified all starting points for the nonlinear search
       Sets = []; %nothing to search in the grid
    elseif(RandomSets == 1), % there is only a single random source to initially search
       % over a grid. We will use the entire Lgrid for this.
       Sets = [1:nLgrid]';
    elseif(RandomSets > 1),
       % create random unique sets in the grid
       iPair = 0;
       Sets = [];
       while(iPair < PAIRS),
          [ignore,next_set] = sort(randn(1,nLgrid));
          nLgrid_trim = nLgrid - rem(nLgrid,RandomSets); % nearest multiple of sets
          next_set = reshape(next_set(1:nLgrid_trim),...
             nLgrid_trim/RandomSets,RandomSets);
          iPair = iPair + size(next_set,1); % next set of pairs.
          % each row of next_set is uniquely paired, no twins.
          Sets = [Sets; next_set]; % no guarantee we haven't duplicated a set, but chance small
       end
    end

    % so now we have Sets, a matrix whose each row gives us index numbers
    %  into the Lgrid matrix and hence the precomputed gain matrix

    % The user has possibly preset a few points

    best(Topo_i).c = 0; %initialize the goodness of fit
    % snipped from nchoosek, adjusted notation
    SingleDim = DIMS(iOrder); % dimension of a single source
    SetDim = SingleDim*RandomSets; % how many columns per random source

    LStartndx = [1:size(LStart,2)]; % index to all points
    LStartndx(isnan(LStart(1,:))) = []; % delete the random points
    LStart = LStart(:,LStartndx); %just the fixed points
    % now form the forward model
    % CHEAT: Use first function only
    FixedG = feval(HeadModel.Function{1},LStart,Channel,HeadModel.Param,GUI.Order{nOrder}); % single source

    % now sweep through the sets

    drawnow
    hwait = waitbar(0,sprintf('Scanning through %.0f sets of %.0f sources',length(Sets),RandomSets));
    drawnow

    % preallocate the gain matrix
    ThisGain = zeros(size(Ggrid,1),NumSources*DIMS(iOrder)); % gain matrix of the model
    ThisGain(:,1:size(FixedG,2)) = FixedG; % map into the gain matrix
    BestSet = []; % initialize
    for i = 1:length(Sets), % may be possibly null

       this_set = Sets(i,:); % next set
       this_ndx = [-(DIMS(iOrder)-1):0]';
       % create indexer into the gain matrix
       this_ndx = this_ndx(:,ones(1,length(this_set))) + this_set(ones(DIMS(iOrder),1),:)*DIMS(iOrder);
       this_gain = Ggrid(:,this_ndx(:)); % the gain matrix for this set
       % Append to the fixed gain, preallocated above
       ThisGain(:,[1:size(this_gain,2)]+size(FixedG,2)) = this_gain;
       % project away any existing constraint
       if(~isempty(Ua))
          ThisGain = ThisGain - Ua*(Ua'*ThisGain);
       end
       % Has user requested column normalization of the grid
       if(GUI.Column_norm),
          % in over-constrained least-squares, the addition of a weighting matrix
          % has no effect, unlike weighted min norms.
          % just mucks up the units, but allows better truncated regularizers
          [ignore,ThisGain] = colnorm(ThisGain);
       end

       [Uthis,Sthis,Vthis] = svd(ThisGain,0); % decompose
       Uthis = regsubspace(GUI,Uthis,Sthis); % apply regularizers as needed
       ThisCost = Uthis'*UfSf; % matrix costs
       ThisCost = sum(ThisCost(:).^2); % this cost

       if(ThisCost > best(Topo_i).c), % new winner
          best(Topo_i).c = ThisCost;
          BestSet = this_set;
       end
       if(~rem(i,250)),
          waitbar(i/length(Sets),hwait);
          drawnow
       end
    end % for all sets

    close(hwait)
    drawnow
    % So we now have a best starting point. If Grid is dense enough, we should be close

    disp(sprintf(...
       'Initial solution: [%s], %.1f%% Explained',...
       str_repeater('[%.1f, %.1f, %.1f]','; ',[Lgrid(:,BestSet) LStart]*1000),best(Topo_i).c/Fr2*100));
    hwait = waitbar(.66,... % synthesize a percentage done already
       sprintf(...
       'Initial solution: [%s], %.1f%% Explained',...
       str_repeater('[%.1f, %.1f, %.1f]','; ',[Lgrid(:,BestSet) LStart]*1000),best(Topo_i).c/Fr2*100));
    drawnow

    % Optimization call
    % optimize the solution found on the grid
    %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);
    % get the display off, and set the tolerance to sub millimeters
    BestSources = fminsearch(@min_least_squares_fit,...
       [Lgrid(:,BestSet) LStart],Options,HeadModel.Function,Channel,HeadModel.Param,GUI.Order{nOrder},...
       GUI,UfSf,Ua);

    close(hwait)
    drawnow

    % we have an improved localization
    new_cost = -min_least_squares_fit(...
       BestSources,HeadModel.Function,Channel,HeadModel.Param,GUI.Order{nOrder},...
       GUI,UfSf,Ua);
    disp(sprintf(...
       'Improved solution: [%s], %.1f%% Var Explained',...
       str_repeater('[%.1f, %.1f, %.1f]','; ',BestSources*1000),new_cost/Fr2*100));

    if(new_cost < best(Topo_i).c),
       % somehown fmins came back more poorly than the grid. probably an error somewhere
       disp(' ')
       disp('Why is the new corr less than the gridded?')
       disp(sprintf('New corr: %.2f%%, gridded: %.2f%%',...
          new_corr, best(Topo_i).c));
       disp('Pausing at keyboard for investigation')
       keyboard
    end
    % assign the new correlation to the best of the results
    best(Topo_i).c = new_cost;
    best(Topo_i).L = BestSources;

    % CHEAT: use the first name only
    % recompute to get the corresponding orientation
    Gfinal = feval(HeadModel.Function{1},BestSources,Channel,HeadModel.Param,GUI.Order{nOrder}); % all sources

    % Extract the time series

    if(~isempty(Ua)),
       PpGfinal = Gfinal - Ua*(Ua'*Gfinal); % project away from existing subspace
    else
       PpGfinal = Gfinal;
    end

    if(GUI.Column_norm),
       % in over-constrained least-squares, the addition of a weighting matrix
       % has no effect, unlike weighted min norms.
       % just mucks up the units, but allows better truncated regularizers
       [ignore,PpGfinal] = colnorm(PpGfinal);
    end

    [Ug,Sg,Vg] = svd(PpGfinal,0);

    Ug = regsubspace(GUI,Ug,Sg); % apply regularizers as needed

    RankG = size(Ug,2); % the rank that survived regularization

    % now use Gfinal at this truncated subspace
    [Ug,Sg,Vg] = svd(Gfinal,0);

    switch GUI.REG % regularization switch
    case {'Condition','Energy','None'}
       % space truncated, directly form time series
       S = Vg(:,1:RankG)*(inv(Sg(1:RankG,1:RankG))*(Ug(:,1:RankG)'*Fr));

    case 'Tikhonov'
       % Tikhonov regularization
       % Note that regularizer was NOT applied in the maximization of cost
       Lambda = Sg(1)/GUI.Tikhonov; % relative condition setting
       if(RankG > 1),
          Sg = diag(Sg(1:RankG,1:RankG));
       else
          Sg = Sg(1); % rank one space
       end
       Fa = Sg ./ (Sg.^2 + Lambda^2); % regularized inverted singular values

       S = Vg(:,1:RankG)*(diag(Fa)*(Ug(:,1:RankG)'*Fr)); % regularized time series

    otherwise

      error(sprintf('Unknown regularization method: %s',GUI.REG))

    end % regularization switch

    % now forward model is Gfinal * S'.
    %  Final mopping up is to break each solution into its orthogonal fixed dipoles and time series
    %  along with proper units.

    % By convention, arrange time series in column
    S = S';

    Results = struct('Comment',[],'Function',[],'StudySubject',[],...
       'DataFlag',[],'GUI',GUI,'Date',[],'Subject',[],'Study',[],'Time',[],...
       'ChannelFlag',[],'NoiseCov',[],'SourceCov',[],...
       'Projector',[],'SourceOrder',[],'SourceCorrelation',[],...
       'SourceOrientation',[],'ModelGain',[],...
       'IndepTopo',[],'TimeSeries',[],'PatchNdx',[],...
       'PatchAmp',[],'ImageGridAmp',[],'ImageGridTime',[],'Fsynth',[]);

    Results.Fsynth = Gfinal*S'; % synthesized model

    Results.SourceLoc = cell(1,size(Gfinal,2)); % one source per column
    Results.SourceOrientation = Results.SourceLoc; % corresponding orientations
    Results.SourceCorrelation = zeros(1,size(Gfinal,2)); % variance explained by model
    Results.SourceOrder = zeros(1,size(Gfinal,2)); % model order of each source
    Results.TimeSeries = zeros(size(S)); % new time series
    Results.IndepTopo = zeros(size(Gfinal)); % "independent" topographies

    iS = 0; % indexer into SourceLoc
    for i = 1:size(BestSources,2), % for each rotating source
       ndx = [-(DIMS(iOrder)-1):0]+i*DIMS(iOrder); % index into model columns
       [Us,Ss,Vs] = svd(S(:,ndx),0); % decompose the associated time series
       % Us*Ss gives the new time series, in directions Vs
       for j = 1:min(DIMS(iOrder),size(Us,2)), % rank may be less than model
          Results.SourceLoc{iS+j} = BestSources(:,i); % same source location
          Results.SourceOrientation{iS+j} = Vs(:,j); % fixed orientation
          Results.SourceOrder(iS+j) = GUI.Order{nOrder}; % model order
          Results.IndepTopo(:,iS+j) = Gfinal(:,ndx)*Vs(:,j); % this orientation
          Results.TimeSeries(:,iS+j) = Us(:,j)*Ss(j,j); % new time series
          % Calculate the percent variance explained by this topography
          Results.SourceCorrelation(iS+j) = ...
             fronorm(Results.IndepTopo(:,iS+j)*Results.TimeSeries(:,iS+j)')^2 /Fr2;
          % now by convention, we actually scale time series by the column norm
          %  That gives a time series that visually gives the forward contribution
          %  The indep topo is not adjusted, such that the process is reversible
          Results.TimeSeries(:,iS+j) = Results.TimeSeries(:,iS+j)*norm(Results.IndepTopo(:,iS+j));
       end
       iS = iS + DIMS(iOrder); % increment block indexer
    end
    % mop up. Some columns may be null if rank was less than model order
    for i = length(Results.SourceLoc):-1:1, % work backwards, since we are deleting
       if(isempty(Results.SourceLoc{i})),
          Results.SourceLoc(i) = [];
          Results.SourceOrientation(i) = [];
          Results.SourceOrder(i) = [];
          Results.IndepTopo(:,i) = [];
          Results.TimeSeries(:,i) = [];
          Results.SourceCorrelation(i) = [];
       end
    end

    Results.Comment = sprintf('Least-Squares, rank %.0f, %.0f sources',GUI.Rank,size(BestSources,2));
    Results.Date = datestr(datenum(now),0);
    Results.Subject = Subject;
    Results.Study = Study;
    Results.Time = [GUI.Segment];
    Results.ChannelFlag = GUI.ChannelFlag;
    Results.NoiseCov = [];
    Results.SourceCov = [];
    Results.Projector = Data.Projector;
    Results.ModelGain = [];
    Results.PatchNdx = [];
    Results.PatchAmp = [];
    Results.ImageGridAmp = [];
    Results.ImageGridTime =[];
    Results.Function = mfilename; % name of this calling routine
    Results.DataFlag = DataFlag;
    Results.StudySubject = StudySubject;
    Results.GUI = GUI;

    if(isfield(GUI,'Results')), %user provided such a string
       if(~isempty(GUI.Results)), % and it's not empty
          if(0), % old way, inflexible
             save(fullfile(Users.STUDIES,GUI.Results),...
                'DataFlag','Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
                'ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder','SourceCorrelation',...
                'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
                'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
          else
             % new way,adapts to Results
             FIELDNAMES = fieldnames(Results); % list of the names
             eval(sprintf('%s = getfield(Results,''%s'');',FIELDNAMES{1},FIELDNAMES{1})); % get first field
             save(fullfile(Users.STUDIES,GUI.Results),FIELDNAMES{1}); % save it
             for i = 2:length(FIELDNAMES), % repeat for all additional
                eval(sprintf('%s = getfield(Results,''%s'');',FIELDNAMES{i},FIELDNAMES{i})); % get first field
                save(fullfile(Users.STUDIES,GUI.Results),FIELDNAMES{i},'-append'); % save it
             end
          end
          disp(sprintf('Wrote results to %s',GUI.Results)); % to the command window
          msgbox(sprintf('Wrote results to %s',GUI.Results),'Least Squares Completed'); % don't force modal
       end
    end

    if(GUI.DisplayGraphics),

       % let user see a solution
       tResults = Results;
       [ignore,tResults.IndepTopo] = colnorm(tResults.IndepTopo); % feature of the results update
       results_update(struct('F',Fr,'Time',Data.Time),tResults);

       drawnow

       if(hfdepthgauge), %is there a figure to draw to?
          depth_point = BestSources;
          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
    end % if user wanted graphics

    return

    function e = min_least_squares_fit(L,Function,Channel,Param,Order,GUI,UfSf,Ua);
    %MIN_LEAST_SQUARES_FIT minimize the least squares error using the a specified gain model
    % LIMITED: ONLY LOOKS AT THE FIRST FUNCTION CALL AND APPLIES EQUALL TO ALL CHANNELS!
    % function e = min_least_squares_fit(L,Function,Channel,Param,Order,GUI,UfSf,Ua);
    % call as fminsearch(...
    %   'min_least_squares_fit',L,[],[],Function,Channel,Param,Order,GUI,UfSf,Ua);
    % Ua must already be orthogonal
    % UfSf is the reduced rank data, USV'=F, with V tossed.
    % 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"

    % ----------------------------- Script History ---------------------------------
    % John C. Mosher, Ph.D.
    % 19-May-2004 JCM  Comments cleaning
    % ----------------------------- Script History ---------------------------------

    % based on min_fgain_gui for rap_music.

    % generate the gain matrix

    %%%%% 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
      G= feval(Function{1},L,Channel,Param,Order);
    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

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

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

    % squared frobenius norm of the fit
    c = Ug'*UfSf;
    c = sum(c(:).^2);

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

    return

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON