
function varargout = mm_model(varargin);
%-  Format varargout = mm_model(varargin);
%- Set sub-space of interest and the related matrix of normalisation. 
%- See subfunction for details.
%====================================================================
% This package was developped by Ferath Kherif primarily
% with some help from Jean-Baptiste Poline.
% 
% Please do not distribute and redirect requests for the toolbox to 
% us. For bugs, remarks, additions, etc, please contact jb at 
% poline@shfj.cea.fr. Some of the ideas implemented in the MM toolbox
% have not yet been published. Please refrain from publishing the  
% original ideas that can be found here ... If you are using this  
% material for publication, please see with us how you can 
% acknowledge our work.
%====================================================================
% WARNING : THIS SOFTWARE IS DISTRIBUTED FREE WITHOUT ANY GUARANTY ! 
%====================================================================


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

%$Author: jcooper $
%$Id: mars_mm_model.m,v 1.1.1.1 2002/10/04 17:40:09 jcooper Exp $
%$ProjectHeader: MLM 0.39 Thu, 22 Mar 2001 11:18:01 +0100 poline $ 

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


[NF,nu,h,d,M12,XG,sXG] = sf_model_mlm(varargin{2:nargin});
varargout 	= {NF,nu,h,d,M12,XG,sXG};
		


	
%===================================================================
function [NF,nu,h,d,M12,XG,sXG] = sf_model_mlm(Xs, V, nROI, xC, erdf);
% Set sub-space of interest and the related matrix of normalisation. 
% FORMAT [NF,nu,h,d,M12,XG] = mm_model();
%- nu, h, d : degrees of freedom
%- NF : matrix of normalisation
%===================================================================


%--------------------------------------------------------------------
%- SET, COMPUTE,NORMALIZE SPACES OF INTEREST
%--------------------------------------------------------------------
%- set X10 and XG
%- XG= X -PG(X), PG projection operator on XG (cf. eq 1, 2)
%--------------------------------------------------------------------
sX1o	= spm_sp('set',spm_FcUtil('X1o',xC,Xs));
sXG	= spm_sp('set',spm_FcUtil('X0',xC,Xs));
X1o 	= spm_sp('oP',sX1o,Xs.X);
XG  	= spm_sp('r',sXG,Xs.X);

%- Compute Normalized effexts : M1/2=X'G*V*XG (cf eq 3)
%--------------------------------------------------------------------
% warning off;
up	= spm_sp('uk',sX1o); ; %- PG=up*up'
qi	= up'*Xs.X;
sigma	= up'*V*up;
M12	= (chol(sigma)*qi)';
M_12	= pinv(M12);

%- Compute NF : normalise factor (cf eq 4)
%--------------------------------------------------------------------
NF	= M_12*spm_sp('X',Xs)'*spm_sp('r',sXG,spm_sp('X',Xs));

%- degrees of freedom
%- nROI : number of ROI (corresponds to the number of Resels) 
%--------------------------------------------------------------------
d	= nROI*(4*log(2)/pi)^(3/2);
h	= sX1o.rk; %-rank of the sub-space of interest.
nu	= erdf;
