|
Dear Sue,
Here is the current version of the code. This is the internal version which hasn't changed since sept last year so it should be in any new downloads.
As to your questions.
I hadn't noticed the theta phi before and will look at it, but all the information you need is in the moment and orientation.
To answer your question about how many dipoles to use, you should start with 1 and keep increasing until the model evidence goes down.
Similarly if you have prior knowledge you can test models consisting of different numbers of dipoles in different locations.
The paper you need to look at is
Variational Bayesian inversion of the equivalent current dipole model in EEG/MEG.
Kiebel SJ, Daunizeau J, Phillips C, Friston KJ.
Neuroimage. 2008 Jan 15;39(2):728-41. Epub 2007 Sep 14.
Best wishes
Gareth
-----Original Message-----
From: Sun Delin [mailto:[log in to unmask]]
Sent: 01 February 2011 13:34
To: Gareth Barnes
Subject: Re: RE: [SPM] Error occurred when doing dipole analysis by using SPM8
Dear Gareth Barnes,
Great thanks for your updated script. I tried it on my ERP data by reversing one of four conditions within the time window of 50 ms length (in fact, 550~600 ms after face presentation) with four single dipoles without prior information of source location or moment. No error occurred. However, the value of "theta-phi or" in the "Dipole orientation & strength" seems not being presented correctly (see attached output). Moreover, I would like to know how many dipoles, and single or symatic dipoles I should employ to run the dipole analysis? What parameters could be used to judge which model is the best? You could see that the four dipoles in my output file were all located deep inside the brain. This is far from my hypothesis that the sources of the LPC component (550~600 ms) should be located at several different brain regions.
Best regards,
Dr. Sun Delin
Post-doctoral Fellow
Laboratory of Neuropsychology and Laboratory of Cognitive Affective Neuroscience, The University of Hong Kong Room 622B, Knowles Building, The University of Hong Kong, Pokfulam Road, Hong Kong Tel. (Mobile) : (852) 5174 1885 Tel. (Office) : (852) 2241 5655
URL: http://www.researcherid.com/rid/A-4154-2010
http://hub.hku.hk/rp/rp00873
email: [log in to unmask]
2011-02-01
======= At 2011-02-01, 19:44:15 you wrote: =======
>Dear Sue
>Please try the current internal version code and see if it solves the
>problem. If it does I'll send to the list.
>Best wishes
>Gareth
>
>
>-----Original Message-----
>From: Vladimir Litvak [mailto:[log in to unmask]]
>Sent: 01 February 2011 11:27
>To: Gareth Barnes
>Subject: Fwd: [SPM] Error occurred when doing dipole analysis by using
>SPM8
>
>---------- Forwarded message ----------
>From: Sun Delin <[log in to unmask]>
>Date: Tue, Feb 1, 2011 at 5:26 AM
>Subject: [SPM] Error occurred when doing dipole analysis by using SPM8
>To: [log in to unmask]
>
>
>Dear SPMers,
>
>? met an error when doing dipole analysis by SPM8 (version 4010) as
>follows:
>??? Error using ==> mtimes
>Inner matrix dimensions must agree.
>
>Error in ==> spm_eeg_inv_vbecd_gui at 478 ?nverse.jmni{ii} =
>orM1*P.post_mu_w; % dipole(s) orient/ampl in mni space
>
>Error in ==> spm_eeg_invert_ui at 51
>??? = spm_eeg_inv_vbecd_gui(D,val);
>
>Error in ==> spm_eeg_inv_imag_api>Inverse_Callback at 94 handles.D =
>spm_eeg_invert_ui(handles.D);
>
>Error in ==> spm_eeg_inv_imag_api at 53 ???eval(varargin{:}); % FEVAL
>switchyard
>
>??? Error while evaluating uicontrol Callback
>
>?his error occurred after 10 iterations of inversion. However, such
>dipole inversion ran smoothly when I use an older version of SPM8, i.e.
>version 3648. Has anyone met such problem?
>
>Bests,
>Sun Delin
= = = = = = = = = = = = = = = = = = = =
function D = spm_eeg_inv_vbecd_gui(D,val)
% GUI function for Bayesian ECD inversion
% - load the necessary data, if not provided
% - fill in all the necessary bits for the VB-ECD inversion routine,
% - launch the B_ECD routine, aka. spm_eeg_inv_vbecd
% - displays the results.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
%
% $Id: spm_eeg_inv_vbecd_gui.m 4071 2010-09-22 13:44:04Z gareth $
%%
% Load data, if necessary
%==========
if nargin<1
D = spm_eeg_load;
end
%%
% Check if the forward model was prepared & handle the other info bits
%========================================
if ~isfield(D,'inv')
error('Data must have been prepared for inversion procedure...')
end
if nargin==2
% check index provided
if val>length(D.inv)
val = length(D.inv);
D.val = val;
end
else
if isfield(D,'val')
val = D.val;
else
% use last one
val = length(D.inv);
D.val = val;
end
end
% Use val to define which is the "current" inv{} to use
% If no inverse solution already calculated (field 'inverse' doesn't exist)
% use that inv{}. Otherwise create a new one by copying the previous
% inv{} structure
if isfield(D.inv{val},'inverse')
% create an extra inv{}
Ninv = length(D.inv);
D.inv{Ninv+1} = D.inv{val};
if isfield(D.inv{Ninv+1},'contrast')
% no contrast field used here !
D.inv{Ninv+1} = rmfield(D.inv{Ninv+1},'contrast');
end
val = Ninv+1;
D.val = val;
end
if ~isfield(D.inv{val}, 'date')
% Set time , date, comments & modality
clck = fix(clock);
if clck(5) < 10
clck = [num2str(clck(4)) ':0' num2str(clck(5))];
else
clck = [num2str(clck(4)) ':' num2str(clck(5))];
end
D.inv{val}.date = strvcat(date,clck); %#ok<VCAT>
end
if ~isfield(D.inv{val}, 'comment'),
D.inv{val}.comment = {spm_input('Comment/Label for this analysis', '+1', 's')};
end
D.inv{val}.method = 'vbecd';
%% Struct that collects the inputs for vbecd code
P = [];
P.modality = spm_eeg_modality_ui(D, 1, 1);
if isfield(D.inv{val}, 'forward') && isfield(D.inv{val}, 'datareg')
for m = 1:numel(D.inv{val}.forward)
if strncmp(P.modality, D.inv{val}.forward(m).modality, 3)
P.forward.vol = D.inv{val}.forward(m).vol;
if ischar(P.forward.vol)
P.forward.vol = ft_read_vol(P.forward.vol);
end
P.forward.sens = D.inv{val}.datareg(m).sensors;
% Channels to use
P.Ic = setdiff(meegchannels(D, P.modality), badchannels(D));
M1 = D.inv{val}.datareg.toMNI;
[U, L, V] = svd(M1(1:3, 1:3));
orM1(1:3,1:3) =U*V'; %% for switching orientation between meg and mni space
% disp('Undoing transformation to Tal space !');
%disp('Fixing sphere centre !');
%P.forward.vol.o=[0 0 28]; P.forward.vol.r=100;
mnivol = ft_transform_vol(M1, P.forward.vol); %% used for inside head calculation
end
end
end
if isempty(P.Ic)
error(['The specified modality (' P.modality ') is missing from file ' D.fname]);
else
P.channels = D.chanlabels(P.Ic);
end
[P.forward.vol, P.forward.sens] = ft_prepare_vol_sens( ...
P.forward.vol, P.forward.sens, 'channel', P.channels);
if ~isfield(P.forward.sens,'prj')
P.forward.sens.prj = D.coor2D(P.Ic);
end
%%
% Deal with data
%===============
% time bin or time window
msg_tb = ['time_bin or average_win [',num2str(round(min(D.time)*1e3)), ...
' ',num2str(round(max(D.time)*1e3)),'] ms'];
ask_tb = 1;
while ask_tb
tb = spm_input(msg_tb,1,'r'); % ! in msec
if length(tb)==1
if tb>=min(D.time([], 'ms')) && tb<=max(D.time([], 'ms'))
ask_tb = 0;
end
elseif length(tb)==2
if all(tb>=floor(min(D.time([], 'ms')))) && all(tb<=ceil(max(D.time([], 'ms')))) && tb(1)<=tb(2)
ask_tb = 0;
end
end
end
if length(tb)==1
[kk,ltb] = min(abs(D.time([], 'ms')-tb)); % round to nearest time bin
else
[kk,ltb(1)] = min(abs(D.time([], 'ms')-tb(1))); % round to nearest time bin
[kk,ltb(2)] = min(abs(D.time([], 'ms')-tb(2)));
ltb = ltb(1):ltb(2); % list of time bins 'tb' to use
end
% trial type
if D.ntrials>1
msg_tr = ['Trial type number [1 ',num2str(D.ntrials),']'];
ltr = spm_input(msg_tr,2,'i',num2str(1:D.ntrials));
tr_q = 1;
else
tr_q = 0;
ltr = 1;
end
% data, averaged over time window considered
EEGscale=1;
%% SORT OUT EEG UNITS AND CONVERT VALUES TO VOLTS
if strcmp(upper(P.modality),'EEG'),
allunits=strvcat('uV','mV','V');
allscales=[1e-6, 1e-3, 1]; %%
EEGscale=0;
eegunits = unique(D.units(D.meegchannels('EEG')));
Neegchans=numel(D.units(D.meegchannels('EEG')));
for j=1:length(allunits),
if strcmp(deblank(allunits(j,:)),deblank(eegunits));
EEGscale=allscales(j);
end; % if
end; % for j
if EEGscale==0,
warning('units unspecified');
if mean(std(D(P.Ic,ltb,ltr)))>1e-2,
guess_ind=[1 2 3];
else
guess_ind=[3 2 1];
end;
msg_str=sprintf('Units of EEG are %s ? (rms=%3.2e)',allunits(guess_ind(1),:),mean(std(D(P.Ic,ltb,ltr))));
dip_ch = sprintf('%s|%s|%s',allunits(guess_ind(1),:),allunits(guess_ind(2),:),allunits(guess_ind(3),:));
dip_val = [1,2,3];
def_opt=1;
unitind= spm_input(msg_str,2,'b',dip_ch,dip_val,def_opt);
%ans=spm_input(msg_str,1,'s','yes');
allunits(guess_ind(unitind),:)
D = units(D, 1:Neegchans, allunits(guess_ind(unitind),:));
EEGscale=allscales(guess_ind(unitind));
D.save; %% Save the new units
end; %if EEGscale==0
end; % if eeg data
dat_y = squeeze(mean(D(P.Ic,ltb,ltr)*EEGscale,2));
%%
% Other bits of the P structure, apart for priors and #dipoles
%==============================
P.ltr = ltr;
P.Nc = length(P.Ic);
% Deal with dipoles number and priors
%====================================
dip_q = 0; % number of dipole 'elements' added (single or pair)
dip_c = 0; % total number of dipoles in the model
adding_dips = 1;
clear dip_pr
priorlocvardefault=[100, 100, 100]; %% location variance default in mm
nopriorlocvardefault=[80*80, 80*80, 80*80];
nopriormomvardefault=[10, 10, 10]*100; %% moment variance in nAM
priormomvardefault=[1, 1, 1]; %%
while adding_dips
if dip_q>0,
msg_dip =['Add dipoles to ',num2str(dip_c),' or stop?'];
dip_ch = 'Single|Symmetric Pair|Stop';
dip_val = [1,2,0];
def_opt=3;
else
msg_dip =['Add dipoles to model'];
def_opt=1;
dip_ch = 'Single|Symmetric Pair';
dip_val = [1,2];
end
a_dip = spm_input(msg_dip,2+tr_q+dip_q,'b',dip_ch,dip_val,def_opt);
if a_dip == 0
adding_dips = 0;
elseif a_dip == 1
% add a single dipole to the model
dip_q = dip_q+1;
dip_pr(dip_q) = struct( 'a_dip',a_dip, ...
'mu_w0',[],'mu_s0',[],'S_s0',eye(3),'S_w0',eye(3));
%% 'ab20',[],'ab30',[]); %% ,'Tw', [],'Ts', []);
% Location prior
spr_q = spm_input('Location prior ?',1+tr_q+dip_q+1,'b', ...
'Informative|Non-info',[1,0],2);
if spr_q
% informative location prior
str = 'Location prior';
while 1
s0mni = spm_input(str, 1+tr_q+dip_q+2,'e',[0 0 0])';
outside = ~ft_inside_vol(s0mni',mnivol);
s0=D.inv{val}.datareg.fromMNI*[s0mni' 1]';
s0=s0(1:3);
str2='Prior location variance (mm2)';
diags_s0_mni = spm_input(str2, 1+tr_q+dip_q+2,'e',priorlocvardefault)';
S_s0_ctf=orM1*diag(diags_s0_mni)*orM1'; %% transform covariance
%% need to leave diags(S0) free
if all(~outside), break, end
str = 'Prior location must be inside head';
end
dip_pr(dip_q).mu_s0 = s0;
else
% no location prior
dip_pr(dip_q).mu_s0 = zeros(3,1);
diags_s0_mni= nopriorlocvardefault';
S_s0_ctf=diag(diags_s0_mni);
end
dip_pr(dip_q).S_s0=S_s0_ctf; %
% Moment prior
wpr_q = spm_input('Moment prior ?',1+tr_q+dip_q+spr_q+2,'b', ...
'Informative|Non-info',[1,0],2);
if wpr_q
% informative moment prior
w0_mni= spm_input('Moment prior', ...
1+tr_q+dip_q+spr_q+3,'e',[0 0 0])';
str2='Prior moment variance (nAm2)';
diags_w0_mni = spm_input(str2, 1+tr_q+dip_q+2,'e',priormomvardefault)';
dip_pr(dip_q).mu_w0 =orM1*w0_mni;
S_w0_ctf=orM1*w0_mni*orM1';
else
% no location prior
dip_pr(dip_q).mu_w0 = zeros(3,1);
S_w0_ctf= diag(nopriormomvardefault);
end
%% set up covariance matrix for orientation with no crosstalk terms (for single
%% dip)
dip_pr(dip_q).S_w0=S_w0_ctf;
dip_c = dip_c+1;
else
% add a pair of symmetric dipoles to the model
dip_q = dip_q+1;
dip_pr(dip_q) = struct( 'a_dip',a_dip, ...
'mu_w0',[],'mu_s0',[],'S_s0',eye(6),'S_w0',eye(6));
%%...
% 'ab20',[],'ab30',[]); %%,'Tw',eye(6),'Ts',eye(6));
% Location prior
spr_q = spm_input('Location prior ?',1+tr_q+dip_q+1,'b', ...
'Informative|Non-info',[1,0],2);
if spr_q
% informative location prior
str = 'Location prior (one side only)';
while 1
s0mni = spm_input(str, 1+tr_q+dip_q+2,'e',[0 0 0])';
syms0mni=s0mni;
syms0mni(1)=-syms0mni(1);
outside = ~ft_inside_vol(s0mni',mnivol);
s0=D.inv{val}.datareg.fromMNI*[s0mni' 1]';
s0sym=D.inv{val}.datareg.fromMNI*[syms0mni' 1]';
str2='Prior location variance (mm2)';
tmp_diags_s0_mni = spm_input(str2, 1+tr_q+dip_q+2,'e',priorlocvardefault)';
tmp_diags_s0_mni= [tmp_diags_s0_mni ; tmp_diags_s0_mni ];
if all(~outside), break, end
str = 'Prior location must be inside head';
end
dip_pr(dip_q).mu_s0 = [s0(1:3);s0sym(1:3)];
else
% no location prior
dip_pr(dip_q).mu_s0 = zeros(6,1);
tmp_diags_s0 = [nopriorlocvardefault';nopriorlocvardefault'];
tmp_diags_s0_mni=[nopriorlocvardefault';nopriorlocvardefault'];
end %% end of if informative prior
%% setting up a covariance matrix where there is covariance between
%% the x parameters negatively coupled, y,z positively.
mni_dip_pr(dip_q).S_s0 = eye(length(tmp_diags_s0_mni)).*repmat(tmp_diags_s0_mni,1,length(tmp_diags_s0_mni));
mni_dip_pr(dip_q).S_s0(4,1)=-mni_dip_pr(dip_q).S_s0(4,4); % reflect in x
mni_dip_pr(dip_q).S_s0(5,2)=mni_dip_pr(dip_q).S_s0(5,5); % maintain y and z
mni_dip_pr(dip_q).S_s0(6,3)=mni_dip_pr(dip_q).S_s0(6,6);
mni_dip_pr(dip_q).S_s0(1,4)=mni_dip_pr(dip_q).S_s0(4,1);
mni_dip_pr(dip_q).S_s0(2,5)=mni_dip_pr(dip_q).S_s0(5,2);
mni_dip_pr(dip_q).S_s0(3,6)=mni_dip_pr(dip_q).S_s0(6,3);
%% transform to MEG space
%dip_pr(dip_q).S_s0(:,1:3)=orM1*mni_dip_pr(dip_q).S_s0(:,1:3)*orM1'; %% NEED TO LOOK AT THIS
%dip_pr(dip_q).S_s0(:,4:6)=orM1*mni_dip_pr(dip_q).S_s0(:,4:6)*orM1';
tmp1=orM1*mni_dip_pr(dip_q).S_s0(1:3,1:3)*orM1'; %% NEED TO LOOK AT THIS
tmp2=orM1*mni_dip_pr(dip_q).S_s0(1:3,4:6)*orM1';
tmp3=orM1*mni_dip_pr(dip_q).S_s0(4:6,4:6)*orM1';
tmp4=orM1*mni_dip_pr(dip_q).S_s0(4:6,1:3)*orM1';
dip_pr(dip_q).S_s0(1:3,1:3)=tmp1;
dip_pr(dip_q).S_s0(1:3,4:6)=tmp2;
dip_pr(dip_q).S_s0(4:6,4:6)=tmp3;
dip_pr(dip_q).S_s0(4:6,1:3)=tmp4;
% Moment prior
wpr_q = spm_input('Moment prior ?',1+tr_q+dip_q+spr_q+2,'b', ...
'Informative|Non-info',[1,0],2);
if wpr_q
% informative moment prior
tmp= spm_input('Moment prior (right only)', ...
1+tr_q+dip_q+spr_q+3,'e',[1 1 1])';
tmp = [tmp ; tmp] ; tmp(4) = tmp(4);
dip_pr(dip_q).mu_w0 = tmp;
str2='Prior moment variance (nAm2)';
diags_w0 = spm_input(str2, 1+tr_q+dip_q+spr_q+3,'e',priormomvardefault)';
tmp_diags_w0=[diags_w0; diags_w0];
else
% no moment prior
dip_pr(dip_q).mu_w0 = zeros(6,1);
tmp_diags_w0 = [nopriormomvardefault'; nopriormomvardefault'];
end
%dip_pr(dip_q).S_w0=eye(length(diags_w0)).*repmat(diags_w0,1,length(diags_w0));
%% couple all orientations, except x, positively or leave for now...
mni_dip_pr(dip_q).S_w0 = eye(length(tmp_diags_w0)).*repmat(tmp_diags_w0,1,length(tmp_diags_w0));
mni_dip_pr(dip_q).S_w0(4,1)=-mni_dip_pr(dip_q).S_w0(4,4); % reflect x orientation
mni_dip_pr(dip_q).S_w0(5,2)=mni_dip_pr(dip_q).S_w0(5,5); %
mni_dip_pr(dip_q).S_w0(6,3)=mni_dip_pr(dip_q).S_w0(6,6); %
mni_dip_pr(dip_q).S_w0(1,4)=-mni_dip_pr(dip_q).S_w0(4,1); %
mni_dip_pr(dip_q).S_w0(2,5)=mni_dip_pr(dip_q).S_w0(5,2); %
mni_dip_pr(dip_q).S_w0(3,6)=mni_dip_pr(dip_q).S_w0(6,3); %
tmp1=orM1*mni_dip_pr(dip_q).S_w0(1:3,1:3)*orM1';
tmp2=orM1*mni_dip_pr(dip_q).S_w0(1:3,4:6)*orM1';
tmp3=orM1*mni_dip_pr(dip_q).S_w0(4:6,4:6)*orM1';
tmp4=orM1*mni_dip_pr(dip_q).S_w0(4:6,1:3)*orM1';
dip_pr(dip_q).S_w0(1:3,1:3)=tmp1;
dip_pr(dip_q).S_w0(1:3,4:6)=tmp2;
dip_pr(dip_q).S_w0(4:6,4:6)=tmp3;
dip_pr(dip_q).S_w0(4:6,1:3)=tmp4;
dip_c = dip_c+2;
end
end
%str2='Data SNR (amp)';
% SNRamp = spm_input(str2, 1+tr_q+dip_q+2+1,'e',5)';
SNRamp=3;
hE=log(SNRamp^2); %% expected log precision of data
hC=1; % variability of the above precision
str2='Number of iterations';
Niter = spm_input(str2, 1+tr_q+dip_q+2+2,'e',10)';
%%
% Get all the priors together and build structure to pass to inv_becd
%============================
priors = struct('mu_w0',cat(1,dip_pr(:).mu_w0), ...
'mu_s0',cat(1,dip_pr(:).mu_s0), ...
'S_w0',blkdiag(dip_pr(:).S_w0),'S_s0',blkdiag(dip_pr(:).S_s0),'hE',hE,'hC',hC);
P.priors = priors;
%%
% Launch inversion !
%===================
% Initialise inverse field
inverse = struct( ...
'F',[], ... % free energy
'pst',D.time, ... % all time points in data epoch
'tb',tb, ... % time window/bin used
'ltb',ltb, ... % list of time points used
'ltr',ltr, ... % list of trial types used
'n_seeds',length(ltr), ... % using this field for multiple reconstruction
'n_dip',dip_c, ... % number of dipoles used
'loc',[], ... % loc of dip (3 x n_dip)
'j',[], ... % dipole(s) orient/ampl, in 1 column
'cov_loc',[], ... % cov matrix of source location
'cov_j',[], ... % cov matrix of source orient/ampl
'Mtb',1, ... % ind of max EEG power in time series, 1 as only 1 tb.
'exitflag',[], ... % Converged (1) or not (0)
'P',[]); % save all kaboodle too.
for ii=1:length(ltr)
P.y = dat_y(:,ii);
P.ii = ii;
%% set up figures
P.handles.hfig = spm_figure('GetWin','Graphics');
spm_clf(P.handles.hfig)
P.handles.SPMdefaults.col = get(P.handles.hfig,'colormap');
P.handles.SPMdefaults.renderer = get(P.handles.hfig,'renderer');
set(P.handles.hfig,'userdata',P)
dip_amp=[];
for j=1:Niter,
Pout(j) = spm_eeg_inv_vbecd(P);
close(gcf);
varresids(j)=var(Pout(j).y-Pout(j).ypost);
pov(j)=100*(1-varresids(j)/var(Pout(j).y)); %% percent variance explained
allF(j)=Pout(j).F;
dip_mom=reshape(Pout(j).post_mu_w,3,length(Pout(j).post_mu_w)/3);
dip_amp(j,:)=sqrt(dot(dip_mom,dip_mom));
% display
megloc=reshape(Pout(j).post_mu_s,3,length(Pout(j).post_mu_s)/3); % loc of dip (3 x n_dip)
mniloc=D.inv{val}.datareg.toMNI*[megloc;ones(1,size(megloc,2))]; %% actual MNI location (with scaling)
megmom=reshape(Pout(j).post_mu_w,3,length(Pout(j).post_mu_w)/3); % moments of dip (3 x n_dip)
megposvar=reshape(diag(Pout(j).post_S_s),3,length(Pout(j).post_mu_s)/3); %% estimate of positional uncertainty in three principal axes
mnimom=orM1*megmom; %% convert moments into mni coordinates through a rotation (no scaling or translation)
mniposvar=(orM1*sqrt(megposvar)).^2; %% convert pos variance into approx mni space by switching axes
displayVBupdate2(Pout(j).y,pov,allF,Niter,dip_amp,mnimom,mniloc(1:3,:),mniposvar,P,j,[],Pout(j).F,Pout(j).ypost,[]);
end; % for j
allF=[Pout.F];
[maxFvals,maxind]=max(allF);
P=Pout(maxind); %% take best F
% Get the results out.
inverse.pst = tb*1e3;
inverse.F(ii) = P.F; % free energy
megloc=reshape(P.post_mu_s,3,length(P.post_mu_s)/3); % loc of dip (3 x n_dip)
meg_w=reshape(P.post_mu_w,3,length(P.post_mu_w)/3); % moments of dip (3 x n_dip)
mni_w=orM1*meg_w; %% orientation in mni space
mniloc=D.inv{val}.datareg.toMNI*[megloc;ones(1,size(megloc,2))]; %% actual MNI location (with scaling)
inverse.mniloc{ii}=mniloc(1:3,:);
inverse.loc{ii} = megloc;
inverse.j{ii} = P.post_mu_w; % dipole(s) orient/ampl, in 1 column in meg space
inverse.jmni{ii} = reshape(mni_w,1,prod(size(mni_w)))'; % dipole(s) orient/ampl in mni space
inverse.cov_loc{ii} = P.post_S_s; % cov matrix of source location
inverse.cov_j{ii} = P.post_S_w; % cov matrix of source orient/ampl
inverse.exitflag(ii) = 1; % Converged (1) or not (0)
inverse.P{ii} = P; % save all kaboodle too.
%% show final result
pause(1);
spm_clf(P.handles.hfig)
megmom=reshape(Pout(maxind).post_mu_w,3,length(Pout(maxind).post_mu_w)/3); % moments of dip (3 x n_dip)
%megposvar=reshape(diag(Pout(maxind).post_S_s),3,length(Pout(maxind).post_mu_s)/3); %% estimate of positional uncertainty in three principal axes
mnimom=orM1*megmom; %% convert moments into mni coordinates through a rotation (no scaling or translation)
longorM1=zeros(size(Pout(maxind).post_S_s,1));
for k1=1:length(Pout(maxind).post_S_s)/3;
longorM1((k1-1)*3+1:k1*3,(k1-1)*3+1:k1*3)=orM1;
end; % for k1
S0_mni=longorM1*Pout(maxind).post_S_s*longorM1';
mniposvar=diag(S0_mni); %% convert pos variance into approx mni space by switching axes
mniposvar=reshape(mniposvar,3,length(Pout(maxind).post_S_s)/3);
displayVBupdate2(Pout(j).y,pov,allF,Niter,dip_amp,mnimom,mniloc(1:3,:),mniposvar,P,j,[],Pout(j).F,Pout(j).ypost,maxind);
%displayVBupdate2(Pout(maxind).y,pov,allF,Niter,dip_amp,mniloc,Pout(maxind).post_mu_s,Pout(maxind).post_S_s,P,j,[],Pout(maxind).F,Pout(maxind).ypost,maxind,D);
%
end
D.inv{val}.inverse = inverse;
%%
% Save results and display
%-------------------------
save(D)
return
function [P] = displayVBupdate2(y,pov_iter,F_iter,maxit,dipamp_iter,mu_w,mu_s,diagS_s,P,it,flag,F,yHat,maxind)
%% yHat is estimate of y based on dipole position
if ~exist('flag','var')
flag = [];
end
if ~exist('maxind','var')
maxind = [];
end
if isempty(flag) || isequal(flag,'ecd')
% plot dipoles
try
opt.ParentAxes = P.handles.axesECD;
opt.hfig = P.handles.hfig;
opt.handles.hp = P.handles.hp;
opt.handles.hq = P.handles.hq;
opt.handles.hs = P.handles.hs;
opt.handles.ht = P.handles.ht;
opt.query = 'replace';
catch
P.handles.axesECD = axes(...
'parent',P.handles.hfig,...
'Position',[0.13 0.55 0.775 0.4],...
'hittest','off',...
'visible','off',...
'deleteFcn',@back2defaults);
opt.ParentAxes = P.handles.axesECD;
opt.hfig = P.handles.hfig;
end
w = reshape(mu_w,3,[]);
s = reshape(mu_s, 3, []);
% mesh.faces=D.inv{D.val}.forward.mesh.face; %% in ctf space
% mesh.vertices=D.inv{D.val}.forward.mesh.vert; %% in ctf space
% [out] = spm_eeg_displayECD_ctf(...
% s,w,reshape(diag(S_s),3,[]),[],mesh,opt);
[out] = spm_eeg_displayECD(...
s,w,diagS_s,[],opt);
P.handles.hp = out.handles.hp;
P.handles.hq = out.handles.hq;
P.handles.hs = out.handles.hs;
P.handles.ht = out.handles.ht;
end
% plot data and predicted data
pos = P.forward.sens.prj;
ChanLabel = P.channels;
in.f = P.handles.hfig;
in.noButtons = 1;
try
P.handles.axesY;
catch
figure(P.handles.hfig)
P.handles.axesY = axes(...
'Position',[0.02 0.3 0.3 0.2],...
'hittest','off');
in.ParentAxes = P.handles.axesY;
spm_eeg_plotScalpData(y,pos,ChanLabel,in);
title(P.handles.axesY,'measured data')
end
if isempty(flag) || isequal(flag,'data') || isequal(flag,'ecd')
%yHat = P.gmn*mu_w;
miY = min([yHat;y]);
maY = max([yHat;y]);
try
P.handles.axesYhat;
d = get(P.handles.axesYhat,'userdata');
yHat = yHat(d.goodChannels);
clim = [min(yHat(:))-( max(yHat(:))-min(yHat(:)) )/63,...
max(yHat(:))];
ZI = griddata(...
d.interp.pos(1,:),d.interp.pos(2,:),full(double(yHat)),...
d.interp.XI,d.interp.YI);
set(d.hi,'Cdata',flipud(ZI));
caxis(P.handles.axesYhat,clim);
delete(d.hc)
[C,d.hc] = contour(P.handles.axesYhat,flipud(ZI),...
'linecolor',0.5.*ones(3,1));
set(P.handles.axesYhat,...
'userdata',d);
catch
figure(P.handles.hfig)
P.handles.axesYhat = axes(...
'Position',[0.37 0.3 0.3 0.2],...
'hittest','off');
in.ParentAxes = P.handles.axesYhat;
spm_eeg_plotScalpData(yHat,pos,ChanLabel,in);
title(P.handles.axesYhat,'predicted data')
end
try
P.handles.axesYhatY;
catch
figure(P.handles.hfig)
P.handles.axesYhatY = axes(...
'Position',[0.72 0.3 0.25 0.2],...
'NextPlot','replace',...
'box','on');
end
plot(P.handles.axesYhatY,y,yHat,'.')
set(P.handles.axesYhatY,...
'nextplot','add')
plot(P.handles.axesYhatY,[miY;maY],[miY;maY],'r')
set(P.handles.axesYhatY,...
'nextplot','replace')
title(P.handles.axesYhatY,'predicted vs measured data')
axis(P.handles.axesYhatY,'square','tight')
grid(P.handles.axesYhatY,'on')
end
if isempty(flag) || isequal(flag,'var')
% plot precision hyperparameters
try
P.handles.axesVar1;
catch
figure(P.handles.hfig)
P.handles.axesVar1 = axes(...
'Position',[0.05 0.05 0.25 0.2],...
'NextPlot','replace',...
'box','on');
end
plot(P.handles.axesVar1,F_iter,'o-');
if ~isempty(maxind),
hold on;
h=plot(P.handles.axesVar1,maxind,F_iter(maxind),'rd');
set(h,'linewidth',4);
end;
set(P.handles.axesVar1,'Xlimmode','manual');
set(P.handles.axesVar1,'Xlim',[1 maxit]);
set(P.handles.axesVar1,'Xtick',1:maxit);
set(P.handles.axesVar1,'Xticklabel',num2str([1:maxit]'));
set(P.handles.axesVar1,'Yticklabel','');
title(P.handles.axesVar1,'Free energy ')
axis(P.handles.axesVar1,'square');
set(P.handles.axesVar1,'Ylimmode','auto'); %,'tight')
grid(P.handles.axesVar1,'on')
try
P.handles.axesVar2;
catch
figure(P.handles.hfig)
P.handles.axesVar2 = axes(...
'Position',[0.37 0.05 0.25 0.2],...
'NextPlot','replace',...
'box','on');
end
plot(P.handles.axesVar2,pov_iter,'*-')
if ~isempty(maxind),
hold on;
h=plot(P.handles.axesVar2,maxind,pov_iter(maxind),'rd');
set(h,'linewidth',4);
end;
set(P.handles.axesVar2,'Xlimmode','manual');
set(P.handles.axesVar2,'Xlim',[1 maxit]);
set(P.handles.axesVar2,'Xtick',1:maxit);
set(P.handles.axesVar2,'Xticklabel',num2str([1:maxit]'));
set(P.handles.axesVar2,'Ylimmode','manual'); %,'tight')
set(P.handles.axesVar2,'Ylim',[0 100]);
set(P.handles.axesVar2,'Ytick',[0:20:100]);
set(P.handles.axesVar2,'Yticklabel',num2str([0:20:100]'));
%set(P.handles.axesVar2,'Yticklabel','');
title(P.handles.axesVar2,'Percent variance explained');
axis(P.handles.axesVar2,'square');
grid(P.handles.axesVar2,'on')
try
P.handles.axesVar3;
catch
figure(P.handles.hfig)
P.handles.axesVar3 = axes(...
'Position',[0.72 0.05 0.25 0.2],...
'NextPlot','replace',...
'box','on');
end
plot(P.handles.axesVar3,1:it,dipamp_iter','o-');
if ~isempty(maxind),
hold on;
h=plot(P.handles.axesVar3,maxind,dipamp_iter(maxind,:)','rd');
set(h,'linewidth',4);
end;
set(P.handles.axesVar3,'Xlimmode','manual');
set(P.handles.axesVar3,'Xlim',[1 maxit]);
set(P.handles.axesVar3,'Xtick',1:maxit);
set(P.handles.axesVar3,'Xticklabel',num2str([1:maxit]'));
set(P.handles.axesVar3,'Yticklabel','');
title(P.handles.axesVar3,'Dipole amp (nAm) ')
axis(P.handles.axesVar3,'square');
set(P.handles.axesVar3,'Ylimmode','auto'); %,'tight')
grid(P.handles.axesVar3,'on')
end
if ~isempty(flag) && (isequal(flag,'ecd') || isequal(flag,'mGN') )
try
P.handles.hte(2);
catch
figure(P.handles.hfig)
P.handles.hte(2) = uicontrol('style','text',...
'units','normalized',...
'position',[0.2,0.91,0.6,0.02],...
'backgroundcolor',[1,1,1]);
end
set(P.handles.hte(2),'string',...
['ECD locations: Modified Gauss-Newton scheme... ',num2str(floor(P.pc)),'%'])
else
try
set(P.handles.hte(2),'string','VB updates on hyperparameters')
end
end
try
P.handles.hte(1);
catch
figure(P.handles.hfig)
P.handles.hte(1) = uicontrol('style','text',...
'units','normalized',...
'position',[0.2,0.94,0.6,0.02],...
'backgroundcolor',[1,1,1]);
end
try
set(P.handles.hte(1),'string',...
['Model evidence: p(y|m) >= ',num2str(F(end),'%10.3e\n')])
end
try
P.handles.hti;
catch
figure(P.handles.hfig)
P.handles.hti = uicontrol('style','text',...
'units','normalized',...
'position',[0.3,0.97,0.4,0.02],...
'backgroundcolor',[1,1,1],...
'string',['VB ECD inversion: trial #',num2str(P.ltr(P.ii))]);
end
drawnow
function back2defaults(e1,e2)
hf = spm_figure('FindWin','Graphics');
P = get(hf,'userdata');
try
set(hf,'colormap',P.handles.SPMdefaults.col);
set(hf,'renderer',P.handles.SPMdefaults.renderer);
end
|
