% compute CP statistic for GLARS estimates
% need all steps k=1,..,J (J = number of groups)
% p=[p1,..,pJ] groups dimensions
% Theta is a matrix: the k-th column contains estimates for the k-th step

function CPindex=CP_GLARS(y,X,Theta,p);
n = length(y);

J=length(p);

e=y-X*Theta(:,J);
sigma=std(e);


for j=1:J,
    Den(j) = norm(Theta(sum(p(1:j-1))+1:sum(p(1:j)),1));
    for l=1:J-1,
        Den(j) = Den(j) + norm(Theta(sum(p(1:j-1))+1:sum(p(1:j)),l+1)-Theta(sum(p(1:j-1))+1:sum(p(1:j)),l));
    end;
end;

for k=1:J,
    df(k)=k;
    for j=1:J,
        Num(j) = norm(Theta(sum(p(1:j-1))+1:sum(p(1:j)),1));
        for l=1:k-1
            Num(j) = Num(j) + norm(Theta(sum(p(1:j-1))+1:sum(p(1:j)),l+1)-Theta(sum(p(1:j-1))+1:sum(p(1:j)),l));
        end;
        Num(j)=Num(j)*(p(j)-1);
        df(k) = df(k) + Num(j)/Den(j);
    end;

    CPindex(k) = norm(y-X*Theta(:,k))^2/(n*sigma^2) - 1 + 2*df(k)/n;
end;
