function [infoTransf_maps,infoTransf_mapsSh] = compute_deltaDI_all_trials(neural,behav,sigTime,roiX,roiY,tmproiX,tmproiY,validXlab,validYlab,params,opts,band_label)
% Delta TE comput. across all trials

% Notice: band_label is not really useful here, I just left it to simplify moving
% to across-frequencies computations

% Structure initialization

for featIdx = 1:numel(params.selected_features)
    featLab = params.selected_features{featIdx};
    infoTransf_maps.(featLab).DI = cell(2,1); % this will contain values of DI at fixed stimulus
    infoTransf_maps.(featLab).deltaDI = [];
    infoTransf_mapsSh.(featLab).DI = cell(2,1);
    infoTransf_mapsSh.(featLab).deltaDI = [];
end


for sessIdx = 1:params.nSess
    
    % Behavioral features
    if strcmp(params.stim_feature,'average')
        stim = (mean(behav.v_cons{sessIdx},2)>0.5)';
    elseif strcmp(params.stim_feature,'first_sample')        
        stim = eqpop(behav.v_cons{sessIdx}(:,1), opts.Sbins)';
    end
    choice = behav.v_choice{sessIdx};

    s_vals = unique(stim);
    c_vals = unique(choice);
    assert(numel(s_vals)==numel(c_vals),'Stimulus and choice should have the same number of values')
    shuffIdxs_S = nan(opts.nShuff,numel(stim));
    shuffIdxs_C = nan(opts.nShuff,numel(choice));
    shuffIdxs_all = nan(opts.nShuff,numel(stim));
    for sIdx = 1:opts.nShuff
        ridx = randperm(numel(stim));
        shuffIdxs_all(sIdx,:) = ridx;
        for Ss = 1:numel(s_vals)
            % Stimulus shuff
            idx = find((stim == s_vals(Ss))==1);
            ridx = randperm(numel(idx));
            shuffIdxs_S(sIdx,idx) = idx(ridx);
            % Choice shuff
            idx = find((choice == c_vals(Ss))==1);
            ridx = randperm(numel(idx));
            shuffIdxs_C(sIdx,idx) = idx(ridx);
        end
    end

    for Ss = 1:numel(s_vals) % loop over stimulus values
        bX = []; bY = []; bhX = []; bhY = []; bXt = []; bYt = [];
        
        stim_idxs = find((stim == s_vals(Ss))==1);
        % Binning neural responses
        for t = 1:params.totTimePoints
            edgs = eqpop(neural.data_subj{sessIdx}(:,t,(2*roiX-1)), opts.n_bins);
            [~,bX.(band_label).left(:,t)] = histc(neural.data_subj{sessIdx}(:,t,(2*roiX-1)), edgs);
            edgs = eqpop(neural.data_subj{sessIdx}(:,t,(2*roiY-1)), opts.n_bins);
            [~,bY.(band_label).left(:,t)] = histc(neural.data_subj{sessIdx}(:,t,(2*roiY-1)), edgs);
            edgs = eqpop(neural.data_subj{sessIdx}(:,t,(2*roiX)), opts.n_bins);
            [~,bX.(band_label).right(:,t)] = histc(neural.data_subj{sessIdx}(:,t,(2*roiX)), edgs);
            edgs = eqpop(neural.data_subj{sessIdx}(:,t,(2*roiY)), opts.n_bins);
            [~,bY.(band_label).right(:,t)] = histc(neural.data_subj{sessIdx}(:,t,(2*roiY)), edgs);
        end
    
        %% Compute DI and FIT
                    
        emitFreq = band_label;
        recFreq = band_label;
    
        for featIdx = 1:numel(params.selected_features)
            featLab = params.selected_features{featIdx};
    
            % Time is w.r.t. receiver Y
            for tidx = 1:numel(sigTime.(recFreq).(validYlab))
                t = sigTime.(recFreq).(validYlab)(tidx);
                % For each time point we select a range of delays so that
                % only the possible infomration-emitting time steps are
                % considered, this increases the computational efficiency
                validDelays = t-sigTime.(emitFreq).(validXlab); 
                validDelays(validDelays<1)=[]; % remove negative delays (emitter in the future)
                validDelays(validDelays>params.maxDelay)=[]; % remove delays longer than params.max_delay
                validDelays = sort(validDelays);
    
                if ~isempty(validDelays) % check that some valid delays exist
                    for d = validDelays
                        
                        % Define hX (h = history --> past), hY, Xt and Yt

                        % Define optim. features
                        % Y present
                        bYt.(recFreq).(featLab) = bY.(recFreq).(featLab)(:,t)';

                        % X past
                        bhX.(emitFreq).(featLab) = bX.(emitFreq).(featLab)(:,t-d)';

                        % Y past
                        bhY.(recFreq).(featLab) = bY.(recFreq).(featLab)(:,t-d)';

                        
                        % Compute DI at fixed stim
                        
                        % DI X to Y
                        [tmpDI] = compute_FIT_TE(stim, bhX.(emitFreq).(featLab), bYt.(recFreq).(featLab), bhY.(recFreq).(featLab));
                        infoTransf_maps.(featLab).DI{Ss}(sessIdx,tidx,d) = tmpDI;
    
                        
                        % Conditional shuffling null hypothesis (for FIT)
                        % and simple shuffling for DI
                        if opts.doShuff
                            XSh_all = zeros(1,numel(bhX.(emitFreq).(featLab))); % shuffled across all trials
                            XSh_s = zeros(1,numel(bhX.(emitFreq).(featLab))); % shuffled at fixed S
                            XSh_c = zeros(1,numel(bhX.(emitFreq).(featLab))); % shuffled at fixed C
                            for sIdx = 1:opts.nShuff
    
                                [~,relatSh] = sort(shuffIdxs_all(sIdx,stim_idxs));% shuff idxs relative to current stim value
                                % All trials shuff
                                XSh_all = bhX.(emitFreq).(featLab)(relatSh);
                                % Stimulus fixed shuff
                                [~,relatSh] = sort(shuffIdxs_all(sIdx,stim_idxs));% shuff idxs relative to current stim value
                                XSh_s = bhX.(emitFreq).(featLab)(relatSh);
                                % Choice fixed shuff
                                [~,relatSh] = sort(shuffIdxs_all(sIdx,stim_idxs));% shuff idxs relative to current stim value
                                XSh_c = bhX.(emitFreq).(featLab)(relatSh);
                            
                                if strcmp(params.null_type,'conditioned')
                                    [tmpDI] = compute_FIT_TE(XSh_all, bhX.(emitFreq).(featLab), bYt.(recFreq).(featLab), bhY.(recFreq).(featLab));
                                    infoTransf_mapsSh.(featLab).DI{Ss}(sessIdx,tidx,d) = tmpDI;
                                elseif strcmp(params.null_type,'simple')
                                    [tmpDI] = compute_FIT_TE(XSh_all, bhX.(emitFreq).(featLab), bYt.(recFreq).(featLab), bhY.(recFreq).(featLab));
                                    infoTransf_mapsSh.(featLab).DI{Ss}(sessIdx,tidx,d) = tmpDI;
                                else
                                    disp('Shuffling option not implemented')
                                    return
                                end
                            end % end number of shufflings 
                        end % end computing null hypothesis
                        
                    end % end delay loop
                end 
            end % end time loop
        end % end feature loop
    end % end stimulus loop

    % Compute deltaDI
    for featIdx = 1:numel(params.selected_features)
        featLab = params.selected_features{featIdx};
        infoTransf_maps.(featLab).deltaDI = infoTransf_maps.(featLab).DI{2}-infoTransf_maps.(featLab).DI{1};
        infoTransf_mapsSh.(featLab).deltaDI = infoTransf_mapsSh.(featLab).DI{2}-infoTransf_mapsSh.(featLab).DI{1};
    end
end % end sess loop
end