function [B,BINT,R,RINT,STATS] = regress2(y,X,logic)
% function [B,BINT,R,RINT,STATS] = regress2(y,X,logic)
%
% Simple extension of Matlab's "regress" to N-Dimensional y.  Does a
% for-loop over additional dimensions of y, excluding nans.
%
% y:        observations
% x:        predictors
% logic:    optional. if string "logic", converts BINT and RINT to logicals
%
% BINT:     1: c.i. of regression coefficient crosses zero. Coefficients
%           are not supposed to cross zero.
% RINT:     1: Possible outlier.
%
% use data(BINT | RINT) = nan; to reject
%
% STATS:      r^2, F, p, error variance (should be rms residual)

dim = size(y);

if length(dim) > 3, error('add reshaping of N-D'), end

B = nan(dim(2),size(X,2)); % regression coefficients; e.g. slope,offset
if nargout > 4
    STATS = nan(dim(2),4); % r^2, F, p, error variance (should be rms residual)
elseif nargout > 3
    RINT = nan([dim 2]) ;% outliers
elseif nargout > 2
    R = nan(dim); % residuals
elseif nargout > 1
    BINT = nan(dim(2),size(X,2),2); % coefficient intervals
end

% this is almost x6 reduction in regression calls:
do = sum(~isnan(y)) > 2;% need at least 2 values for a linear regression.  Otherwise errors and nans
j = find(do);
J = length(j);

warning('off','MATLAB:divideByZero')
warning('off','stats:regress:RankDefDesignMat')
warning('off','stats:regress:NoConst')

for n = j
%     disp(sprintf('regressing: %f%%',100*find(j==n)/J))
    
    if nargout > 4
        [B(n,:),BINT(n,:,:),R(:,n),RINT(:,n,:),STATS(n,:)] = regress(   y(:,n) ,  X  );
    elseif nargout > 3
        [B(n,:),BINT(n,:,:),R(:,n),RINT(:,n,:)] = regress(   y(:,n) ,  X  );
    elseif nargout > 2
        [B(n,:),BINT(n,:,:),R(:,n)] = regress(   y(:,n) ,  X  );
    elseif nargout > 1
        [B(n,:),BINT(n,:,:)] = regress(   y(:,n) ,  X  );
    else
        B(n,:) = regress(   y(:,n) ,  X  );
    end
%     chk = [bint(:); stats(:)];
%     if any(isnan(chk) | isinf(chk)), error('wtf'), end
end

warning('on','MATLAB:divideByZero')
warning('on','stats:regress:RankDefDesignMat')
warning('on','stats:regress:NoConst')

if nargin > 2 && strcmp(logic,'logic')
    if nargout > 3
        RINT = RINT(:,:,1) > 0 | 0 > RINT(:,:,2);
    end
    if nargout > 1
        BINT = BINT(:,:,1) <= 0 & 0 <= BINT(:,:,2);
    end
end