123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737 |
- function varargout = spm_SpUtil(varargin)
- % Space matrix utilities
- % FORMAT varargout = spm_SpUtil(action,varargin)
- %
- %_______________________________________________________________________
- %
- % spm_SpUtil is a multi-function function containing various utilities
- % for Design matrix and contrast construction and manipulation. In
- % general, it accepts design matrices as plain matrices or as space
- % structures setup by spm_sp.
- %
- % Many of the space utilities are computed using an SVD of the design
- % matrix. The advantage of using space structures is that the svd of
- % the design matrix is stored in the space structure, thereby saving
- % unnecessary repeated computation of the SVD. This presents a
- % considerable efficiency gain for large design matrices.
- %
- % Note that when space structures are passed as arguments is is
- % assummed that their basic fields are filled in. See spm_sp for
- % details of (design) space structures and their manipulation.
- %
- % Quick Reference :
- %---------------------
- % ('isCon',x,c) :
- % ('allCon',x,c) :
- % ('ConR',x,c) :
- % ('ConO',x,c) :
- % ('size',x,dim) :
- % ('iX0check',i0,sL) :
- %---------------------
- % ('i0->c',x,i0) : Out : c
- % ('c->Tsp',x,c) : Out : [X1o [X0]]
- % ('+c->Tsp',x,c) : Out : [ukX1o [ukX0]]
- % ('i0->x1o',x,i0) : Use ('i0->c',x,i0) and ('c->Tsp',X,c)
- % ('+i0->x1o',x,i0) : Use ('i0->c',x,i0) and ('+c->Tsp',X,c)
- % ('X0->c',x,X0) :~
- % ('+X0->c',x,cukX0) :~
- %---------------------
- % ('trRV',x[,V]) :
- % ('trMV',x[,V]) :
- % ('i0->edf',x,i0,V) :
- %
- %---------------------
- %
- % Improvement compared to the spm99 beta version :
- %
- % Improvements in df computation using spm_SpUtil('trRV',x[,V]) and
- % spm_SpUtil('trMV',sX [,V]). The degrees of freedom computation requires
- % in general that the trace of RV and of RVRV be computed, where R is a
- % projector onto either a sub space of the design space or the residual
- % space, namely the space that is orthogonal to the design space. V is
- % the (estimated or assumed) variance covariance matrix and is a number
- % of scans by number of scans matrix which can be huge in some cases. We
- % have (thanks to S Rouquette and JB) speed up this computation
- % by using matlab built in functions of the frobenius norm and some theorems
- % on trace computations.
- %
- % ======================================================================
- %
- % FORMAT i = spm_SpUtil('isCon',x,c)
- % Tests whether weight vectors specify contrasts
- % x - Design matrix X, or space structure of X
- % c - contrast matrix (I.e. matrix of contrast weights, contrasts in columns)
- % Must have column dimension matching that of X
- % [defaults to eye(size(X,2)) to test uniqueness of parameter estimates]
- % i - logical row vector indiciating estimability of contrasts in c
- %
- % A linear combination of the parameter estimates is a contrast if and
- % only if the weight vector is in the space spanned by the rows of X.
- %
- % The algorithm works by regressing the contrast weight vectors using
- % design matrix X' (X transposed). Any contrast weight vectors will be
- % fitted exactly by this procedure, leaving zero residual. Parameter
- % tol is the tolerance applied when searching for zero residuals.
- %
- % Christensen R (1996)
- % "Plane Answers to Complex Questions"
- % 2nd Ed. Springer-Verlag, New York
- %
- % Andrade A, Paradis AL, Rouquette S and Poline JB, NeuroImage 9, 1999
- % ----------------
- %
- % FORMAT i = spm_SpUtil('allCon',x,c)
- % Tests whether all weight vectors specify contrasts:
- % Same as all(spm_SpUtil('isCon',x,c)).
- %
- % ----------------
- %
- % FORMAT r = spm_SpUtil('ConR',x,c)
- % Assess orthogonality of contrasts (wirit the data)
- % x - Design matrix X, or space structure of X
- % c - contrast matrix (I.e. matrix of contrast weights, contrasts in columns)
- % Must have column dimension matching that of X
- % defaults to eye(size(X,2)) to test independence of parameter estimates
- % r - Contrast correlation matrix, of dimension the number of contrasts.
- %
- % For the general linear model Y = X*B + E, a contrast weight vector c
- % defines a contrast c*B. This is estimated by c*b, where b are the
- % least squares estimates of B given by b=pinv(X)*Y. Thus, c*b = w*Y,
- % where weight vector w is given by w=c*pinv(X); Since the data are
- % assummed independent, two contrasts are indpendent if the
- % corresponding weight vectors are orthogonal.
- %
- % r is the matrix of normalised inner products between the weight
- % vectors corresponding to the contrasts. For iid E, r is the
- % correlation matrix of the contrasts.
- %
- % The logical matrix ~r will be true for orthogonal pairs of contrasts.
- %
- % ----------------
- %
- % FORMAT r = spm_SpUtil('ConO',x,c)
- % Assess orthogonality of contrasts (wirit the data)
- % x - Design matrix X, or space structure of X
- % c - contrast matrix (I.e. matrix of contrast weights, contrasts in columns)
- % Must have column dimension matching that of X
- % [defaults to eye(size(X,2)) to test uniqueness of parameter estimates]
- % r - Contrast orthogonality matrix, of dimension the number of contrasts.
- %
- % This is the same as ~spm_SpUtil('ConR',X,c), but uses a quicker
- % algorithm by looking at the orthogonality of the subspaces of the
- % design space which are implied by the contrasts:
- % r = abs(c*X'*X*c')<tol
- %
- % ----------------
- %
- % FORMAT c = spm_SpUtil('i0->c',x,i0)
- % Return F-contrast for specified design matrix partition
- % x - Design matrix X, or space structure of X
- % i0 - column indices of null hypothesis design matrix
- %
- % This functionality returns a rank n mxp matrix of contrasts suitable
- % for an extra-sum-of-squares F-test comparing the design X, with a
- % reduced design. The design matrix for the reduced design is X0 =
- % X(:,i0), a reduction of n degrees of freedom.
- %
- % The algorithm, due to J-B, and derived from Christensen, computes the
- % contrasts as an orthonormal basis set for the rows of the
- % hypothesised redundant columns of the design matrix, after
- % orthogonalisation with respect to X0. For non-unique designs, there
- % are a variety of ways to produce equivalent F-contrasts. This method
- % produces contrasts with non-zero weights only for the hypothesised
- % redundant columns.
- %
- % ----------------
- %
- % case {'x0->c'} %-
- % FORMAT c = spm_SpUtil('X0->c',sX,X0)
- % ----------------
- %
- % FORMAT [X1,X0] = spm_SpUtil('c->TSp',X,c)
- % Orthogonalised partitioning of design space implied by F-contrast
- % x - Design matrix X, or space structure of X
- % c - contrast matrix (I.e. matrix of contrast weights, contrasts in columns)
- % Must have column dimension matching that of X
- % X1o - contrast space - design matrix corresponding according to contrast
- % (orthogonalised wirit X0)
- % X0 - matrix reduced according to null hypothesis
- % (of same size as X but rank deficient)
- % FORMAT [uX1,uX0] = spm_SpUtil('c->TSp+',X,c)
- % + version to deal with the X1o and X0 partitions in the "uk basis"
- %
- % ( Note that unless X0 is reduced to a set of linearely independant )
- % ( vectors, c will only be contained in the null space of X0. If X0 )
- % ( is "reduced", then the "parent" space of c must be reduced as well )
- % ( for c to be the actual null space of X0. )
- %
- % This functionality returns a design matrix subpartition whose columns
- % span the hypothesised null design space of a given contrast. Note
- % that X1 is orthogonal(ised) to X0, reflecting the situation when an
- % F-contrast is tested using the extra sum-of-squares principle (when
- % the extra distance in the hypothesised null space is measured
- % orthogonal to the space of X0).
- %
- % Note that the null space design matrix will probably not be a simple
- % sub-partition of the full design matrix, although the space spanned
- % will be the same.
- %
- % ----------------
- %
- % FORMAT X1 = spm_SpUtil('i0->x1o',X,i0)
- % x - Design matrix X, or space structure of X
- % i0 - Columns of X that make up X0 - the reduced model (Ho:B1=0)
- % X1 - Hypothesised null design space, i.e. that part of X orthogonal to X0
- % This offers the same functionality as the 'c->TSp' option, but for
- % simple reduced models formed from the columns of X.
- %
- % FORMAT X1 = spm_SpUtil('i0->x1o+',X,i0)
- % + version to deal with the X1o and X0 partitions in the "uk basis"
- %
- % ----------------
- %
- % FORMAT [trRV,trRVRV] = spm_SpUtil('trRV',x[,V])
- % trace(RV) & trace(RVRV) - used in df calculation
- % x - Design matrix X, or space structure of X
- % V - V matrix [defult eye] (trRV == trRVRV if V==eye, since R idempotent)
- % trRV - trace(R*V), computed efficiently
- % trRVRV - trace(R*V*R*V), computed efficiently
- % This uses the Karl's cunning understanding of the trace:
- % (tr(A*B) = sum(sum(A'*B)).
- % If the space of X is set, then algorithm uses x.u to avoid extra computation.
- %
- % ----------------
- %
- % FORMAT [trMV, trMVMV]] = spm_SpUtil('trMV',x[,V])
- % trace(MV) & trace(MVMV) if two ouput arguments.
- % x - Design matrix X, or space structure of X
- % V - V matrix [defult eye] (trMV == trMVMV if V==eye, since M idempotent)
- % trMV - trace(M*V), computed efficiently
- % trMVMV - trace(M*V*M*V), computed efficiently
- % Again, this uses the Karl's cunning understanding of the trace:
- % (tr(A*B) = sum(sum(A'.*B)).
- % If the space of X is set, then algorithm uses x.u to avoid extra computation.
- %
- % ----------------
- %
- % OBSOLETE use FcUtil('H') for spm_SpUtil('c->H',x,c)
- % Extra sum of squares matrix O for beta's from contrast
- % x - Design matrix X, or space structure of X
- % c - contrast matrix (I.e. matrix of contrast weights, contrasts in columns)
- % Must have column dimension matching that of X
- % O - Matrix such that b'*O*b = extra sum of squares for F-test of contrast c
- %
- % ----------------
- %
- % OBSOLETE use spm_sp('=='...) for spm_SpUtil('c==X1o',x,c) {or 'cxpequi'}
- % x - Design matrix X, or space structure of X
- % c - contrast matrix (I.e. matrix of contrast weights, contrasts in columns)
- % Must have column dimension matching that of X
- % b - True is c is a spanning set for space of X
- % (I.e. if contrast and space test the same thing)
- %
- % ----------------
- %
- % FORMAT [df1,df2] = spm_SpUtil('i0->edf',x,i0,V) {or 'edf'}
- % (effective) df1 and df2 the residual df for the projector onto the
- % null space of x' (residual forming projector) and the numerator of
- % the F-test where i0 are the columns for the null hypothesis model.
- % x - Design matrix X, or space structure of X
- % i0 - Columns of X corresponding to X0 partition X = [X1,X0] & with
- % parameters B = [B1;B0]. Ho:B1=0
- % V - V matrix
- %
- % ----------------
- %
- % FORMAT sz = spm_SpUtil('size',x,dim)
- % FORMAT [sz1,sz2,...] = spm_SpUtil('size',x)
- % Returns size of design matrix
- % (Like MatLab's `size`, but copes with design matrices inside structures.)
- % x - Design matrix X, or structure containing design matrix in field X
- % (Structure needn't be a space structure.)
- % dim - dimension which to size
- % sz - size
- %
- %_______________________________________________________________________
- % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
- % Andrew Holmes Jean-Baptiste Poline
- % $Id: spm_SpUtil.m 4137 2010-12-15 17:18:32Z guillaume $
- % (frobenius norm trick by S. Rouquette)
- %-Format arguments
- %-----------------------------------------------------------------------
- if nargin==0, error('do what? no arguments given...')
- else action = varargin{1}; end
- switch lower(action),
- case {'iscon','allcon','conr','cono'}
- %=======================================================================
- % i = spm_SpUtil('isCon',x,c)
- if nargin==0, varargout={[]}; error('isCon : no argument specified'), end;
- if nargin==1,
- varargout={[]}; warning('isCon : no contrast specified'); return;
- end;
- if ~spm_sp('isspc',varargin{2})
- sX = spm_sp('Set',varargin{2});
- else sX = varargin{2}; end
- if nargin==2, c=eye(spm_sp('size',sX,2)); else c=varargin{3}; end;
- if isempty(c), varargout={[]}; return, end
- switch lower(action)
- case 'iscon'
- varargout = { spm_sp('eachinspp',sX,c) };
- case 'allcon'
- varargout = {spm_sp('isinspp',sX,c)};
- case 'conr'
- if size(c,1) ~= spm_sp('size',sX,2)
- error('Contrast not of the right size'), end
- %-Compute inner products of data weight vectors
- % (c'b = c'*pinv(X)*Y = w'*Y
- % (=> w*w' = c'*pinv(X)*pinv(X)'*c == c'*pinv(X'*X)*c
- r = c'*spm_sp('XpX-',sX)*c;
- %-normalize by "cov(r)" to get correlations
- r = r./(sqrt(diag(r))*sqrt(diag(r))');
- r(abs(r) < sX.tol)=0; %-set near-zeros to zero
- varargout = {r}; %-return r
- case 'cono'
- %-This is the same as ~spm_SpUtil('ConR',x,c), and so returns
- % the contrast orthogonality (though not their corelations).
- varargout = { abs(c'* spm_sp('XpX',sX) *c) < sX.tol};
- end
- case {'+c->tsp','c->tsp'} %- Ortho. partitioning implied by F-contrast
- %=======================================================================
- % spm_SpUtil('c->Tsp',sX,c)
- % + version of 'c->tsp'.
- % The + version returns the same in the base u(:,1:r).
- %--------- begin argument check ------------------------------
- if nargin ~= 3, error(['Wrong number of arguments in ' action])
- else sX = varargin{2}; c = varargin{3}; end;
- if nargout > 2, error(['Too many output arguments in ' action]), end;
- if ~spm_sp('isspc',sX), sX = spm_sp('set',varargin{2}); end;
- if sX.rk == 0, error('c->Tsp null rank sX == 0'); end;
- if ~isempty(c) && spm_sp('size',sX,2) ~= size(c,1),
- error(' c->TSp matrix & contrast dimensions don''t match');
- end
- %--------- end argument check ---------------------------------
- %- project c onto the space of X' if needed
- %-------------------------------------------
- if ~isempty(c) && ~spm_sp('isinspp',sX,c),
- warning([sprintf('\n') 'c is not a proper contrast in ' action ...
- ' in ' mfilename sprintf('\n') '!!! projecting...' ]);
- disp('from'), c, disp('to'), c = spm_sp('oPp:',sX,c)
- end
- cukFlag = strcmp(lower(action),'+c->tsp');
- switch nargout
- % case 0
- % warning(['no output demanded in ' mfilename ' ' action])
- case {0,1}
- if ~isempty(c) && any(any(c)) %- c not empty & not null
- if cukFlag, varargout = { spm_sp('cukxp-:',sX,c) };
- else varargout = { spm_sp('xp-:',sX,c) };
- end
- else if isempty(c), varargout = { [] }; %- c empty
- else %- c null
- if cukFlag, varargout = { spm_sp('cukx',sX,c) };
- else varargout = { spm_sp('x',sX)*c };
- end
- end
- end
- case 2
- if ~isempty(c) && any(any(c)) %- not empty and not null
- if cukFlag,
- varargout = {
- spm_sp('cukxp-:',sX,c), ... %- X1o
- spm_sp('cukx',sX,spm_sp('r',spm_sp('set',c))) }; %- X0
- else
- varargout = {
- spm_sp(':',sX, spm_sp('xp-:',sX,c)), ... %- X1o
- spm_sp(':',sX, ...
- spm_sp('x',sX)*spm_sp('r',spm_sp('set',c))) }; %- X0
- end
- else
- if isempty(c), %- empty
- if cukFlag, varargout = { [], spm_sp('cukx',sX) };
- else varargout = { [], spm_sp('x',sX) };
- end
- else %- null
- if cukFlag,
- varargout = { spm_sp(':',sX,spm_sp('cukx',sX,c)), ...
- spm_sp(':',sX,spm_sp('cukx',sX)) };
- else
- varargout = { spm_sp('x',sX)*c, spm_sp('x',sX)};
- end
- end;
- end
- otherwise
- error(['wrong number of output argument in ' action]);
- end
- case {'i0->x1o','+i0->x1o'} %- Space tested whilst keeping size of X(i0)
- %=======================================================================
- % X1o = spm_SpUtil('i0->X1o',sX,i0)
- % arguments are checked in calls to spm_Util
- %--------------------------------------------
- if nargin<3, error('Insufficient arguments'),
- else sX = varargin{2}; i0 = varargin{3}; end;
- cukFlag = strcmp(lower(action),'+i0->x1o');
- c = spm_SpUtil('i0->c',sX,i0);
- if cukFlag,
- varargout = { spm_SpUtil('+c->TSp',sX,c) };
- else
- varargout = { spm_SpUtil('c->TSp',sX,c) };
- end
- case {'i0->c'} %-
- %=======================================================================
- % c = spm_SpUtil('i0->c',sX,i0)
- %
- % if i0 == [] : returns a proper contrast
- % if i0 == [1:size(sX.X,2)] : returns [];
- %
- %- Algorithm : avoids the pinv(X0) and insure consistency
- %- Get the estimable parts of c0 and c1
- %- remove from c1_estimable the estimable part of c0.
- %- Use the rotation making X1o orthog. to X0.
- %- i0 is checked when Fc is created
- %- If i0 defines a space that is the space of X but with
- %- fewer vectors, c is null.
- %--------- begin argument check --------------------------------
- if nargin<3, error('Insufficient arguments'),
- else sX = varargin{2}; i0 = varargin{3}; end;
- if ~spm_sp('isspc',sX), sX = spm_sp('set',varargin{2}); end;
- if spm_sp('rk',sX) == 0, error('i0->c null rank sX == 0'); end;
- sL = spm_sp('size',sX,2);
- i0 = sf_check_i0(i0,sL);
- %--------- end argument check ----------------------------------
- c0 = eye(sL); c0 = c0(:,i0);
- c1 = eye(sL); c1 = c1(:,setdiff(1:sL,i0));
- %- try to avoid the matlab error when doing svd of matrices with
- %- high multiplicities. (svd convergence pb)
- if ~ spm_sp('isinspp',sX,c0), c0 = spm_sp('oPp:',sX,c0); end;
- if ~ spm_sp('isinspp',sX,c1), c1 = spm_sp('oPp:',sX,c1); end;
- if ~isempty(c1)
- if ~isempty(c0)
- %- varargout = { spm_sp('res',spm_sp('set',opp*c0),opp*c1) };
- %- varargout = { c1 - c0*pinv(c0)*c1 }; NB: matlab pinv uses
- %- svd: will fail if spm_sp('set') fails.
- varargout = { spm_sp('r:',spm_sp('set',c0),c1) };
- else varargout = { spm_sp('xpx',sX) }; end;
- else
- varargout = { [] }; %- not zeros(sL,1) : this is return when
- %- appropriate
- end
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- case {'+x0->c','x0->c'} %-
- %=======================================================================
- % c = spm_SpUtil('X0->c',sX,X0)
- % c = spm_SpUtil('+X0->c',sX,cukX0)
- % + version of 'x0->c'.
- % The + version returns the same in the base u(:,1:r).
- warning('Not tested for release - provided for completeness');
- cukFlag = strcmp(lower(action),'+x0->c');
- %--------- begin argument check ---------
- if nargin<3, error('Insufficient arguments'),
- else
- sX = varargin{2};
- if cukFlag, cukX0 = varargin{3}; else X0 = varargin{3}; end
- end
- if ~spm_sp('isspc',sX), sX = spm_sp('set',varargin{2}); end
- if spm_sp('rk',sX) == 0, error(['null rank sX == 0 in ' action]); end
- if cukFlag
- if ~isempty(cukX0) && spm_sp('rk',sX) ~= size(cukX0,1),
- cukX0, spm_sp('rk',sX),
- error(['cukX0 of wrong size ' mfilename ' ' action]), end
- else
- if ~isempty(X0) && spm_sp('size',sX,1) ~= size(X0,1),
- X0, spm_sp('size',sX,1),
- error(['X0 of wrong size ' mfilename ' ' action]),X0, end
- end
- %--------- end argument check ---------
- if cukFlag
- if isempty(cukX0), X0 = []; else X0 = spm_sp('ox',sX)*cukX0; end
- end
- varargout = { sf_X0_2_c(X0,sX) };
- case {'c->h','betarc'} %-Extra sum of squares matrix for beta's from
- %- contrast : use F-contrast if possible
- %=======================================================================
- % H = spm_SpUtil('c->H',sX,c)
- error(' Obsolete : Use F-contrast utilities ''H'' or ''Hsqr''... ');
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- %=======================================================================
- %=======================================================================
- % trace part
- %=======================================================================
- %=======================================================================
- %
- case 'trrv' %-Traces for (effective) df calculation
- %=======================================================================
- % [trRV,trRVRV]= spm_SpUtil('trRV',x[,V])
- if nargin == 1, error('insufficient arguments');
- else sX = varargin{2}; end;
- if ~spm_sp('isspc',sX), sX = spm_sp('Set',sX); end;
- rk = spm_sp('rk',sX);
- sL = spm_sp('size',sX,1);
- if sL == 0,
- warning('space with no dimension ');
- if nargout==1, varargout = {[]};
- else varargout = {[], []}; end
- else
- if nargin > 2 && ~isempty(varargin{3})
- V = varargin{3};
- u = sX.u(:,1:rk);
- clear sX;
- if nargout==1
- %-only trRV needed
- if rk==0 || isempty(rk), trMV = 0;
- else trMV = sum(sum( u .* (V*u) ));
- end
- varargout = { trace(V) - trMV};
- else
- %-trRVRV is needed as well
- if rk==0 || isempty(rk),
- trMV = 0;
- trRVRV = (norm(V,'fro'))^2;
- trV = trace(V);
- clear V u
- else
- Vu = V*u;
- trV = trace(V);
- trRVRV = (norm(V,'fro'))^2;
- clear V;
- trRVRV = trRVRV - 2*(norm(Vu,'fro'))^2;
- trRVRV = trRVRV + (norm(u'*Vu,'fro'))^2;
- trMV = sum(sum( u .* Vu ));
- clear u Vu
- end
- varargout = {(trV - trMV), trRVRV};
- end
-
- else %- nargin == 2 | isempty(varargin{3})
- if nargout==1
- if rk==0 || isempty(rk), varargout = {sL};
- else varargout = {sL - rk};
- end
- else
- if rk==0 || isempty(rk), varargout = {sL,sL};
- else varargout = {sL - rk, sL - rk};
- end
- end
- end
- end
- case 'trmv' %-Traces for (effective) Fdf calculation
- %=======================================================================
- % [trMV, trMVMV]] = spm_SpUtil('trMV',sX [,V])
- %
- % NB : When V is given empty, the routine asssumes it's identity
- % This is used in spm_FcUtil.
- if nargin == 1, error('insufficient arguments');
- else sX = varargin{2}; end;
- if ~spm_sp('isspc',sX), sX = spm_sp('Set',sX); end;
- rk = spm_sp('rk',sX);
- if isempty(rk)
- warning('Rank is empty');
- if nargout==1, varargout = {[]};
- else varargout = {[], []}; end
- return;
- elseif rk==0, warning('Rank is null in spm_SpUtil trMV ');
- if nargout==1, varargout = {0};
- else varargout = {0, 0}; end
- return;
- end;
- if nargin > 2 && ~isempty(varargin{3}) %- V provided, and assumed correct !
- V = varargin{3};
- u = sX.u(:,1:rk);
- clear sX;
- if nargout==1
- %-only trMV needed
- trMV = sum(sum(u' .* (u'*V) ));
- varargout = {trMV};
- else
- %-trMVMV is needed as well
- Vu = V*u;
- clear V
- trMV = sum(sum( u .* Vu ));
- trMVMV = (norm(u'*Vu,'fro'))^2;
- clear u Vu
- varargout = {trMV, trMVMV};
- end
- else % nargin == 2 | isempty(varargin{3}) %-no V specified: trMV == trMVMV
- if nargout==1
- varargout = {rk};
- else
- varargout = {rk, rk};
- end
- end
-
- case {'i0->edf','edf'} %-Effective F degrees of freedom
- %=======================================================================
- % [df1,df2] = spm_SpUtil('i0->edf',sX,i0,V)
- %-----------------------------------------------------------------------
- %--------- begin argument check ----------------------------------------
- if nargin<3, error('insufficient arguments'),
- else i0 = varargin{3}; sX = varargin{2}; end
- if ~spm_sp('isspc',sX), sX = spm_sp('Set',sX); end;
- i0 = sf_check_i0(i0,spm_sp('size',sX,2));
- if nargin == 4, V=varargin{4}; else V = eye(spm_sp('size',sX,1)); end;
- if nargin>4, error('Too many input arguments'), end;
- %--------- end argument check ------------------------------------------
- warning(' Use F-contrast utilities if possible ... ');
- [trRV,trRVRV] = spm_SpUtil('trRV', sX, V);
- [trMpV,trMpVMpV] = spm_SpUtil('trMV',spm_SpUtil('i0->x1o',sX, i0),V);
- varargout = {trMpV^2/trMpVMpV, trRV^2/trRVRV};
- %=======================================================================
- %=======================================================================
- % Utilities
- %=======================================================================
- %=======================================================================
- case 'size' %-Size of design matrix
- %=======================================================================
- % sz = spm_SpUtil('size',x,dim)
- if nargin<3, dim=[]; else dim = varargin{3}; end
- if nargin<2, error('insufficient arguments'), end
- if isstruct(varargin{2})
- if isfield(varargin{2},'X')
- sz = size(varargin{2}.X);
- else error('no X field'); end;
- else
- sz = size(varargin{2});
- end
- if ~isempty(dim)
- if dim>length(sz), sz = 1; else sz = sz(dim); end
- varargout = {sz};
- elseif nargout>1
- varargout = cell(1,min(nargout,length(sz)));
- for i=1:min(nargout,length(sz)), varargout{i} = sz(i); end
- else
- varargout = {sz};
- end
- case 'ix0check' %-
- %=======================================================================
- % i0c = spm_SpUtil('iX0check',i0,sL)
- if nargin<3, error('insufficient arguments'),
- else i0 = varargin{2}; sL = varargin{3}; end;
- varargout = {sf_check_i0(i0,sL)};
- otherwise
- %=======================================================================
- error('Unknown action string in spm_SpUtil')
- %=======================================================================
- end
- %=======================================================================
- function i0c = sf_check_i0(i0,sL)
- % NB : [] = sf_check_i0([],SL);
- %
- if all(ismember(i0,[0,1])) && length(i0(:))==sL, i0c=find(i0);
- elseif ~isempty(i0) && any(floor(i0)~=i0) || any(i0<1) || any(i0>sL)
- error('logical mask or vector of column indices required')
- else i0c = i0; end
- %=======================================================================
- function c = sf_X0_2_c(X0,sX)
- %
- %- Algorithm to avoids the pinv(X0) and insure consistency
- %- Get a contrast that span the space of X0 and is estimable
- %- Get the orthogonal complement and project onto the estimable space
- %- Strip zeros columns and use the rotation making X1o orthog. to X0
- % !!! tolerance dealing ?
- if ~isempty(X0)
- sc0 = spm_sp('set',spm_sp('x-',sX,X0));
- if sc0.rk
- c = spm_sp('oPp:',sX,spm_sp('r',sc0));
- else
- c = spm_sp('oPp',sX);
- end;
- c = c(:,any(c));
- sL = spm_sp('size',sX,2);
- %- why the "& size(X0,2) ~= sL" !!!?
- if isempty(c) && size(X0,2) ~= sL
- c = zeros(sL,1);
- end
- else
- c = spm_sp('xpx',sX);
- end
- %- c = spm_sp('r',sc0,spm_sp('oxp',sX)); would also works.
|