tic;
%warning off
[rows,cols]=size(x1);
fvalues = zeros(cols,1);
pvalues = zeros(cols,1);
r2values = zeros(cols,1);
adjr2values = zeros(cols,1);
cpvalues = zeros(cols,1);
aicvalues = zeros(cols,1);
bicvalues = zeros(cols,1);
rssvalues = zeros(cols,1);
selvars = zeros(cols,cols);
summary = zeros(cols,8);
for ijk = 1:15  %1:cols-1;
    c = ijk;
    [sel, b, rss, r2, cp, rss2, f, adj, pvalue, aic, bic] = bestsub(x1, y, c);
    fvalues(ijk) = f; pvalues(ijk) = pvalue; aicvalues(ijk) = aic;
    r2values(ijk) = r2; cpvalues(ijk) = cp; bicvalues(ijk) = bic;
    adjr2values(ijk) = adj; rssvalues(ijk) = rss;
    for iii = 1:c
        selvars(c,iii) = sel(iii);
    end
end
[m,d]=size(x1);
mu = mean(y); ydiff = (y-mu).*(y-mu); sstot = sum(ydiff);
rss = rss2 .* (m-d-1);
rssvalues(cols) = rss;
cpvalues(cols,1) = rss./rss2 - m + 2*(d+1);
cc = d+1; mse = rss./(m-cc); msr = (sstot-rss)./d;
fvalues(cols) = msr./mse; v1 = d; v2 = m - cc; X = v2./(v2+v1.*f);
pvalues(cols) = betainc(X,v2./2,v1./2);
LL = -(m./2).*log(2.*pi.*(rss./m))-m./2;
aicvalues(cols) = -2.*(LL-cc); bicvalues(cols) = LL-(1./2).*cc.*log(m);
r2values(cols) = (sstot-rss)./sstot;
adjr2values(cols) = 1 - (1-r2).*((m-1)./(m-cc));
summary(:,1) = 1:cols; summary(:,2) = r2values; summary(:,3) = adjr2values;
summary(:,4) = fvalues; summary(:,5) = pvalues; summary(:,6) = cpvalues;
summary(:,7) = aicvalues; summary(:,8) = bicvalues;
toc