JiscMail Logo
Email discussion lists for the UK Education and Research communities

Help for SPM Archives


SPM Archives

SPM Archives


SPM@JISCMAIL.AC.UK


View:

Message:

[

First

|

Previous

|

Next

|

Last

]

By Topic:

[

First

|

Previous

|

Next

|

Last

]

By Author:

[

First

|

Previous

|

Next

|

Last

]

Font:

Proportional Font

LISTSERV Archives

LISTSERV Archives

SPM Home

SPM Home

SPM  February 2011

SPM February 2011

Options

Subscribe or Unsubscribe

Subscribe or Unsubscribe

Log In

Log In

Get Password

Get Password

Subject:

Re: Error occurred when doing dipole analysis by using SPM8

From:

Gareth Barnes <[log in to unmask]>

Reply-To:

Gareth Barnes <[log in to unmask]>

Date:

Tue, 1 Feb 2011 14:22:02 -0000

Content-Type:

multipart/mixed

Parts/Attachments:

Parts/Attachments

text/plain (93 lines) , spm_eeg_inv_vbecd_gui.m (814 lines)

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

Top of Message | Previous Page | Permalink

JiscMail Tools


RSS Feeds and Sharing


Advanced Options


Archives

July 2018
June 2018
May 2018
April 2018
March 2018
February 2018
January 2018
December 2017
November 2017
October 2017
September 2017
August 2017
July 2017
June 2017
May 2017
April 2017
March 2017
February 2017
January 2017
December 2016
November 2016
October 2016
September 2016
August 2016
July 2016
June 2016
May 2016
April 2016
March 2016
February 2016
January 2016
December 2015
November 2015
October 2015
September 2015
August 2015
July 2015
June 2015
May 2015
April 2015
March 2015
February 2015
January 2015
December 2014
November 2014
October 2014
September 2014
August 2014
July 2014
June 2014
May 2014
April 2014
March 2014
February 2014
January 2014
December 2013
November 2013
October 2013
September 2013
August 2013
July 2013
June 2013
May 2013
April 2013
March 2013
February 2013
January 2013
December 2012
November 2012
October 2012
September 2012
August 2012
July 2012
June 2012
May 2012
April 2012
March 2012
February 2012
January 2012
December 2011
November 2011
October 2011
September 2011
August 2011
July 2011
June 2011
May 2011
April 2011
March 2011
February 2011
January 2011
December 2010
November 2010
October 2010
September 2010
August 2010
July 2010
June 2010
May 2010
April 2010
March 2010
February 2010
January 2010
December 2009
November 2009
October 2009
September 2009
August 2009
July 2009
June 2009
May 2009
April 2009
March 2009
February 2009
January 2009
December 2008
November 2008
October 2008
September 2008
August 2008
July 2008
June 2008
May 2008
April 2008
March 2008
February 2008
January 2008
December 2007
November 2007
October 2007
September 2007
August 2007
July 2007
June 2007
May 2007
April 2007
March 2007
February 2007
January 2007
2006
2005
2004
2003
2002
2001
2000
1999
1998


JiscMail is a Jisc service.

View our service policies at https://www.jiscmail.ac.uk/policyandsecurity/ and Jisc's privacy policy at https://www.jisc.ac.uk/website/privacy-notice

Secured by F-Secure Anti-Virus CataList Email List Search Powered by the LISTSERV Email List Manager