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);
rperm = zeros(1,cols); X = x1(:,1:dd); x1 = []; smark = zeros(1,dd);
vsel = []; vset = 1:dd;
for h = 1:dd
    ssbest = -999999999999;
    for i = 1:dd
        if smark(i) == 1
            continue
        end
        vtrial = setdiff(vset,union(vsel,i)); vlen = length(vtrial);
        x = [ones(size(y))]; icol = 1;
        for j = 1:vlen
            icol = icol + 1; vv = vtrial(j);
            x = [x X(:,vv)];
        end
        for j = 1:vlen
            icol = icol + 1; vv = vtrial(j);
            x = [x X(:,vv).^2];
        end
        for j = 1:vlen-1
            for k = 1+j:vlen
                icol = icol + 1; vv = vtrial(j); uu = vtrial(k);
                x = [x X(:,vv).*X(:,uu)];     
            end
        end
        v = x\y; yhat = x*v; R = sum((y-yhat).^2);
        if R > ssbest
            ssbest = R; isel = i;
        end
    end
    vsel(h) = isel; smark(isel) = 1;
end
C = vsel
x1 = X(:,C);
for i = 1:dd
    x1(:,i+dd) = x1(:,i).^2;
    rperm(i) = C(i);
    rperm(i+dd) = C(i) + dd;
end
nn = 2.*dd;
for i = 1:dd-1
    for j = i+1:dd
        nn = nn+1;
        x1(:,nn) = x1(:,i).*x1(:,j);
        nnn = 2.*dd;
        for ii = 1:dd-1
            for jj = ii+1:dd
                nnn = nnn+1;
                if (ii==C(i)&&jj==C(j)) || (ii==C(j) && jj == C(i))
                    rperm(nn) = nnn;
                end
            end
        end
    end
end

r2best = 0;
for ijk = 1:15 %cols-1;
    c = ijk
    [sel, b, rss, r2, cp, rss2, f, adj, pvalue, aic, bic] = bestsub_h(x1, y, c, dd, r2best);
    fvalues(ijk) = f; pvalues(ijk) = pvalue; aicvalues(ijk) = aic;
    r2values(ijk) = r2; cpvalues(ijk) = cp; bicvalues(ijk) = bic;
    adjr2values(ijk) = adj; rssvalues(ijk) = rss;
    sel = rperm(sel); sel = sort(sel);
    for iii = 1:c
        selvars(c,iii) = sel(iii);
    end
    r2best = r2;
end
[m,d]=size(x1);
mu = mean(y); ydiff = (y-mu).*(y-mu); sstot = sum(ydiff);
rss = rss2 .* (m-d-1); 
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