function [R] = spm_dem_ERP(varargin)
% simulated electrophysiological response based on conditional estimates
% FORMAT [R] = spm_dem_ERP(qU,qU,...)
% qU - conditional estimates of states
% R  - summed response over peri-stimulus time
%
% These simulated response assume that LFPs are generated by superficial
% pyramidal cells that correspond to units encoding prediction error.
% Peristimulus time histograms (PSTH) assume that states (U) are encoded
% using a non-negative firing rate that is proportional to exp(U), using
% an opponent system
%__________________________________________________________________________


% loop over ERPs
%--------------------------------------------------------------------------
clf
N     = length(varargin);
color = {'r','b','g','y','c','m'};
for i = 1:N
    
    
    % loop over ERPs
    %----------------------------------------------------------------------
    qU = varargin{i};
    n  = length(qU.z);
    for  j = 1:n
        
        % PST (assuming 8ms times bins)
        %------------------------------------------------------------------
        try
            EEG = qU.Z{j};
        catch
            EEG = qU.z{j};
        end
        T   = (1:length(EEG))*32;
        
        % ERPs
        %------------------------------------------------------------------
        subplot(n,2,2*(j - 1) + 1)
        plot(T,mean(EEG,1),'Color',color{i}),hold on
        plot(T,EEG,':','Color',color{i})
        title(sprintf('LFPs (Causal): level %i',j),'FontSize',16)
        xlabel('pst (ms)')
        axis square
        
        
        % PSTH
        %------------------------------------------------------------------
        PSTH    = mean(exp(EEG) + exp(-EEG),1)/2;
        R{j}(i) = mean(PSTH);
        
        %         subplot(n,2,2*j)
        %         plot(T,PSTH,color{i},'LineWidth',1)
        %         title(sprintf('PSTH: level %i',j),'FontSize',16)
        %         xlabel('pst (ms)')
        %         axis square
        %         hold on
        
    end
    for  j = 1:(n - 1)
        
        % PST (assuming 8ms times bins)
        %------------------------------------------------------------------
        try
            EEG = qU.W{j};
        catch
            EEG = qU.w{j};
        end
        T   = (1:length(EEG))*32;
        
        % ERPs
        %------------------------------------------------------------------
        subplot(n,2,2*(j - 1) + 2)
        plot(T,mean(EEG,1),'Color',color{i}),hold on
        plot(T,EEG,':','Color',color{i})
        title(sprintf('LFPs (Hidden): level %i',j),'FontSize',16)
        xlabel('pst (ms)')
        axis square
        
    end
    
end