function [out] = PPIRunRegression(regs,timec,settings)
% runs for each subject, for each time point a separate regression with all
% regressors (regs)

% some lengths extracted here to make code below more compact
ntimep = size(timec.Region(1).v{1},2);   % number of timepoints in each epoch
%nregion= size(timec.Region,2);           % number of regions
nregres= length(regs.Labels);              % number of regressors

% get the index for the target region
targRegion = strcmp(settings.targetRegion,{timec.Region.name});
seedRegion = strcmp(settings.seedRegion,{timec.Region.name});
globalTC   = strcmp(settings.globalTC,{timec.Region.name});
psychInd   = strcmp(settings.psychReg,regs.Labels);
 
% % remove the contrasts if they exist
% if isfield(regs,'contrast') 
%     regs(1).l=regs(1).l(1:size(regs(1).v,2))
% end
extraRegs = sum(globalTC) + 2; % how many extra regressors to do the ppi
% run the regressions
%if numel(unique(regs(1).v(1:80,psychInd)))==3
%    betas = nan(25,1,ntimep,nregres+1+ extraRegs+1); %subs, timepoints, regressors
%else
    betas = nan( length(regs.DesignMatrix),1,ntimep,nregres+1+ extraRegs);
    betas2 = nan(length(regs.DesignMatrix) ,1,ntimep,nregres+1+ extraRegs);
%end
%partCorr = nan(25,1,ntimep,nregres+1+ extraRegs);% having two for final dimension to keep matrix size consistent
for is = 1:length(regs.DesignMatrix)%[1:18 20:25] % for each participant
    for ip = 1:ntimep % for each time point
       % keyboard
        ydat  = timec.Region(targRegion).v{is}(:,ip);   % ydat for the regression, i.e. the intrapolated BOLD values
        indtr = size(ydat,1); % how many fmri trials we have for each person
        
        % normalise the predictors:
        temp=regs.DesignMatrix{is}(1:indtr,:);
        preds=  zeros(size(regs.DesignMatrix{is}(1:indtr,:))); % Predictors that we'll use in the glm
         
        for icol=1:size(preds,2)
            preds(~isnan(temp(:,icol)),icol) = zscore(temp(~isnan(temp(:,icol)),icol));
        end
        % add the time course of interest, the psy x phys interaction
        % and the global time course as predictors
        phys = timec.Region(seedRegion).v{is}(:,ip);
        physXpsy = zscore(zscore(phys).* preds(:,psychInd));

        % instead of main effect and interaction: split the physiological
        % by psychological - first check whether we have exactly three
        % conditions, if so, split data into 3
 
            physSplit = zeros(length(phys),2);
            psychSplitInd = [preds(:,psychInd)> median(preds(:,psychInd)) preds(:,psychInd)<median(preds(:,psychInd))];
            physSplit(psychSplitInd(:,1),1) = zscore(phys(psychSplitInd(:,1)));
            physSplit(psychSplitInd(:,2),2) = zscore(phys(psychSplitInd(:,2)));
         
        globalControl = zscore(timec.Region(globalTC).v{is}(:,ip)); 
        preds2 = [preds globalControl phys physXpsy];
        %preds3 = [preds globalControl physSplit];
        %betas(is,1,ip,:) = glmfit(preds3,ydat);
        betas2(is,1,ip,:) = glmfit(preds2,ydat); %running fitglm instead would give more information, but it takes longer to run (and we only need the regression weights
       % partCorr(is,1,ip,end) = partialcorr(ydat,physXpsy,[preds globalControl phys]);
    end
end
%keyboard
%out.betas = betas;  %partCorr;% betas;
out.betas2= betas2;
if ~isempty(globalTC)
    %if exist('uniques','var')
    %     out.labels= ['Const' regs(1).l 'GlobalTC' [settings.seedRegion 'xHigh' regs(1).l{psychInd}]  [settings.seedRegion 'xMid' regs(1).l{psychInd}] [settings.seedRegion 'xLow' regs(1).l{psychInd}]];
    %else
        %out.labels= ['Const' regs.Labels{:} 'GlobalTC' [settings.seedRegion 'xHigh' regs.Labels{psychInd}] [settings.seedRegion 'xLow' regs.Labels{psychInd}]];%settings.seedRegion [settings.seedRegion ' x ' settings.psychReg]];
        %out.labels2= ['Const' regs.Labels 'GlobalTC' [settings.seedRegion] [settings.seedRegion 'x' regs.Labels{psychInd}]];
        out.labels{1}='Const';
        for r=1:length(regs.Labels)
            out.labels{1+r}=regs.Labels{r};
        end
        out.labels{1+length(regs.Labels)+1}='GlobalTC';
        out.labels2=out.labels;
        out.labels{1+length(regs.Labels)+2}=[settings.seedRegion 'xHigh' regs.Labels{psychInd}];
        out.labels{1+length(regs.Labels)+3}=[settings.seedRegion 'xLow' regs.Labels{psychInd}];
        out.labels2{1+length(regs.Labels)+2}=[settings.seedRegion];
        out.labels2{1+length(regs.Labels)+3}=[settings.seedRegion 'x' regs.Labels{psychInd}];



    %end
else
    out.labels= ['Const' regs.Labels settings.seedRegion [settings.seedRegion ' x ' settings.psychReg]];
    keyboard
end
% if strcmp(settings.psychReg,'occludedvssharedMainDiff')
%    keyboard 
% end

out.descr = 'subject x region x time point x coefficients';
out.regions= {timec.Region.name};
out.ppiInfo = settings;
