%----------------------------------------------------------------------------------
%- set the methods

COV_estim = 'assumed';           % covariance is assumed to be imposed by filter K
GLM_resol = 'OLS'                % ordinary least square 
 


%----------------------------------------------------------------------------------
%- get the design structure, the filter and the images

load(spm_get(1,'SPMcfg.mat','Select SPMcfg.mat'));
if ~isfield(xX,'X'), error('Please load an SPMcfg containing an X'), end
Xs = spm_sp('set',xX.X);

%----------------------------------------------------------------------------------
%- get the ROIs and the data : Y

load Y

nROI = size(Y,2);  %- Y is a time by nROI matrix


%----------------------------------------------------------------------------------
%- Estimation of the covariance structure of the Ys

switch COV_estim

	case {'AR(p)'}
		%- compute the temporal cov of Y (V) with AR(p)
		if ~isfield(xX,'xVi')
		   xX.xVi = struct(	'Vi', speye(size(xX.X,1)),...
					'Form',	'AR(p)'); 
		end
		% xX.xVi = estimate_cov(Y,xX);


	case {'assumed'}
		if ~isfield(xX,'xVi')
		   xX.xVi = struct(	'Vi', speye(size(xX.X,1)),...
					'Form',	'none'); 
		end
		%- else, the covariance structure is supposed to be
		%- stored in xX.xVi

	otherwise
		warning('COV_estim does not exist');
end


switch GLM_resol

	case {'OLS'}
		'OLS'
		%- no filter already defined 
		if ~isfield(xX,'K')
		   xX.K  = speye(size(xX.X,1));
		end
		% else assume that the filter is xX.K
	
		KVi    = spm_filter('apply', xX.K, xX.xVi.Vi); 
		V      = spm_filter('apply', xX.K, KVi'); 
		Y      = spm_filter('apply', xX.K, Y);
		clear KVi;

	case {'MaxLik'} 
		'MaxLik' 
		%- compute the inverse filter -  put it in K ?
		%- filter data and design
		%- V = speye(size(xX.X,1));

	otherwise
		warning('GLM_resol does not exist');
end



%- compute GLM 

[betas RMS erdf] = mars_GLM_estim(Xs, Y, V); 


%- compute stat
%- get xCon, then get some contrasts



PxCon = spm_get(1,'xCon.mat','Select xCon.mat');
load(PxCon);

% [Ic,xCon] = spm_conman(xX,xCon,'T|F',Inf,' Select contrasts ',' ',1);
Ic = 1:length(xCon);

%- results

[StatResu, StatResuP, StatResuPc, MVResu] = mars_stat_compute(xCon(Ic), Xs, V, betas, RMS);
MVres = mars_stat_computeMV(xCon, Xs, V, betas, RMS, Y);

spmResDir = spm_str_manip(PxCon,'H');

save MarsBarStat MVres StatResu StatResuP StatResuPc MVResu betas RMS V Ic spmResDir


