[Master Index]
[Index for Toolbox]
rapmusic_gui
(Toolbox/rapmusic_gui.m in BrainStorm 2.0 (Alpha))
Function Synopsis
Results = rapmusic_gui(StudySubject,GUI);
Help Text
RAPMUSIC_GUI - Execute RAPMUSIC using GUI inputs, BrainStorm MMII version
function Results = rapmusic_gui(StudySubject,GUI);
StudySubject is the conventional structure of strings,
StudySubject = find_brainstorm_structure(<studyname>).
GUI is structure of parameters needed to run this function,
see below.
All strings represent .mat files that are selectively loaded from disk.
The string names should represent fully qualified filenames.
Parameters in GUI are:
.DataName, string of actual data file used
.Results, optional file to write results
.Segment, an index vector giving the column index numbers to use in
the spatiotemporal data matrix.
.Order, ordered vector of multipolar orders to process, e.g. [-1 0 1]
.OrderHeadModel, one-one with .Order, the referential filename to the
corresponding HeadModel to use with each Order
.Rank, rank to use, e.g. 10
.Corr, correlation to use in processing, e.g. 0.98
.ChannelFlag, a vector of Flags in the data file to process, See Data.ChannelFlag
.DisplayGraphics, if non-zero, plots information as the algorithm runs
The following three structure fields are optional for regularization purposes.
They may be null or missing to represent unused. If multiple fields are given,
then precedence is given in the order given below.
.Condition, condition number to use in truncation, e.g. 100
.Energy, fractional energy to use in truncation, e.g. .95
.Tikhonov, Tikhonov condition number to use in regularization
If all are null, no regularization is performed in the RAP-MUSIC loops.
Another flag that effects regularization is
.Column_norm, if non-zero then columns of gain matrix are normalized in the
original model
If GUI.Results is non-null,then writes results to that file.
Cross-Reference Information
This function calls
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\rapmusic_gui.m
function Results = rapmusic_gui(StudySubject,GUI);
%RAPMUSIC_GUI - Execute RAPMUSIC using GUI inputs, BrainStorm MMII version
% function Results = rapmusic_gui(StudySubject,GUI);
% StudySubject is the conventional structure of strings,
% StudySubject = find_brainstorm_structure(<studyname>).
% GUI is structure of parameters needed to run this function,
% see below.
% All strings represent .mat files that are selectively loaded from disk.
% The string names should represent fully qualified filenames.
% Parameters in GUI are:
% .DataName, string of actual data file used
% .Results, optional file to write results
% .Segment, an index vector giving the column index numbers to use in
% the spatiotemporal data matrix.
% .Order, ordered vector of multipolar orders to process, e.g. [-1 0 1]
% .OrderHeadModel, one-one with .Order, the referential filename to the
% corresponding HeadModel to use with each Order
% .Rank, rank to use, e.g. 10
% .Corr, correlation to use in processing, e.g. 0.98
% .ChannelFlag, a vector of Flags in the data file to process, See Data.ChannelFlag
% .DisplayGraphics, if non-zero, plots information as the algorithm runs
% The following three structure fields are optional for regularization purposes.
% They may be null or missing to represent unused. If multiple fields are given,
% then precedence is given in the order given below.
% .Condition, condition number to use in truncation, e.g. 100
% .Energy, fractional energy to use in truncation, e.g. .95
% .Tikhonov, Tikhonov condition number to use in regularization
% If all are null, no regularization is performed in the RAP-MUSIC loops.
% Another flag that effects regularization is
% .Column_norm, if non-zero then columns of gain matrix are normalized in the
% original model
% If GUI.Results is non-null,then writes results to that file.
%<autobegin> ---------------------- 14-Jun-2004 17:12:32 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Inverse Modeling
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\bst_color_scheme.m
% toolbox\bst_imaging_display.m
% toolbox\bst_layout.m
% toolbox\bst_message_window.m
% toolbox\bst_mriviewer.m
% toolbox\bst_static_taskbar.m
% toolbox\bst_wavedata_display.m
% toolbox\colnorm.m
% toolbox\extract_channels.m
% toolbox\get_user_directory.m
% toolbox\good_channel.m
% toolbox\load_brainstorm_file.m
% toolbox\norcol.m
% toolbox\regcheck.m
% toolbox\regsubspace.m
% toolbox\results_visualization.m
% toolbox\save_fieldnames.m
% toolbox\str_repeater.m
% toolbox\subcorr_gui.m
%
% Subfunctions in this file, in order of occurrence in file:
% P = combs(v,m)
% e = min_fgain_gui(L,Function,Channel,Param,Order,GUI,Uf,Ua,Whitener,imegsens,irefsens,ieegsens);
%
% Application data and their calls in this file:
% 'TileType'
% 'visfig'
%
% setappdata(CorticEnvFigure,'TileType','T')
% setappdata(hTimeSeries,'TileType','T')
% setappdata(hTimeSeries,'visfig',visfig);
%
% visfig = getappdata(hTimeSeries,'visfig');
%
% At Check-in: $Author: Mosher $ $Revision: 44 $ $Date: 6/14/04 3:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-2004
%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
% University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
% Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
% CNRS, Hopital de la Salpetriere, Paris, France
%
% See BrainStorm website at http://neuroimage.usc.edu for further information.
%
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 14-Jun-2004 17:12:32 -----------------------
% John C. Mosher, Ph.D.
% Sylvain Baillet, Ph.D.
% ---------------------------------------------------------------------------------
% 6/03/99 JCM, consolidating the regularization functions
% 6/09/99 JCM, updating for .DataName and .Results
% 6/24/99 JCM, updating for visualization of music metric
% 10/18/99 JCM, adding the 2 synchronous dipole model
% 10/22/99 SB, updating tessellation file management fo visualization
% JCM 10/26 updating 2 synchronous model handling
% JCM 11/9/99 updating referential filenames, dropping the SearchTess stuff.
% JCM 12/2/99 replacing fmins with fminsearch, fmins now outdated.
% SB Dec-Jan CTF gradiometer stuff.
% JCM 1/31/00 synchronizing versions
% JCM 2/16/00 Segment is now an index vector, display graphics handled by
% results_update
% JCM 2/23/00 some cleanup
% JCM 3/19/00 Added in Baillet's fixes for CTF reference sensors
% JCM 3/20/00 Compatibility with results_update
% SB 4/27/01 Little clean-up and update on bad-channel referencing with
% CTF channel structure
% JCM 2/25/2002 used function handles for subfunction for fminsearch, replaced "disp" with bst_message_window
% JCM 10-May-2002 major changes to headmodel, to conform to MMII approach
% i.e. no more "SearchGain", just "Gain", etc.
% Have to rethink the different gain matrices as well. The
% magnetic dipole is on hiatus, we have just the current dipole
% and current multipole, orders -1 and 1.
% JCM 13-May-2002 Added GUI.OrderHeadModel, which is a cell-array of
% referential strings to the HeadModel corresponding to GUI.Order
% JCM 14-May-2002 Moving HeadModel around to allow that each SourceOrder is called
% separately. SB and JCM also updated the HeadModel to have fields
% SourceOrder,HeadModelType. RAP-MUSIC doesn't actually care which
% HeadModelType search grid is called, but the SourceOrder is crucial.
% 24-Jun-2002 JCM Comments changes, upped Pairs to 5000, fixed bug with Data Projector
% and the displaying of results. Some of A is now DP (Data Projector)
% to distinguish those components of the new model from those of DP.
% 21-Oct-2002 JCM Major: Changed projection operation. Rather than project the signal
% subspace away from the found sources, we project away the *data*, then
% recompute the signal subspace. In ProjOption = 1 (hardwired just
% below these comments), we decrement the signal subspace rank by one for
% each recursion. In ProjOption = 2, we keep the rank the same. In
% theory, the full rank subspace will become garbage anyhow.
% In ProjOption = 0, now deprecated, we keep the old scheme, for
% legacy testing.
% JCM 31-Oct-2002 Bug handling for missing UDP, (orthogonal data projector matrix).
% SB 24-Oct-2003 Changed call to bst_wave_display
% SB 05-Mar-2004 Changed call to bst_wave_display (again)
% JCM 03-May-2004 Cleaning comments. bst_mriviewer is used to display the
% final solution. bst_imaging_dispay and
% bst_wavedata_display are also called.
% ---------------------------------------------------------------------------------
% ------ HARDWIRED PARAMETERS ---------------
% maximum number of random pairs to search in the initial synchronous
% grid. Set based on speed performance.
PAIRS = 5000;
% Projection Options, see comments above from 21-Oct-2002 JCM.
% Set ProjOption to 0, 1, or 2
% 0 deprecated, the old style of projecting just the signal subspace:
% 1 recommended, project the data, then find the signal subspace as the new
% dominant subspace vectors, where the rank of the new space has been
% decremented each time:
% 2 experimental, project the data, then find the signal subspace, but
% keep the rank the same:
ProjOption = 1;
% DIMS is probably deprecated, now that each head model is separate.
DIMS = [3 3 12]; % the number of linear parameters per multipolar oders -1, 0, and 1.
% --------- END HARDWIRED ------------------
% keep the user from getting bored.
bst_message_window('Loading information to run RAP-MUSIC . . .');
Users = get_user_directory; % the default search paths
% load the head model information
if isempty(GUI.OrderHeadModel),
% User must compute a head model prior to running RAP_MUSIC
errordlg('Please compute a head model prior to running RAP_MUSIC')
return
end
% which subject is this, informational
Subject = load_brainstorm_file(StudySubject.Subject);
if(0), % JCM 10-May-2002, won't display tesselated results anymore
% Subject Information
if isempty(Subject.Tesselation) % User must assign a Tesselation file to the subject -> a dummy one if necessary
errordlg('Please assign first a Tesselation file to the Subject file')
return
end
if(GUI.DisplayGraphics)
% Tesselation information, setup up FV for viewing the source solutions.
% Every subject must have a tesselation for use, even if surrogate.
SubjectTess = load_brainstorm_file(Subject.Tesselation);
%Find cortical tesselations
% Look for tessellations beguinning with 'cortex' as a keyword -> update documentation
icrtx = strmatch('cortex',lower(SubjectTess.Comment));
if(isempty(icrtx)),
errordlg('Need a tesslation that begins with ''cortex''. Return to Head Modeler');
return
end
% If multiple cortical files, Look for the smallest tessellation.
if(length(icrtx) > 1)
k = 0;
for i = icrtx(:)',
k = k+1;
siz(k) = size(SubjectTess.Faces{i},1);
end
[ignore, imin] = min(siz);
else
imin = 1; % only one cortical tesselation
end
% assign the cortical tesselation to the structure
FV = struct('faces',SubjectTess.Faces{icrtx(imin)},...
'vertices',SubjectTess.Vertices{icrtx(imin)}');
clear SubjectTess % don't need the rest of the Tesselations
% Create a figure for viewing the source results
hfdepthgauge = figure; % new figure, used to show depthgauges into cortex
set(hfdepthgauge,'units','normalized');
set(hfdepthgauge,'position',[.675 .55 .31 .31]); % upper right corner
set(hfdepthgauge,'Name',GUI.Results);
[hfdepthgauge,hpatch] = view_surface(GUI.Results,FV.faces,FV.vertices);
set(hpatch,'FaceLighting','flat'); % much faster
uimenu('Label','Toggle Face','callback','toggleface(get(gcbo,''UserData''))','UserData',hpatch);
uimenu('Label','Toggle Lighting','callback','togglelight(get(gcbo,''UserData''))','UserData',hpatch);
rotate3d on
axis tight vis3d
drawnow
end
end % no more viewing cortical tesselation
% regularization check
GUI = regcheck(GUI); % returns GUI.REG set to the appropriate field
Study = load_brainstorm_file(StudySubject.Study); % and the study information
% don't need SubjectImage, SubjectTess, SubjectTriConn, SubjectVertConn
% now trim the data to just the desired channels
% Use the first head model to establish desired channels
[HeadModel, Channel, GoodChannel, DataFlag] = ...
extract_channels(GUI.OrderHeadModel{1},StudySubject,GUI);
% GoodChannel represents that subset of channels desired for analysis
% Channel and HeadModel have already been trimmed to this set
% find indices of all MEG and EEG channels
imegsens = good_channel(Channel,[],'MEG');
irefsens = good_channel(Channel,[],'MEG REF'); % may return empty
ieegsens = good_channel(Channel,[],'EEG');
% Load the data set
Data = load_brainstorm_file(GUI.DataName);
if ~isempty(Data.Projector),
if(size(Data.Projector,1) ~= length(GoodChannel)),
bst_message_window('wrap',{...
'Projector not same row dimension as GoodChannels.',...
'Projector removed'});
Data.Projector = [];
cd(Users.STUDIES)
save_fieldnames(Data,GUI.DataName); % save the empty projector
end
end
% Trim the data to the requested time segments
F = Data.F(:,GUI.Segment);
F = F(GoodChannel,:); % good channels only
Time = Data.Time(GUI.Segment);
% done altering the data and parameters
% did the user provide an existing control. Decompose it
% JCM 24-Jun-2002, switched this A to DP (Data Projector),
% to keep clear what is DP and what is the Independent Topographies
% generated by this fit to the data.
if isfield(Data,'Projector'),
if(isempty(Data.Projector)),
DP = [];
else
DP = Data.Projector; % handy shorthand
end
else % user did not give a projector
DP = [];
end
% Is there a noise covariance for a whitener?
if ~isempty(Data.NoiseCov),
Whitener = Data.NoiseCov(GoodChannel,:);
Whitener = pinv(Whitener);
% now whiten the data and the precomputed gain
F = Whitener*F;
HeadModel.Gain{1} = Whitener*double(HeadModel.Gain{1});
else
Whitener = []; %
end
% The data could have already been placed into a subspace
% away from the projector.
% . Noise: Should already be in a lower dimension
% . Other solution: Projector represents prior solution, so
% Data may not be projected.
% Let's ensure.
if(~isempty(DP)), % was there a Data Projector?
UDP = orth(DP); % make sure that DP is orthogonalized
F = F - UDP*(UDP'*F); % project data away from projector
else
UDP = []; % orthog DP is also empty
end
% Now the data exist in a lower dimensional subspace away from the
% projector. Proceed as normal, building a topographies matrix.
%
% Initialize the Independent Topographies matrix separate
% from the DP matrix. In the code that follows, Ua effectively represents
% orth([DP A]), but A starts off as null.
A = []; % A initialized for the loops
Ua = orth([DP A]); % DP and A orthogonalized, works for null case as well
Aorg = []; % the final topographies, no spatial filters applied
bst_message_window('overwrite','Head model loaded, loading data . . .');
% Preprocess the Data
% User has asked for an rank RANK subspace search.
% Form the signal subspace
if(size(F,1) >= size(F,2)), % skinny matrix
[Uf,Sf,Vf] = svd(F,0); % economy
else % wide matrix
[Vf,Sf,Uf] = svd(F.',0); % direct transpose
end
% form the reduced rank data
Fr = Uf(:,1:GUI.Rank)*(Sf(1:GUI.Rank,1:GUI.Rank)*Vf(:,1:GUI.Rank)');
% form UPF as the projected signal subspace
UPF = Uf(:,1:GUI.Rank);
% form the structure for the answers in the loops below.
% pre-allocate some space. No problem if exceeded, just
% slows down some. Will remap these to Results when done.
[best(1:100)] = deal(struct('L',[],'Order',[],...
'c',[],'O',[],'rank',[]));
% Entering into this loop, we have the predefined topology matrix A and Ua, possibly null
% if the user gave no initial projector. We will build on this matrix.
nOrder = 1; % order indexer in GUI.Order. Possible to search in any order, but
% logically we only anticipate increasing modeling order, such as [-1 0 1];
% New MMII code, load particular headmodel here
% The first headmodel was already loaded above
Topo_i = 1; % topology indexer, which source we are processing
bst_message_window('append',{'','',...
sprintf('************* Run at %s *************',datestr(now))});
while((nOrder <= length(GUI.Order)) & (Topo_i <= GUI.Rank)), % multipolar order maximum
Order = GUI.Order(nOrder); % The source order
% handle synchronous source flags.
% Orders -1, 0, and 1 are current, magnetic, and 1st-order multipolar
% Orders 2, 3, 4 are synchronous pairs for -1, 0, and 1.
% Not elegant, but workable
% iOrder is used for indexing the search grids and orders, 1:6 one-to-one with -1:4
switch Order
case {-1,0,1}
% we're fine, but set Sync to 1
Sync = 1; % single sources
iOrder = Order + 2; % Order -1 is 1, Order 1 is 3, array indexer
case 2 % synchronous current dipole
Order = -1; % they are really dipoles
Sync = 2; % want two synchronous sources
% FIX, we really need a two dipole grid, not the denser single dipole grid
iOrder = 1; % Use the synchronous current dipole grid
case 3 % synchronous magnetic dipole
Order = 0; % magnetic dipoles
Sync = 2;
iOrder = 2; % Use the synchronous magnetic dipole grid
case 4 % synchronous multipolar
Order = 1; % 1st-order multipoles
Sync = 2;
iOrder = 3; % Use the multipolar dipole grid
otherwise % -1, 0, 1 orders
errormsg(sprintf('Don''t know this order: %.0f',Order));
return
end
bst_message_window('append',{'',sprintf(...
'PROCESSING Source %.0f at Order %.0f (%.0f)',Topo_i,Order,Sync),...
sprintf('with regularization of %s',GUI.REG),''});
% these are not that big, and useful to remap
if(0), % old calls
Lgrid = double(HeadModel.GridLoc{iOrder});
Ggrid = double(HeadModel.Gain{iOrder});
else % new MMII calls
Lgrid = double(HeadModel.GridLoc{1});
Ggrid = double(HeadModel.Gain{1});
end
nLgrid = size(Lgrid,2); % how many are there
% Has user requested column normalization of the grid
if(GUI.Column_norm),
% CHEAT, May 2002, to reimplement this, we need to store the "ignored" data
[ignore,Ggrid] = colnorm(Ggrid);
end
% first scan the grid looking for a maximum
% what is the rank of the signal subspace
best(Topo_i).rank = size(UPF,2);
% project the grid away from the existing model subspace
if(~isempty(Ua)), % may be empty on the first iteration
PGgrid = Ggrid - Ua*(Ua'*Ggrid);
else % ua is empty
PGgrid = Ggrid; % grid remains unprojected
end
if(Sync == 1), % searching for single sources
% form grid of decomposed gain matrices
UGrid = zeros(size(Ggrid)); % decomposed grid
% now must form and possibly truncate each grid point
Dim = DIMS(iOrder); % how many columns per source
for i = 1:nLgrid, % for each source location
ndx = i*Dim+[(1-Dim):0];
[Utemp,Stemp,Vtemp] = svd(PGgrid(:,ndx),0);
Utemp = regsubspace(GUI,Utemp,Stemp); % trim to conditioned subspace
% keep the null space columns to zero
UGrid(:,ndx(1:size(Utemp,2))) = Utemp;
end % for each source location
% We will bypass calling subcorr_gui here, for speed. As of 6/13/99, all
% forms of regularization use the unaltered subspaces formed above.
% UPF is the signal subspace.
% UPF is already orthogonal and possibly rank reduced.
C = UPF'*UGrid; % matrix inner product against the grid
c = zeros(1,nLgrid); % the correlations
for i = 1:nLgrid, % for each source location
cndx = i*Dim+[(1-Dim):0]; % which submatrix
tmp = svd(C(:,cndx)); % decompose it
c(i) = tmp(1); % maximum correlation between UPF and this source
end
[max_c,Di] = max(c); % maximum correlation over the grid
best(Topo_i).L = Lgrid(:,Di); % location of best source
best(Topo_i).c = max_c; % corresponding correlation
best(Topo_i).Order = Order; % model order
bst_message_window('append',sprintf(...
'Initial solution on Grid: [%s], %.1f%%',...
str_repeater('%.1f',',',best(Topo_i).L*1000),best(Topo_i).c*100));
bst_message_window('Optimizing in a nonlinear directed search . . .');
% Optimization call
% optimize the solution found on the grid
if(0), % outdated fmins call, as of Matlab 5.3
best(Topo_i).L = fmins(@min_fgain_gui,...
best(Topo_i).L,[],[],HeadModel.Function,Channel,HeadModel.Param,Order,...
GUI,UPF,Ua,Whitener);
else % the newer call
%See OPTIMSET for details. FMINSEARCH uses
%these options: Display, TolX, TolFun, MaxFunEvals, and MaxIter.
% We're dealing in meters, and 1e-4 is tenths of a millimeter.
% Let's go for smaller. TolFun is correlation keep to six places
Options = optimset('Display','off','TolX',1e-6,'TolFun',1e-6);
% get the display off, and set the tolerance to sub millimeters
best(Topo_i).L = fminsearch(@min_fgain_gui,...
best(Topo_i).L,Options,HeadModel.Function,Channel,HeadModel.Param,Order,...
GUI,UPF,Ua,Whitener,imegsens,irefsens, ieegsens);
end
elseif(Sync > 1), % multiple sources
best(Topo_i).c = 0; %initialize the correlation
% snipped from nchoosek, adjusted notation
SingleDim = DIMS(iOrder); % dimension of a single source
SetDim = SingleDim*Sync; % how many columns per source
% FIX: area to research, better gridding.
% we have two means of initializing the directed search, either exhaustively
% testing every combination of grid points, or seeking random pairs
% Exhaustive can be too demanding if there are too many grid points.
% For now, stick to random pairs.
% pick a switch with a comment line
% switch 'Exhaustive' % exhaustively test all combinations
switch 'RandomPair' % randomly select pairs for testing
case 'Exhaustive'
bst_message_window('append',sprintf('Searching through %.0f exhaustive combinations. . .',...
nchoosek(nLgrid,Sync)));
vIndices = [1:nLgrid]; % index to all dipole numbers, row vec
for Comb_idx = 1:nLgrid-Sync+1, % indexing through the combinations
Q = combs(vIndices(Comb_idx+1:nLgrid),Sync-1); % local function call
next_set = [vIndices(ones(size(Q,1),1),Comb_idx) Q];
% code below virtually identical to RandomPair
% bst_message_window('append',sprintf('Another %.0f sets',size(next_set,1)));
% next_set has the next set of possible combinations in the grid
% form grid of decomposed gain matrices
UGrid = zeros(size(Ggrid,1),SetDim*size(next_set,1)); % decomposed grid
for i = 1:size(next_set,1), % for each source combination
ndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
sndx = zeros(1,SetDim); % set indexer
for j = 1:length(next_set(i,:)), % each item in the set
sndx([(1-SingleDim):0]+j*SingleDim) = [(1-SingleDim):0] + next_set(i,j)*SingleDim;
end
Utemp = PGgrid(:,sndx); % extract the desired columns
[Utemp,Stemp,Vtemp] = svd(Utemp,0); %decompose
Utemp = regsubspace(GUI,Utemp,Stemp); % trim to conditioned subspace
% keep the null space columns to zero
UGrid(:,ndx(1:size(Utemp,2))) = Utemp;
end % for each source combination
C = UPF'*UGrid; % inner product
c = zeros(1,size(next_set,1));
for i = 1:size(next_set,1),
cndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
tmp = svd(C(:,cndx));
c(i) = tmp(1);
end
[max_c,Di] = max(c); % maximum correlation in this set
if(max_c > best(Topo_i).c), % better combination
best(Topo_i).c = max_c;
best(Topo_i).L = Lgrid(:,next_set(Di,:));
best(Topo_i).Order = GUI.Order(nOrder);
end
end % next set of combination indices
case 'RandomPair'
bst_message_window('append',sprintf('Searching through at least %.0f random combinations. . .',PAIRS));
iPair = 0; % indexer into the pairs
while (iPair < PAIRS);
% create random unique pairs in the grid
[ignore,next_set] = sort(randn(1,nLgrid));
nLgrid_trim = nLgrid - rem(nLgrid,Sync); % nearest multiple
next_set = reshape(next_set(1:nLgrid_trim),nLgrid_trim/Sync,Sync);
iPair = iPair + size(next_set,1); % next set of pairs.
% each row of next_set is uniquely paired, no twins.
% code below virtually identical to Exhaustive
% bst_message_window('append',sprintf('Another %.0f sets',size(next_set,1)));
% next_set has the next set of possible combinations in the grid
% form grid of decomposed gain matrices
UGrid = zeros(size(Ggrid,1),SetDim*size(next_set,1)); % decomposed grid
for i = 1:size(next_set,1), % for each source combination
ndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
sndx = zeros(1,SetDim); % set indexer
for j = 1:length(next_set(i,:)), % each item in the set
sndx([(1-SingleDim):0]+j*SingleDim) = [(1-SingleDim):0] + next_set(i,j)*SingleDim;
end
Utemp = PGgrid(:,sndx); % extract the desired columns
[Utemp,Stemp,Vtemp] = svd(Utemp,0); %decompose
Utemp = regsubspace(GUI,Utemp,Stemp); % trim to conditioned subspace
% keep the null space columns to zero
UGrid(:,ndx(1:size(Utemp,2))) = Utemp;
end % for each source combination
C = UPF'*UGrid; % inner product
c = zeros(1,size(next_set,1));
for i = 1:size(next_set,1),
cndx = i*SetDim+[(1-SetDim):0]; % index in the sets matrix
tmp = svd(C(:,cndx));
c(i) = tmp(1);
end
[max_c,Di] = max(c); % maximum correlation in this set
if(max_c > best(Topo_i).c), % better combination
best(Topo_i).c = max_c;
best(Topo_i).L = Lgrid(:,next_set(Di,:));
best(Topo_i).Order = GUI.Order(nOrder);
end
end % next set of combination indices
end % case switch of random vs. exhaustive searching
bst_message_window('overwrite',sprintf(...
'Initial solution on Grid: [%s], %.1f%%',...
str_repeater('%.1f',',',best(Topo_i).L*1000),best(Topo_i).c*100));
bst_message_window('Optimizing in a nonlinear directed search . . .');
% Optimization call
% optimize the solution found on the grid
if(0), % outdated fmins call, as of Matlab 5.3
best(Topo_i).L = fmins(@min_fgain_gui,...
best(Topo_i).L,[],[],HeadModel.Function,Channel,HeadModel.Param,Order,...
GUI,UPF,Ua,Whitener);
else % the newer call
%See OPTIMSET for details. FMINSEARCH uses
%these options: Display, TolX, TolFun, MaxFunEvals, and MaxIter.
% We're dealing in meters, and 1e-4 is tenths of a millimeter.
% Let's go for smaller.
Options = optimset('Display','off','TolX',1e-6,'TolFun',1e-6);
% get the display off, and set the tolerance to sub millimeters
best(Topo_i).L = fminsearch(@min_fgain_gui,...
best(Topo_i).L,Options,HeadModel.Function,Channel,HeadModel.Param,Order,...
GUI,UPF,Ua,Whitener,imegsens,irefsens,ieegsens);
end
end % handling the synchronization cases
% we have completed the localization, ready to see if it's any good
% we have an improved localization
new_corr = -min_fgain_gui(best(Topo_i).L,HeadModel.Function,...
Channel,HeadModel.Param,Order,GUI,UPF,Ua,Whitener,imegsens,irefsens, ieegsens);
bst_message_window('overwrite',sprintf(...
'Improved solution: [%s], %.1f%%',...
str_repeater('%.1f',',',best(Topo_i).L*1000),new_corr*100));
if(new_corr < best(Topo_i).c-100*eps),
% somehown fmins came back more poorly than the grid. probably an error somewhere
bst_message_window('append',' ')
bst_message_window('append','Why is the new corr less than the gridded?')
bst_message_window('append',sprintf('New corr: %.2f%%, gridded: %.2f%%',...
100*new_corr, 100*best(Topo_i).c));
bst_message_window('append','Pausing at keyboard for investigation, ''dbcont'' to continue.')
keyboard
end
% assign the new correlation to the best of the results
best(Topo_i).c = new_corr;
if(best(Topo_i).c >= GUI.Corr) % we have an acceptable solution!
% CHEAT: use the first name only
% recompute to get the corresponding orientation
G1 = feval(HeadModel.Function{1},best(Topo_i).L,Channel,HeadModel.Param,Order,imegsens,irefsens, ieegsens); % single source
G1org = G1; % original topography, no whitener, no projector
if(GUI.Column_norm), % column normalization
[ignore,G1] = colnorm(G1);
end
if ~isempty(Whitener),
G1 = Whitener*G1;
end
if ~isempty(Ua), % if there already existing projectors
PpG1 = G1 - Ua*(Ua'*G1); % project away from existing subspace
else
PpG1 = G1;
end
% Note that here we get the orientation in the reduced subspace.
% get the associated orientation, regularized if need be.
[Co,X,Y] = subcorr_gui(PpG1,UPF,GUI);
best(Topo_i).O = X(:,1)/norm(X(:,1)); % best unit norm orientation
bst_message_window('append',{sprintf(...
'Orientation [%s]',...
str_repeater('%.2f',',',best(Topo_i).O)),...
sprintf('Rank %.0f Source Order %.0f',...
best(Topo_i).rank,Order)});
% need to form the current topography, append to existing
% and project the signal subspace away from it
A = [A G1*best(Topo_i).O];
Aorg = [Aorg G1org*best(Topo_i).O]; % unfiltered for visualization
Ua = orth([DP A]); % orthogonalize the existing space with the projector
% re-form UPF as the projected signal subspace
% 21 October 2002 Projection Options
% See Comments at top, plus hardwired section at top
switch ProjOption
case 0
% legacy, now a deprecated method
UPF = Uf(:,1:GUI.Rank); % get the original signal subspace
UPF = UPF - Ua*(Ua'*UPF); % project away existing subspace
% orthogonalize the reduced subspace
[UPF,SPF,VPF] = svd(UPF,0);
% %%%%%%%%%% NOT A CHEAT BUT LOOK AT THIS %%%%
% we've had a lot of discussions about the next line. Do we trim off the weakest space
% in the next recursion. Earlier concerns were about regularization
% JCM thoughts 2/23/00, why not run regsubspace?
switch 'reduce'
case 'reduce'
UPF = UPF(:,1:(end-Topo_i)); % remove the weakest vector spaces
case 'retain'
% do nothing, don't truncate
end
case {1,2}
% recommended 21 Oct 2002
% User has asked for an rank RANK subspace search.
% Form the signal subspace
ReducedF = F - Ua*(Ua'*F); % The data in the new subspace
if(size(ReducedF,1) >= size(ReducedF,2)), % skinny matrix
[UPF,ignore,ignore] = svd(ReducedF,0); % economy
else % wide matrix
[ignore,ignore,UPF] = svd(ReducedF.',0); % direct transpose
end
% form UPF from the signal subspace
UPF = UPF(:,1:GUI.Rank);
if(ProjOption == 1), % remove the weakest vector spaces from the subspace
UPF = UPF(:,1:(end-Topo_i)); % remove the weakest vector spaces
end
otherwise
error('Unknown Projection Method hardwired into code.');
bst_message_window('Contact brainstorm@sipi.usc.edu about error');
end
if(GUI.DisplayGraphics),
% setup the GUIs for the display of time series
global data % Get display parameters from Viewer
if Topo_i == 1
hTimeSeries = figure;
setappdata(hTimeSeries,'TileType','T')
bst_color_scheme(hTimeSeries);
bst_layout('align',hTimeSeries,2,2,4)
[ignore,ResultsFileName] = fileparts(GUI.Results);
set(hTimeSeries,'Name',ResultsFileName);
data.Results.Display.Handles.Figures.TimeSeries = hTimeSeries;
end
% Graphics option, could be a lot better:
% let user see a solution
% need to regularize
if(~isempty(UDP)),
Ap = A - UDP*(UDP'*A); % project model away
else
Ap = A; % no projector
end
% Data already projected away
S = [pinv(Ap)*Fr]'; % time series associated with the reduced rank data
Frsyn = Ap*S'; % the data synthesized in the whitened, projected space
Frem = Fr-Frsyn; % remainder
SourceLoc = cell(1,Topo_i);
SourceOrientation = cell(1,Topo_i);
for i = 1:Topo_i,
SourceLoc{i} = best(i).L; % map each solution into a cell.
SourceOrientation{i} = best(i).O;
end
% temporary results structure
if(0) % no longer displaying this way
tResults = struct('IndepTopo',Ap,...
'TimeSeries',S,'Time',GUI.Segment,'SourceOrder',[best(1:Topo_i).Order],...
'Comment',sprintf('RAP-MUSIC, rank %.0f, corr %.3f',GUI.Rank,GUI.Corr),...
'Subject',Subject,'Study',Study,'SourceCorrelation',[best(1:Topo_i).c],...
'SourceLoc',[],'SourceOrientation',[],'GUI',GUI,'Function',[]);
tResults.SourceLoc = SourceLoc; % would otherwise DEAL
tResults.SourceOrientation = SourceOrientation;
tResults.Function = mfilename; % name of this calling routine
% send the reduced rank time series. User threw away the rest
results_time_series('update',hTimeSeries,struct('F',Fr,'Time',Data.Time),tResults);
end
% -> Display source time series
D = num2cell(S',2);
for src = 1:length(D)
SrcName{src} = sprintf('Source %d',src);
end
GoodTypes = {Channel(1).Type};
chans = good_channel(Channel,[],GoodTypes{1});
DisplayStruct = data.Measures.Display;
%bst_layout('align',data.Results.Display.Figure(1),2,2,4);
OPTIONS = struct(...
'ModalityLabel',{GoodTypes},...
'WaveLabel',{SrcName},...
'Channel', data.Measures.Channel(data.Measures.Display.SelectedChannels{DisplayStruct.Modality}),...
'IndepTopo',Aorg,...
'Time', Time, ...
'MoveCursor','on'...
);
OPTIONS.SelectedChannels = data.Measures.Display.SelectedChannels(data.Measures.Display.Modality);
OPTIONS.DataType = 'source';
figure(hTimeSeries), clf
figs = data.Results.Display.Handles.Figures;
[data.Results.Display.Handles,OPTIONS] = ...
bst_wavedata_display(D,OPTIONS);
data.Results.Display.Handles.Figures = figs; clear figs
% <- DONE
% Now display original and modeled data + residuals
DisplayStruct = data.Results.Display;
dOPTIONS = struct(...
'ModalityLabel',{GoodTypes},...
'WaveLabel',{{'Original data','Modeled',sprintf('Residuals: %3.2f%%',...
100*mean(norcol(Frem))/mean(norcol(Fr)))}},...
'Channel', data.Measures.Channel(data.Measures.Display.SelectedChannels{data.Measures.Display.Modality}),...
'Time', Time, ...
'MoveCursor','on'...
);
dOPTIONS.SelectedChannels = cell(1,3);
dOPTIONS.SelectedChannels{data.Measures.Display.Modality} = chans;
% -> Fill out parameters before call to specific display subfunction
DefaultLabels = {'MEG','EEG','OTHER'}; % Default modality labels
figure(data.Measures.Display.Handles.Figures.TimeSeries), clf
figs = data.Measures.Display.Handles.Figures;
% was Time(GUI.Segment)`
dOPTIONS.DataType = 'surface';
[data.Measures.Display.Handles,dOPTIONS] = bst_wavedata_display({Fr,Frsyn,Frem},dOPTIONS);
data.Measures.Display.Flags = struct('SensorMarkers',OPTIONS.SensorMarkers,...
'SensorLabels',OPTIONS.SensorLabels, 'ShowContours', OPTIONS.ShowContours,...
'AbsoluteColormap',OPTIONS.AbsoluteColormap, 'Colorbar',OPTIONS.Colorbar);
data.Measures.Display.Handles.Figures = figs;
drawnow
if(0), % 10-May-2002 no more spatial display
if(hfdepthgauge), %is there a figure to draw to?
depth_point = best(Topo_i).L;
figure(hfdepthgauge)
hlinedepth = zeros(size(depth_point,2),1); % for each source location
for ihline = 1:length(hlinedepth),
[ignore,ignore,ignore,ignore,hlinedepth(ihline)] = ...
depthgauge(depth_point(:,ihline),...
depth_point(:,ihline)/norm(depth_point(:,ihline)),...
100/1000,5/1000,1000,gca);
end
% set the color of the depthgauge to that of the default line colors in the axes
tempcolor = get(gca,'ColorOrder');
mod_i = mod(Topo_i-1,size(tempcolor,1))+1; % modulo the length
set(hlinedepth,'Color',tempcolor(mod_i,:),'MarkerFaceColor',tempcolor(mod_i,:))
end
drawnow
pause(.25) % see if this helps force a screen image to appear
end % no more spatial display
% Here is the new spatial display (Nov-11 2003 - SB)
% JCM mods 13 Nov to account for missing tesselations
TessFile = bst_static_taskbar('GET','TESS');
existsTessFile = 0; % assume no
if ~isempty(TessFile) & exist(TessFile,'file')
existsTessFile = 1; % does exist
end
if existsTessFile,
if Topo_i == 1 % Fisrt source found %data.Results.Display.CorticalEnvelope
CorticEnvFigure = figure;
set(CorticEnvFigure,'visible','off' ); drawnow
setappdata(CorticEnvFigure,'TileType','T')
bst_layout('align',CorticEnvFigure,1,2,1);
bst_layout('align',data.Results.Display.Handles.Figures.TimeSeries,2,2,4);
bst_layout('align',data.Measures.Display.Handles.Figures.TimeSeries,2,2,2)
ImageGrid.GridLoc = {TessFile};
bst_message_window('wrap',{...
' ',...
sprintf('Loading %s tessellation file',TessFile)});
Users = get_user_directory;
load(fullfile(Users.SUBJECTS,TessFile),'Comment')
ImageGrid.iGrid = 1;
if length(Comment)>1 % Several surfaces in this tessellation file
[ImageGrid.iGrid,OK] = listdlg('PromptString','Select an envelope from list:',...
'SelectionMode','single',...
'ListString',Comment);
if OK == 0
return
end
end
envOPTIONS.Dipoles.SourceLoc = {best(1).L};
envOPTIONS.Dipoles.Orientation = {best(1).O};
envOPTIONS.Dipoles.TimeSeries = [D{:}]';
envOPTIONS.Time = 1;
envOPTIONS.currDensityTransparency = 1;
envOPTIONS.Dipoles.Color = data.Results.Display.Handles.LineColor;
[PlotHandles,envOPTIONS] = bst_imaging_display(...
CorticEnvFigure,ImageGrid,...
0,envOPTIONS);
envOPTIONS.hPatch = PlotHandles.hPatch;
else % Update display with next sources
envOPTIONS.Dipoles.SourceLoc = {best(1:Topo_i).L};
envOPTIONS.Dipoles.Orientation = {best(1:Topo_i).O};
envOPTIONS.Dipoles.TimeSeries = reshape([D{:}],length(D{1}),Topo_i);
envOPTIONS.Dipoles.Color = data.Results.Display.Handles.LineColor;
[PlotHandles,envOPTIONS] = bst_imaging_display(...
CorticEnvFigure,ImageGrid,...
0,envOPTIONS);
%envOPTIONS.hPatch = envPlotHandles.hPatch;
end % if first topo or subsequent topographies
else
CorticEnvFigure = []; % no tesselation figure
end % if a tesselation file exists
end % if user wanted graphics
% Note that the search grid will always be projected at the beginning of the loop
Topo_i = Topo_i + 1; % increment the counter to the next source topography
else % correlation too low, increment order
bst_message_window('append',{'',...
sprintf('Correlation too low at Order %.0f. Next Order',...
GUI.Order(nOrder))})
nOrder = nOrder + 1; % order indexer
if(nOrder <= length(GUI.Order)), % more orders left to search
% then load the next head model associated with the new order
% The Channel, GoodChannel, DataFlag don't change, just the HeadModel call
HeadModel = extract_channels(...
GUI.OrderHeadModel{nOrder},StudySubject,GUI); % next HeadModel
end
end % checking cross correlation value
end % while valid orders
if GUI.DisplayGraphics & existsTessFile,
data.Results.Display.Figure(2) = CorticEnvFigure;
data.Results.Display.OPTIONS = envOPTIONS;
end
% done with processing, now mop up.
best(Topo_i:end) = []; % remove last failed source and pre-allocated zeros
if(isempty(A)),% nothing, honey
S = [];
msgbox('No Sources Found','Notice','modal');
return
end
if(~isempty(UDP)),
Ap = A - UDP*(UDP'*A); % project model away
else
Ap = A; % no projector
end
% The data have already been projected away from the projector
S = [pinv(Ap)*Fr]'; % reduced rank time series, Hermitian transpose
% Now synthesize the data without the projector applied
Fsynth = Aorg * S'; % the synthesized data, in the unwhitened unprojected space
% now map to results
Comment = sprintf('RAP-MUSIC, rank %.0f, corr %.3f',GUI.Rank,GUI.Corr);
Date = datestr(datenum(now),0);
Subject = Subject;
Study = Study;
Time = [GUI.Segment];
ChannelFlag = GUI.ChannelFlag;
NoiseCov = [];
SourceCov = [];
Projector = Data.Projector;
SourceOrder = [best.Order];
SourceCorrelation = [best.c]; % the source correlations
SourceLoc = cell(1,length(best));
SourceOrientation = cell(1,length(best));
Aall = zeros(size(Ap));
for i = 1:length(best),
SourceLoc{i} = best(i).L; % map each solution into a cell.
SourceOrientation{i} = best(i).O;
end
ModelGain = [];
IndepTopo = Aorg; % Fixed orientation gain matrix
TimeSeries = S; % time series normalized by gain columns
PatchNdx = [];
PatchAmp = [];
ImageGridAmp = [];
ImageGridTime =[];
Function = mfilename; % name of this calling routine
Results = struct('Comment',Comment,'Function',mfilename,'StudySubject',StudySubject,...
'DataFlag',DataFlag,'GUI',GUI,'Date',Date,'Subject',Subject,'Study',Study,'Time',Time,...
'ChannelFlag',ChannelFlag,'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
'Projector',Projector,'SourceOrder',SourceOrder,'SourceCorrelation',SourceCorrelation,...
'SourceOrientation',[],'ModelGain',ModelGain,...
'IndepTopo',IndepTopo,'TimeSeries',TimeSeries,'PatchNdx',PatchNdx,...
'PatchAmp',PatchAmp,'ImageGridAmp',ImageGridAmp,'ImageGridTime',...
ImageGridTime,'Fsynth',Fsynth);
Results.SourceOrientation = SourceOrientation; % struct command would deal out this variable
Results.SourceLoc = SourceLoc; % ditto
if(isfield(GUI,'Results')), %user provided such a string
if(~isempty(GUI.Results)), % and it's not empty
save(fullfile(Users.STUDIES,GUI.Results),...
'Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
'DataFlag','ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder','SourceCorrelation',...
'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
bst_message_window('append',{'',sprintf('Wrote results to %s',GUI.Results)}); % to the command window
end
end
bst_message_window('append',{'',sprintf('RAP-MUSIC ALGORITHM DONE')}); % console screen
if GUI.DisplayGraphics,
figure(hTimeSeries); % bring to the front
drawnow
visfig = getappdata(hTimeSeries,'visfig');
if isempty(visfig) | ~ishandle(visfig), % not created or gone
visfig = results_visualization;
setappdata(hTimeSeries,'visfig',visfig); % save it to prevent duplication
end
vishandles = guidata(visfig);
results_visualization('loadparametric_Callback',...
vishandles.loadparametric, [], vishandles, Results.GUI.Results);
end
if GUI.FinalDisplayGraphics,
bst_mriviewer('LoadSourceMap_Callback','','','',Results.GUI.Results)
% visfig = results_visualization;
% vishandles = guidata(visfig);
% results_visualization('loadparametric_Callback',...
% vishandles.loadparametric, [], vishandles, Results.GUI.Results);
end
%---------------------------------------------------------------------
% snipped from nchoosek function in Matlab
function P = combs(v,m)
%COMBS All possible combinations.
% COMBS(1:N,M) or COMBS(V,M) where V is a row vector of length N,
% creates a matrix with N!/((N-M)! M!) rows and M columns containing
% all possible combinations of N elements taken M at a time.
%
% This function is only practical for situations where M is less
% than about 15.
if nargin~=2, error('Requires 2 input arguments.'); end
v = v(:).'; % Make sure v is a row vector.
n = length(v);
if n == m
P = v;
elseif m == 1
P = v.';
else
P = [];
if m < n & m > 1
for k = 1:n-m+1
Q = combs(v(k+1:n),m-1); % recursively calling
P = [P; [v(ones(size(Q,1),1),k) Q]];
end
end
end
% --------------- THE FMINS FUNCTION, USE WITH FUNCTION HANDLE --------------------
function e = min_fgain_gui(L,Function,Channel,Param,Order,GUI,Uf,Ua,Whitener,imegsens,irefsens,ieegsens);
%MIN_FGAIN_GUI minimize the error using the a specified gain model
% LIMITED: ONLY LOOKS AT THE FIRST FUNCTION CALL AND APPLIES EQUALL TO ALL CHANNELS!
% function e = min_fgain_gui(L,Function,Channel,Param,Order,GUI,Uf,Ua);
% call as fmins('min_fgain_gui',L,[],[],Function,Channel,Param,Order,GUI,Uf,Ua);
% Uf and Ua must already be orthogonal
% Uf is the subspace with which to correlate
% Ua is an existing subspace to project away from
% Function,Channel,Param,Order
% Function is the name of the gain matrix call, e.g. "os_meg"
% John C. Mosher, Ph.D.
% See Copyright.m file for information
% $Date: 6/14/04 3:38p $ $Revision: 44 $
% generate the gain matrix
global GoodChannelMEG GoodChannelEEG
%%%%% CHEAT %%%%%%%%
% Fix is to look for all calls to the same function and call in clusters.
if(0) % too expensive too call one at a time.
DIMS = [3 3 12];
Dims= DIMS(Order+2);
G = zeros(length(Channel),Dims);
for i =1:length(Channel),
G(i,:) = feval(Function{i},L,Channel(i),Param(i),Order);
end
else
% CHEAT: CALL ONLY THE FIRST FUNCTION
global ChannelFlag
ChannelFlag = GUI.ChannelFlag; % Need info on bad channels for CTF gradient computation.
G= feval(Function{1},L,Channel,Param,Order,imegsens,irefsens,ieegsens);
end
if ~isempty(Whitener),
G = Whitener*G; % data were already whitened
end
if(0), % disabled by JCM is irrelevant
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following is relevant for BEM models only
if ~isempty(findstr(lower(Function{1}),'meg'))
if length(GoodChannelMEG) < size(G,1) % There are bad MEG Channels
G = G(GoodChannelMEG-GoodChannelMEG(1)+1,:);
end
else
if ~isempty(findstr(lower(Function{1}),'eeg'))
if length(GoodChannelEEG) < size(G,1) % There are bad EEG Channels
G = G(GoodChannelEEG-GoodChannelEEG(1)+1,:);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
if(GUI.Column_norm), % scale the model
[ignore,G] = colnorm(G);
end
% if the user provided previous recursions
% project away from it
if(exist('Ua','var')),
if(~isempty(Ua))
G = G - Ua*(Ua'*G);
end
end
% since this function is repeatedly called in directed minimization, reluctant
% to call subcorr_gui here. More efficient here because we know that Uf is
% already orthogonal.
% Orthogonalize
[Ug,Sg,Vg] = svd(G,0);
% now set Ug to the reduced rank called for by the regularization
Ug = regsubspace(GUI,Ug,Sg);
% subspace correlations
c = svd(Ug'*Uf);
if(isempty(c)), % nothing valid
e = 0;
else
e = -c(1); % negative maximum for minimizing
end
Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON