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
