function varargout = spm_FcUtil(varargin)

% Contrast utilities

% FORMAT varargout = spm_FcUtil(action,varargin)

%_______________________________________________________________________

%

% spm_FcUtil is a multi-function function containing various utilities

% for contrast construction and manipulation. In general, it accepts

% design matrices as plain matrices or as space structures setup by

% spm_sp (that is preferable in general).

%

% The use of spm_FcUtil should help with robustness issues and

% maintainability of SPM. % Note that when space structures are passed

% as arguments is is assumed that their basic fields are filled in.

% See spm_sp for details of (design) space structures and their

% manipulation.

%

%

% ======================================================================

% case 'fconfields' %- fields of F contrast

% Fc = spm_FcUtil('FconFields')

%

%- simply returns the fields of a contrast structure.

%

%=======================================================================

% case 'set' %- Create an F contrast

% Fc = spm_FcUtil('Set',name, STAT, set_action, value, sX)

%

%- Set will fill in the contrast structure, in particular

%- c (in the contrast space), X1o (the space actually tested) and

%- X0 (the space left untested), such that space([X1o X0]) == sX.

%- STAT is either 'F' or 'T';

%- name is a string describing the contrast.

%

%- There are three ways to set a contrast :

%- set_action is 'c','c ' : value can then be zeros.

%- dimensions are in X',

%- f c is used, value is projected onto sX';

%- iX0 is set to 'c' or 'c ';

%- set_action is 'iX0' : defines the indices of the columns

%- that will not be tested. Can be empty.

%- set_action is 'X0' : defines the space that will remain

%- unchanged. The orthogonal complement is

%- tested; iX0 is set to 'X0';

%-

%=======================================================================

% case 'isfcon' %- Is it an F contrast ?

% b = spm_FcUtil('IsFcon',Fc)

%

%=======================================================================

% case 'fconedf' %- F contrast edf

% [edf_tsp edf_Xsp] = spm_FcUtil('FconEdf', Fc, sX [, V])

%

%- compute the effective degrees of freedom of the numerator edf_tsp

%- and (optionally) the denominator edf_Xsp of the contrast.

%

%=======================================================================

% case 'hsqr' %-Extra sum of squares sqr matrix for beta's from contrast

% hsqr = spm_FcUtil('Hsqr',Fc, sX)

%

%- This computes the matrix hsqr such that a the numerator of an F test

%- will be beta'*hsqr'*hsqr*beta

%

%=======================================================================

% case 'h' %-Extra sum of squares matrix for beta's from contrast

% H = spm_FcUtil('H',Fc, sX)

%

%- This computes the matrix H such that a the numerator of an F test

%- will be beta'*H*beta

%-

%=======================================================================

% case 'yc' %- Fitted data corrected for confounds defined by Fc

% Yc = spm_FcUtil('Yc',Fc, sX, b)

%

%- Input : b : the betas

%- Returns the corrected data Yc for given contrast. Y = Yc Y0 error

%

%=======================================================================

% case 'y0' %- Confounds data defined by Fc

% Y0 = spm_FcUtil('Y0',Fc, sX, b)

%

%- Input : b : the betas

%- Returns the confound data Y0 for a given contrast. Y = Yc Y0 error

%

%=======================================================================

% case {'|_'} %- Fc orthogonalisation

% Fc = spm_FcUtil('|_',Fc1, sX, Fc2)

%

%- Orthogonolise a (list of) contrasts Fc1 wrt a (list of) contrast Fc2

%- such that the space these contrasts test are orthogonal.

%- If contrasts are not estimable contrasts, works with the estimable

%- part. In any case, returns estimable contrasts.

%

%=======================================================================

% case {'|_?'} %- Are contrasts orthogonals

% b = spm_FcUtil('|_?',Fc1, sX [, Fc2])

%

%- Tests whether a (list of) contrast is orthogonal. Works with the

%- estimable part if they are not estimable. With only one argument,

%- tests whether the list is made of orthogonal contrasts. With Fc2

%- provided, tests whether the two (list of) contrast are orthogonal.

%

%=======================================================================

% case 'in' %- Fc1 is in list of contrasts Fc2

% [iFc2 iFc1] = spm_FcUtil('In', Fc1, sX, Fc2)

%

%- Tests whether a (list of) contrast Fc1 is in a list of contrast Fc2.

%- returns the indices iFc2 where element of Fc1 have been found

%- in Fc2 and the indices iFc1 of the element of Fc1 found in Fc2.

%- These indices are not necessarily unique.

%

%=======================================================================

% case '~unique' %- Fc list unique

% idx = spm_FcUtil('~unique', Fc, sX)

%

%- returns indices ofredundant contrasts in Fc

%- such that Fc(idx) = [] makes Fc unique.

%

%=======================================================================

% case {'0|[]','[]|0'} %- Fc is null or empty

% b = spm_FcUtil('0|[]', Fc, sX)

%

%- NB : for the "null" part, checks if the contrast is in the null space

%- of sX (completely non estimable !)

%=======================================================================

%_______________________________________________________________________

% Jean-Baptiste Poline

% Copyright (C) 1999-2022 Wellcome Centre for Human Neuroimaging

%-Format arguments

%-----------------------------------------------------------------------

if nargin==0, error('do what? no arguments given...')

else action = varargin{1}; end

switch lower(action),

case 'fconfields' %- fields of F contrast

%=======================================================================

% Fc = spm_FcUtil('FconFields')

if nargout > 1, error('Too many output arguments: FconFields'), end

if nargin > 1, error('Too many input arguments: FconFields'), end

varargout = {sf_FconFields};

case {'set','v1set'} %- Create an F contrast

%=======================================================================

% Fc = spm_FcUtil('Set',name, STAT, set_action, value, sX)

%

% Sets the contrast structure with set_action either 'c', 'X0' or 'iX0'

% resp. for a contrast, the null hyp. space or the indices of which.

% STAT can be 'T' or 'F'.

%

% If not set by iX0 (in which case field .iX0 contains the indices),

% field .iX0 is set as a string containing the set_action: {'X0','c','c ','ukX0'}

%

% if STAT is T, then set_action should be 'c' or 'c '

% (at the moment, just a warning...)

% if STAT is T and set_action is 'c' or 'c ', then

% checks whether it is a real T.

%

% 'v1set' is NOT provided for backward compatibility so far ...

%-check # arguments...

%--------------------------------------------------------------------------

if nargin<6, error('insufficient arguments'), end

if nargout > 1, error('Too many output arguments Set'), end

%-check arguments...

%--------------------------------------------------------------------------

if ~ischar(varargin{2}), error('~ischar(name)'), end

if ~(varargin{3}=='F'||varargin{3}=='T'||varargin{3}=='P'),

error('~(STAT==F|STAT==T|STAT==P)'), end

if ~ischar(varargin{4}), error('~ischar(varargin{4})');

else set_action = varargin{4}; end

sX = varargin{6};

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end

if isempty(sX.X), error('Empty space X in Set'); end

Fc = sf_FconFields;

%- use the name as a flag to insure that F-contrast has been

%- properly created;

Fc.name = varargin{2};

Fc.STAT = varargin{3};

if Fc.STAT=='T' && ~(any(strcmp(set_action,{'c ','c'})))

warning('enter T stat with contrast - here no check rank == 1');

end

[sC,sL] = spm_sp('size',sX);

%- allow to define the contrast the old (version 1) way ?

%- NO. v1 = strcmp(action,'v1set');

switch set_action,

case {'c','c '}

Fc.iX0 = set_action;

c = spm_sp(':', sX, varargin{5});

if isempty(c)

[Fc.X1o.ukX1o,Fc.X0.ukX0] = spm_SpUtil(' c->Tsp',sX,[]);

%- v1 [Fc.X1o Fc.X0] = spm_SpUtil('c->Tsp',sX,[]);

Fc.c = c;

elseif size(c,1) ~= sL,

error(['not contrast dim. in ' mfilename ' ' set_action]);

else

if strcmp(set_action,'c ')

if ~spm_sp('isinspp',sX,c), c = spm_sp('oPp:',sX,c); end

end;

if Fc.STAT=='T' && ~sf_is_T(sX,c)

%- Could be make more self-correcting by giving back an F

error('trying to define a t that looks like an F');

end

Fc.c = c;

[Fc.X1o.ukX1o,Fc.X0.ukX0] = spm_SpUtil(' c->Tsp',sX,c);

%- v1 [Fc.X1o Fc.X0] = spm_SpUtil('c->Tsp',sX,c);

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% 'option given for completeness - not for SPM use'

case {'X0'}

warning('option given for completeness - not for SPM use');

Fc.iX0 = set_action;

X0 = spm_sp(':', sX, varargin{5});

if isempty(X0),

Fc.c = spm_sp('xpx',sX);

Fc.X1o.ukX1o = spm_sp('cukx',sX);

Fc.X0.ukX0 = [];

elseif size(X0,1) ~= sC,

error('dimension of X0 wrong in Set');

else

Fc.c = spm_SpUtil('X0->c',sX,X0);

Fc.X0.ukX0 = spm_sp('ox',sX)'*X0;

Fc.X1o.ukX1o = spm_SpUtil(' c->Tsp',sX,Fc.c);

end

case 'ukX0'

warning('option given for completeness - not for SPM use');

Fc.iX0 = set_action;

if isempty(ukX0),

Fc.c = spm_sp('xpx',sX);

Fc.X1o.ukX1o = spm_sp('cukx',sX);

Fc.X0.ukX0 = [];

elseif size(ukX0,1) ~= spm_sp('rk',sX),

error('dimension of cukX0 wrong in Set');

else

Fc.c = spm_SpUtil(' X0->c',sX,ukX0);

Fc.X0.ukX0 = ukX0;

Fc.X1o.ukX1o = spm_SpUtil(' c->Tsp',sX,Fc.c);

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

case 'iX0'

iX0 = varargin{5};

iX0 = spm_SpUtil('iX0check',iX0,sL);

Fc.iX0 = iX0;

Fc.X0.ukX0 = spm_sp('ox',sX)' * spm_sp('Xi',sX,iX0);

if isempty(iX0),

Fc.c = spm_sp('xpx',sX);

Fc.X1o.ukX1o = spm_sp('cukx',sX);

else

Fc.c = spm_SpUtil('i0->c',sX,iX0);

Fc.X1o.ukX1o = spm_SpUtil(' c->Tsp',sX,Fc.c);

end

otherwise

error('wrong action in Set ');

end

varargout = {Fc};

case 'x0' % spm_FcUtil('X0',Fc,sX)

%=======================================================================

if nargin ~= 3,

error('too few/many input arguments - need 2');

else

Fc = varargin{2}; sX = varargin{3};

end

if nargout ~= 1, error('too few/many output arguments - need 1'), end

if ~sf_IsFcon(Fc), error('argument is not a contrast struct'), end

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end

varargout = {sf_X0(Fc,sX)};

case 'x1o' % spm_FcUtil('X1o',Fc,sX)

%=======================================================================

if nargin ~= 3,

error('too few/many input arguments - need 2');

else

Fc = varargin{2}; sX = varargin{3};

end

if nargout ~= 1, error('too few/many output arguments - need 1'), end

if ~sf_IsFcon(Fc), error('argument is not a contrast struct'), end

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end

varargout = {sf_X1o(Fc,sX)};

case 'isfcon' %- Is it an F contrast ?

%=======================================================================

% yes_no = spm_FcUtil('IsFcon',Fc)

if nargin~=2, error('too few/many input arguments - need 2'), end

if ~isstruct(varargin{2}), varargout={0};

else varargout = {sf_IsFcon(varargin{2})};

end

case 'fconedf' %- F contrast edf

%=======================================================================

% [edf_tsp edf_Xsp] = spm_FcUtil('FconEdf', Fc, sX [, V])

if nargin<3, error('Insufficient arguments'), end

if nargout >= 3, error('Too many output argument.'), end

Fc = varargin{2};

sX = varargin{3};

if nargin == 4, V = varargin{4}; else V = []; end

if ~sf_IsFcon(Fc), error('Fc must be Fcon'), end

if ~spm_sp('isspc',sX)

sX = spm_sp('set',sX); end

if ~sf_isempty_X1o(Fc)

[trMV, trMVMV] = spm_SpUtil('trMV',sf_X1o(Fc,sX),V);

else

trMV = 0;

trMVMV = 0;

end

if ~trMVMV, edf_tsp = 0; warning('edf_tsp = 0'),

else edf_tsp = trMV^2/trMVMV; end;

if nargout == 2

[trRV, trRVRV] = spm_SpUtil('trRV',sX,V);

if ~trRVRV, edf_Xsp = 0; warning('edf_Xsp = 0'),

else edf_Xsp = trRV^2/trRVRV; end;

varargout = {edf_tsp, edf_Xsp};

else

varargout = {edf_tsp};

end

%=======================================================================

%=======================================================================

% parts that use F contrast

%=======================================================================

%=======================================================================

%

% Quick reference : L : can be lists of ...

%-------------------------

% ('Hsqr',Fc, sX) : Out: Hsqr / ESS = b' * Hsqr' * Hsqr * b

% ('H',Fc, sX) : Out: H / ESS = b' * H * b

% ('Yc',Fc, sX, b) : Out: Y corrected = X*b - X0*X0- *Y

% ('Y0',Fc, sX, b) : Out: Y0 = X0*X0- *Y

% ('|_',LFc1, sX, LFc2) : Out: Fc1 orthog. wrt Fc2

% ('|_?',LFc1,sX [,LFc2]): Out: is Fc2 ortho to Fc1 or is Fc1 ortho ?

% ('In', LFc1, sX, LFc2) : Out: indices of Fc2 if "in", 0 otherwise

% ('~unique', LFc, sX) : Out: indices of redundant contrasts

% ('0|[]', Fc, sX) : Out: 1 if Fc is zero or empty, 0 otherwise

case 'hsqr' %-Extra sum of squares matrix for beta's from contrast

%=======================================================================

% hsqr = spm_FcUtil('Hsqr',Fc, sX)

if nargin<3, error('Insufficient arguments'), end

if nargout>1, error('Too many output argument.'), end

Fc = varargin{2};

sX = varargin{3};

if ~sf_IsFcon(Fc), error('Fc must be F-contrast'), end

if ~sf_IsSet(Fc), error('Fcon must be set'); end; %-

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

if sf_isempty_X1o(Fc)

if ~sf_isempty_X0(Fc)

%- assumes that X0 is sX.X

%- warning(' Empty X1o in spm_FcUtil(''Hsqr'',Fc,sX) ');

varargout = { zeros(1,spm_sp('size',sX,2)) };

else

error(' Fc must be set ');

end

else

varargout = { sf_Hsqr(Fc,sX) };

end

case 'h' %-Extra sum of squares matrix for beta's from contrast

%=======================================================================

% H = spm_FcUtil('H',Fc, sX)

% Empty and zeros dealing :

% This routine never returns an empty matrix.

% If sf_isempty_X1o(Fc) | isempty(Fc.c) it explicitly

% returns a zeros projection matrix.

if nargin<2, error('Insufficient arguments'), end

if nargout>1, error('Too many output argument.'), end

Fc = varargin{2};

sX = varargin{3};

if ~sf_IsFcon(Fc), error('Fc must be F-contrast'), end

if ~sf_IsSet(Fc), error('Fcon must be set'); end; %-

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

if sf_isempty_X1o(Fc)

if ~sf_isempty_X0(Fc)

%- assumes that X0 is sX.X

%- warning(' Empty X1o in spm_FcUtil(''H'',Fc,sX) ');

varargout = { zeros(spm_sp('size',sX,2)) };

else

error(' Fc must be set ');

end

else

varargout = { sf_H(Fc,sX) };

end

case 'yc' %- Fitted data corrected for confounds defined by Fc

%=======================================================================

% Yc = spm_FcUtil('Yc',Fc, sX, b)

if nargin < 4, error('Insufficient arguments'), end

if nargout > 1, error('Too many output argument.'), end

Fc = varargin{2}; sX = varargin{3}; b = varargin{4};

if ~sf_IsFcon(Fc), error('Fc must be F-contrast'), end

if ~sf_IsSet(Fc), error('Fcon must be set'); end;

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

% if ~spm_FcUtil('Rcompatible',Fc,sX), ...

% error('sX and Fc must be compatible'), end;

if spm_sp('size',sX,2) ~= size(b,1),

error('sX and b must be compatible'), end;

if sf_isempty_X1o(Fc)

if ~sf_isempty_X0(Fc)

%- if space of interest empty or null, returns zeros !

varargout = { zeros(spm_sp('size',sX,1),size(b,2)) };

else

error(' Fc must be set ');

end

else

varargout = { sf_Yc(Fc,sX,b) };

end

case 'y0' %- Fitted data corrected for confounds defined by Fc

%=======================================================================

% Y0 = spm_FcUtil('Y0',Fc, sX, b)

if nargin < 4, error('Insufficient arguments'), end

if nargout > 1, error('Too many output argument.'), end

Fc = varargin{2}; sX = varargin{3}; b = varargin{4};

if ~sf_IsFcon(Fc), error('Fc must be F-contrast'), end

if ~sf_IsSet(Fc), error('Fcon must be set'); end;

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

if spm_sp('size',sX,2) ~= size(b,1),

error('sX and b must be compatible'), end;

if sf_isempty_X1o(Fc)

if ~sf_isempty_X0(Fc)

%- if space of interest empty or null, returns zeros !

varargout = { sX.X*b };

else

error(' Fc must be set ');

end

else

varargout = { sf_Y0(Fc,sX,b) };

end

case {'|_'} %- Fc orthogonalisation

%=======================================================================

% Fc = spm_FcUtil('|_',Fc1, sX, Fc2)

% returns Fc1 orthogonolised wrt Fc2

if nargin < 4, error('Insufficient arguments'), end

if nargout > 1, error('Too many output argument.'), end

Fc1 = varargin{2}; sX = varargin{3}; Fc2 = varargin{4};

%-check arguments

%-----------------------------------------------------------------------

L1 = length(Fc1);

if ~L1, warning('no contrast given to |_'); varargout = {[]}; return; end

for i=1:L1

if ~sf_IsFcon(Fc1(i)), error('Fc1(i) must be a contrast'), end

end

L2 = length(Fc2);

if ~L2, error('must have at least a contrast in Fc2'); end

for i=1:L2

if ~sf_IsFcon(Fc2(i)), error('Fc2(i) must be a contrast'), end

end

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

%-create an F-contrast for all the Fc2

%--------------------------------------------------------------------------

str = Fc2(1).name; for i=2:L2, str = [str ' ' Fc2(i).name]; end;

Fc2 = spm_FcUtil('Set',str,'F','c ',cat(2,Fc2(:).c),sX);

if sf_isempty_X1o(Fc2) || sf_isnull(Fc2,sX)

varargout = {Fc1};

else

for i=1:L1

if sf_isempty_X1o(Fc1(i)) || sf_isnull(Fc1(i),sX)

%- Fc1(i) is an [] or 0 contrast : ortho to anything;

out(i) = Fc1(i);

else

out(i) = sf_fcortho(Fc1(i), sX, Fc2);

end

end

varargout = {out};

end

case {'|_?'} %- Are contrasts orthogonals

%=======================================================================

% b = spm_FcUtil('|_?',Fc1, sX [, Fc2])

if nargin < 3, error('Insufficient arguments'), end

Fc1 = varargin{2}; sX = varargin{3};

if nargin > 3, Fc2 = varargin{4}; else Fc2 = []; end;

if isempty(Fc1), error('give at least one non empty contrast'), end;

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

for i=1:length(Fc1)

if ~sf_IsFcon(Fc1(i)), error('Fc1(i) must be a contrast'), end

end

for i=1:length(Fc2)

if ~sf_IsFcon(Fc2(i)), error('Fc2(i) must be a contrast'), end

end

varargout = { sf_Rortho(Fc1,sX,Fc2) };

case 'in' %- Fc1 is in list of contrasts Fc2

%=======================================================================

% [iFc2 iFc1] = spm_FcUtil('In', Fc1, sX, Fc2)

% returns indice of Fc2 if "in", 0 otherwise

% NB : If T- stat, the routine checks whether Fc.c is of

% size one. This is ensure if contrast is set

% or manipulated (ortho ..) with spm_FcUtil

% note that the algorithm works \emph{only because} Fc2(?).c

% and Fc1.c are in space(X')

if nargin < 4, error('Insufficient arguments'), end

if nargout > 2, error('Too many output argument.'), end

Fc1 = varargin{2}; Fc2 = varargin{4}; sX = varargin{3};

L1 = length(Fc1);

if ~L1, warning('no contrast given to in');

if nargout == 2, varargout = {[] []};

else varargout = {[]}; end;

return;

end

for i=1:L1

if ~sf_IsFcon(Fc1(i)), error('Fc1(i) must be a contrast'), end

end

L2 = length(Fc2);

if ~L2, error('must have at least a contrast in Fc2'); end

for i=1:L2

if ~sf_IsFcon(Fc2(i)), error('Fc2(i) must be F-contrast'), end

end

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

[idxFc2,idxFc1] = sf_in(Fc1, sX, Fc2);

if isempty(idxFc2), idxFc2 = 0; end

if isempty(idxFc1), idxFc1 = 0; end

switch nargout

case {0,1}

varargout = { idxFc2 };

case 2

varargout = { idxFc2 idxFc1 };

otherwise

error('Too many or not enough output arguments');

end

case '~unique' %- Fc list unique

%=======================================================================

% idx = spm_FcUtil('~unique', Fc, sX)

%- returns indices of redundant contrasts in Fc

%- such that Fc(idx) = [] makes Fc unique.

%- if already unique returns []

if nargin ~= 3, error('Insufficient/too many arguments'), end

Fc = varargin{2}; sX = varargin{3};

%----------------------------

L1 = length(Fc);

if ~L1, warning('no contrast given '); varargout = {[]}; return; end

for i=1:L1

if ~sf_IsFcon(Fc(i)), error('Fc(i) must be a contrast'), end

end

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

%----------------------------

varargout = { unique(sf_notunique(Fc, sX))};

case {'0|[]','[]|0'} %- Fc is null or empty

%=======================================================================

% b = spm_FcUtil('0|[]', Fc, sX)

% returns 1 if F-contrast is empty or null; assumes the contrast is set.

if nargin ~= 3, error('Insufficient/too many arguments'), end

Fc = varargin{2}; sX = varargin{3};

%----------------------------

L1 = length(Fc);

if ~L1, warning('no contrast given to |_'); varargout = {[]}; return; end

for i=1:L1

if ~sf_IsFcon(Fc(i)), error('Fc(i) must be a contrast'), end

end

if ~spm_sp('isspc',sX), sX = spm_sp('set',sX); end;

%----------------------------

idx = [];

for i=1:L1

if sf_isempty_X1o(Fc(i)) || sf_isnull(Fc(i),sX), idx = [idx i]; end

end

if isempty(idx)

varargout = {0};

else

varargout = {idx};

end

%=======================================================================

otherwise

%=======================================================================

error('Unknown action string in spm_FcUtil')

end; %---- switch lower(action),

%=======================================================================

%=======================================================================

% Sub Functions

%=======================================================================

%=======================================================================

%=======================================================================

% Fcon = spm_FcUtil('FconFields')

function Fc = sf_FconFields

Fc = struct(...

'name', '',...

'STAT', '',...

'c', [],...

'X0', struct('ukX0',[]),...

'iX0', [],...

'X1o', struct('ukX1o',[]),...

'eidf', [],...

'Vcon', [],...

'Vspm', []);

%=======================================================================

% used internally. Minimum contrast structure

function minFc = sf_MinFcFields

minFc = struct(...

'name', '',...

'STAT', '',...

'c', [],...

'X0', [],...

'X1o', []);

%=======================================================================

% yes_no = spm_FcUtil('IsFcon',Fc)

function b = sf_IsFcon(Fc)

%- check that minimum fields of a contrast are in Fc

b = 1;

minnames = fieldnames(sf_MinFcFields);

FCnames = fieldnames(Fc);

for str = minnames'

b = b & any(strcmp(str,FCnames));

if ~b, break, end

end

%=======================================================================

% used internally; To be set, a contrast structure should have

% either X1o or X0 non empty. X1o can be non empty because c

% is non empty.

function b = sf_IsSet(Fc)

b = ~sf_isempty_X0(Fc) | ~sf_isempty_X1o(Fc);

%=======================================================================

% used internally

%

function v = sf_ver(Fc)

if isstruct(Fc.X0), v = 2; else v = 1; end

%=======================================================================

% used internally

function b = sf_isempty_X1o(Fc)

if sf_ver(Fc) > 1,

b = isempty(Fc.X1o.ukX1o);

%- consistency check

if b ~= isempty(Fc.c),

Fc.c, Fc.X1o.ukX1o, error('Contrast internally not consistent');

end

else

b = isempty(Fc.X1o);

%- consistency check

if b ~= isempty(Fc.c),

Fc.c, Fc.X1o, error('Contrast internally not consistent');

end

end

%=======================================================================

% used internally

function b = sf_X1o(Fc,sX)

if sf_ver(Fc) > 1,

b = spm_sp('ox',sX)*Fc.X1o.ukX1o;

else

b = Fc.X1o;

end

%=======================================================================

% used internally

function b = sf_X0(Fc,sX)

if sf_ver(Fc) > 1,

b = spm_sp('ox',sX)*Fc.X0.ukX0;

else

b = Fc.X0;

end

%=======================================================================

% used internally

function b = sf_isempty_X0(Fc)

if sf_ver(Fc) > 1,

b = isempty(Fc.X0.ukX0);

else

b = isempty(Fc.X0);

end

%=======================================================================

% Hsqr = spm_Fcutil('Hsqr',Fc,sX)

function hsqr = sf_Hsqr(Fc,sX)

%

% Notations : H equiv to X1o, H = uk*a1, X = uk*ax, r = rk(sX),

% sX.X is (n,p), H is (n,q), a1 is (r,q), ax is (r,p)

% oxa1 is an orthonormal basis for a1, oxa1 is (r,qo<=q)

% Algorithm :

% v1 : Y'*H*(H'*H)-*H'*Y = b'*X'*H*(H'*H)-*H'*X*b

% = b'*X'*oxH*oxH'*X*b

% so hsrq is set to oxH'*X, a (q,n)*(n,p) op. computation of oxH

% v2 : X'*H*(H'*H)-*H'*X = ax'*uk'*uk*a1*(a1'*uk'*uk*a1)-*a1'*uk'*uk*ax

% = ax'*a1*(a1'*a1)-*a1'*ax

% = ax'*oxa1*oxa1'*ax

%

% so hsrq is set to oxa1'*ax : a (qo,r)*(r,p) operation! -:))

% computation of oxa1.

%-**** fprintf('v%d\n',sf_ver(Fc));

if sf_ver(Fc) > 1,

hsqr = spm_sp('ox',spm_sp('set',Fc.X1o.ukX1o))' * spm_sp('cukx',sX);

else

hsqr = spm_sp('ox',spm_sp('set',Fc.X1o))'*spm_sp('x',sX);

end

%=======================================================================

% H = spm_FcUtil('H',Fc)

function H = sf_H(Fc,sX)

%

% Notations : H equiv to X1o, H = uk*a1, X = uk*ax

% Algorithm :

% v1 : Y'*H*(H'*H)-*H'*Y = b'*X'*H*(H'*H)-*H'*X*b

% = b'*c*(H'*H)-*c'*b

% = b'*c*(c'*(X'*X)-*c)-*c'*b

%- v1 : Note that pinv(Fc.X1o' * Fc.X1o) is not too bad

%- because dimensions are only (q,q). See sf_hsqr for notations.

%- Fc.c and Fc.X1o should match. This is ensure by using FcUtil.

%-**** fprintf('v%d\n',sf_ver(Fc));

if sf_ver(Fc) > 1,

hsqr = sf_Hsqr(Fc,sX);

H = hsqr' * hsqr;

else

H = Fc.c * pinv(Fc.X1o' * Fc.X1o) * Fc.c';

% H = {c*spm_sp('x-',spm_sp('Set',c'*spm_sp('xpx-',sX)*c) )*c'}

end

%=======================================================================

% Yc = spm_FcUtil('Yc',Fc,sX,b)

function Yc = sf_Yc(Fc,sX,b)

Yc = sX.X*spm_sp('xpx-',sX)*sf_H(Fc,sX)*b;

%=======================================================================

% Y0 = spm_FcUtil('Y0',Fc,sX,b)

function Y0 = sf_Y0(Fc,sX,b)

Y0 = sX.X*(eye(spm_sp('size',sX,2)) - spm_sp('xpx-',sX)*sf_H(Fc,sX))*b;

%=======================================================================

% Fc = spm_FcUtil('|_',Fc1, sX, Fc2)

function Fc1o = sf_fcortho(Fc1, sX, Fc2)

%--- use the space facility to ensure the proper tolerance dealing...

c1_2 = Fc1.c - sf_H(Fc2,sX)*spm_sp('xpx-:',sX,Fc1.c);

Fc1o = spm_FcUtil('Set',['(' Fc1.name ' |_ (' Fc2.name '))'], ...

Fc1.STAT, 'c ',c1_2,sX);

%- In the large (scans) dimension :

%- c = sX.X'*spm_sp('r:',spm_sp('set',Fc2.X1o),Fc1.X1o);

%- Fc1o = spm_FcUtil('Set',['(' Fc1.name ' |_ (' Fc2.name '))'], ...

%- Fc1.STAT, 'c',c,sX);

%=======================================================================

function b = sf_Rortho(Fc1,sX,Fc2)

if isempty(Fc2)

if length(Fc1) <= 1, b = 0;

else

c1 = cat(2,Fc1(:).c);

b = ~any(any( abs(triu(c1'*spm_sp('xpx-:',sX,c1), 1)) > sX.tol));

end

else

c1 = cat(2,Fc1(:).c); c2 = cat(2,Fc2(:).c);

b = ~any(any( abs(c1'*spm_sp('xpx-:',sX,c2)) > sX.tol ));

end

%=======================================================================

% b = spm_FcUtil('0|[]', Fc, sX)

%- returns 1 if F-contrast is empty or null; assumes the contrast is set.

%- Assumes that if Fc.c contains only zeros, so does Fc.X1o.

%- this is ensured if spm_FcUtil is used

function boul = sf_isnull(Fc,sX)

%

boul = ~any(any(spm_sp('oPp:',sX,Fc.c)));

%=======================================================================

% Fc = spm_FcUtil('Set',name, STAT, set_action, value, sX)

function boul = sf_is_T(sX,c)

%- assumes that the dimensions are OK

%- assumes c is not empty

%- Does NOT assumes that c is space of sX'

%- A rank of zero can be defined

%- if the rank == 1, checks whether same directions

boul = 1;

if ~spm_sp('isinspp',sX,c), c = spm_sp('oPp:',sX,c); end;

if rank(c) > 1 || any(any(c'*c < 0)), boul = 0; end;

%=======================================================================

function [idxFc2, idxFc1] = sf_in(Fc1, sX, Fc2)

L2 = length(Fc2);

L1 = length(Fc1);

idxFc1 = []; idxFc2 = [];

for j=1:L1

%- project Fc1(j).c if not estimable

if ~spm_sp('isinspp',sX,Fc1(j).c), %- warning ?

c1 = spm_sp('oPp:',sX,Fc1(j).c);

else

c1 = Fc1(j).c;

end

sc1 = spm_sp('Set',c1);

S = Fc1(j).STAT;

boul = 0;

for i =1:L2

if Fc2(i).STAT == S

%- the same statistics. else just go on to the next contrast

boul = spm_sp('==',sc1,spm_sp('oPp',sX,Fc2(i).c));

%- if they are the same space and T stat (same direction),

%- then check whether they are in the same ORIENTATION

%- works because size(X1o,2) == 1, else .* with (Fc1(j).c'*Fc2(i).c)

if boul && S == 'T'

atmp = sf_X1o(Fc1(j),sX);

btmp = sf_X1o(Fc2(i),sX);

boul = ~any(any( (atmp' * btmp) < 0 ));

end

%- note the indices

if boul, idxFc1 = [idxFc1 j]; idxFc2 = [idxFc2 i]; end

end

end

end %- for j=1:L1

%=======================================================================

function idx = sf_notunique(Fc, sX)

%- works recursively ...

%- and use the fact that [] i == []

%- quite long for large sets ...

l = length(Fc);

if l == 1, idx = [];

else

idx = [ (1 sf_in(Fc(1),sX,Fc(2:l))) (1 sf_notunique(Fc(2:l), sX))];

end