function [idxSelected,count]=glmRoutine(subsetX,Y,n,selecCriterion,intercept)

count=0;
if selecCriterion==0
    if intercept==0
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial','constant','off');
        pVals=stats.p;
    else
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial');
        pVals=stats.p;
        pVals(1,:)=[];
       if all(pVals>0.99)|all(pVals==0)
           count=count+1;
           sTmp=size(subsetX,2);
           sB1=floor(sTmp/2);
           block1=subsetX(:,1:sB1);
           block2=subsetX(:,sB1+1:end);
           [~,~,stats]=glmfit(block1,[Y,ones(n,1)],'binomial');
           pVals1=stats.p;
           pVals1(1,:)=[];
           if all(pVals1>0.99)
               pVals1=0.0001*ones(size(pVals1));
           end
           [~,~,stats]=glmfit(block2,[Y,ones(n,1)],'binomial');
           pVals2=stats.p;
           pVals2(1,:)=[];
           if all(pVals2>0.99)
               pVals2=0.0001*ones(size(pVals2));
           end
           pVals=[pVals1;pVals2];
       end
    end
    idxSelected=find(pVals<=0.005);
elseif selecCriterion==1
    if intercept==0
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial','constant','off');
        pVals=stats.p;
    else
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial');
        pVals=stats.p;
        pVals(1,:)=[];
        if all(pVals>0.99)|all(pVals==0)
                       count=count+1;
            sTmp=size(subsetX,2);
            sB1=floor(sTmp/2);
            block1=subsetX(:,1:sB1);
            block2=subsetX(:,sB1+1:end);
            [~,~,stats]=glmfit(block1,[Y,ones(n,1)],'binomial');
            pVals1=stats.p;
            pVals1(1,:)=[];
            if all(pVals1>0.99)
                pVals1=0.0001*ones(size(pVals1));
            end
            [~,~,stats]=glmfit(block2,[Y,ones(n,1)],'binomial');
            pVals2=stats.p;
            pVals2(1,:)=[];
            if all(pVals2>0.99)
                pVals2=0.0001*ones(size(pVals2));
            end
            pVals=[pVals1;pVals2];
        end
    end
    idxSelected=find(pVals<=0.01);
elseif selecCriterion==-1
    if intercept==0
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial','constant','off');
        pVals=stats.p;
    else
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial');
        pVals=stats.p;
        pVals(1,:)=[];
        if all(pVals>0.99)|all(pVals==0)
                       count=count+1;
            sTmp=size(subsetX,2);
            sB1=floor(sTmp/2);
            block1=subsetX(:,1:sB1);
            block2=subsetX(:,sB1+1:end);
            [~,~,stats]=glmfit(block1,[Y,ones(n,1)],'binomial');
            pVals1=stats.p;
            pVals1(1,:)=[];
            if all(pVals1>0.99)
                pVals1=0.0001*ones(size(pVals1));
            end
            [~,~,stats]=glmfit(block2,[Y,ones(n,1)],'binomial');
            pVals2=stats.p;
            pVals2(1,:)=[];
            if all(pVals2>0.99)
                pVals2=0.0001*ones(size(pVals2));
            end
            pVals=[pVals1;pVals2];
        end
    end
    idxSelected=find(pVals<=0.05);
else
    if intercept==0
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial','constant','off');
        pVals=stats.p;
    else
        [~,~,stats]=glmfit(subsetX,[Y,ones(n,1)],'binomial');
        pVals=stats.p;
        pVals(1,:)=[];
        if all(pVals>0.99)|all(pVals==0)
                       count=count+1;
            sTmp=size(subsetX,2);
            sB1=floor(sTmp/2);
            block1=subsetX(:,1:sB1);
            block2=subsetX(:,sB1+1:end);
            [~,~,stats]=glmfit(block1,[Y,ones(n,1)],'binomial');
            pVals1=stats.p;
            pVals1(1,:)=[];
            if all(pVals1>0.99)
                pVals1=0.0001*ones(size(pVals1));
            end
            [~,~,stats]=glmfit(block2,[Y,ones(n,1)],'binomial');
            pVals2=stats.p;
            pVals2(1,:)=[];
            if all(pVals2>0.99)
                pVals2=0.0001*ones(size(pVals2));
            end
            pVals=[pVals1;pVals2];
        end
    end
    [~,sortIndex]=sort(pVals,'ascend');
    if length(sortIndex)>1
        idxSelected=sortIndex(1:2);
    else
        idxSelected=sortIndex;
    end
end

