[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
This function is called by
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 Oct 12 12:05:14 2004
Cross-Directory links are: ON