function cvestr = fastcv(y,t,par,K,s,o)

%function cvestr = fastcv(x,t,par,K,rs,o)
%
%Fast cross-validation prediction errors for Self-Modeling Registration
%
%INPUT:
%   x:    data matrix (n x m) (subjects in rows)
%   t:    time grid (1 x m)
%   par:  tentative parameters (r x 2)
%         Each row corresponds to a duplet (q,p)
%   K:    number of test sets (1 x nk; K = [n 5 10] is suggested)
%   rs:   number of random starts (scalar; optional argument, default = 20)
%   o:    order of splines (scalar; optional argument, default = 3)
%
%OUTPUT:
%   cvestr: struct with the following fields:
%   	cve: individual cross-validation prediction errors (n x r x nk)
%            To get the total CVE, take mean over first array dimension.
%       K:   same as input
%       par: same as input
%
%External subrutines called: smreg.m, bspl.m

if nargin<4
    disp('Not enough input arguments')
    cvestr = [];
    return
end
if nargin<6
    o = 3;
end
if nargin<5
    s = 20;
end

[n,m] = size(y);
if size(t,1)>1 & size(t,1)==m
    t = t';
end
if size(t,2)~=m
    disp('Dimensions of X and T do not match')
    cvestr = [];
    return
end

if size(par,2)~=2
    disp('Rows of PAR must be duplets')
    cvestr = [];
    return
end

nk = length(K);
for i = 1:nk
    if K(i)>n
        disp('K must be less than or equal to N')
        cvestr = [];
    	return
    elseif K(i)==1
        disp('Are you sure you want this? Leave-one-out CV corresponds to K=N, not K=1')
    end
end

% Prealocation
r = size(par,1);
h1 = ones(n,m);
g = t;
mu = zeros(1,m);
mug = zeros(1,m);
tt = [t(1),(t(1:m-1)+t(2:m))/2,t(m)];
tdif = tt(2:m+1)-tt(1:m);
cve = zeros(n,r,nk);

% Computing CVE for each (q,p)
for i = 1:r
    disp([sprintf('\n') '----> Fitting full model: parameter ' num2str(i)])
    rand('state',0)
    regstr = smreg(y,t,par(i,1),par(i,2),s,o);
    knots = linspace(t(1),t(m),par(i,2)-o+2);
    b1 = bspl(t,o,knots,0)';
    h1 = ones(n,m) + regstr.s*regstr.C*b1;
    for l = 1:nk
	% Defining test samples
	n1 = round(n/K(l));
	for j = 1:K(l)-1
	    isub{j} = n1*(j-1)+1:n1*j;
	end
	isub{K(l)} = n1*(K(l)-1)+1:n;
	for j = 1:K(l)
	    mu = regstr.a(isub{j})'*(regstr.xw(isub{j},:).*h1(isub{j},:))...
				./(regstr.a(isub{j}).^2'*h1(isub{j},:));
	    for ii = 1:length(isub{j})
	        g(2:m-1) = interp1q(regstr.w(isub{j}(ii),:)',t',t(2:m-1)')';
		mug = interp1q(t',mu',g')';
		cve(isub{j}(ii),i,l) = (y(isub{j}(ii),:)-regstr.a(isub{j}(ii)).*mug).^2*tdif';
	    end
	end
    end
end

cvestr = struct('cve',cve,'K',K,'par',par);