function  Hstats = Hebb_stats(params,N,chev)

%HEBB_STATS -- Estimated PLM performance under Hebbian-learned weights.
%
%  Hstats = Hebb_stats(params)
%  Hstats = Hebb_stats(params,N)
%  Hstats = Hebb_stats(params,N,chev)
%
%  PARAMS is a parameter structure produced by PLM_PARAMS. 
%  The relevant fields are 'rep_gain', 'rep_noise', 'out_noise', and 
%  'W_init' (plus 'rep_size' and 'max_act' which have trivial roles only).
%  HEBB_STATS depends on the global variable PLM_CACHE (see CACHE_INPUTS).
%  N is a positive integer indicating the number of replications.
%  When not supplied, it defaults to PLM_CACHE(1).N.
%  Only the conditions with positive background orientations (arg(3)=2)
%  are used; the negative conditions are inferred by symmetry.
%  CHEV defaults to {randn(N,N_rep),randn(N,N_rep)} and contains the
%  Gaussian noise deviates for congruent and incongruent stimuli, resp.
%  It is relevant only when PARAMS.rep_noise > 0.
%
%  HSTATS is a structure with the following fields:
%   - N, rep_gain, rep_noise, out_noise, W_init -- copy of input params.
%   - con_idx -- PLM_CACHE indices used for the "congruent" conditions.
%   - incon_idx -- PLM_CACHE indices used for the "incongruent" conditns.
%   - mnXcon -- [3x35] mean activation over the representational layer
%       for congruent ("right") stimuli and contrast Cp=.23, .15, .10.
%       Both params.rep_noise and params.rep_gain are taken into account.
%   - mnXincon -- [3x35] mean representations, incongruent ("left") stim.
%   - sdXcon -- [3x35] the standard deviation of the congruent Xs.
%   - sdXincon -- [3x35] std.dev, incongruent Xs.  (Cf. AGGREG_CACHE)
%   - Wsame:  a structure array with fields
%     - weights -- [7x5] = mnXcon-mnXincon
%     - mnYcon -- [1x3] mean output, congruent stimuli, Cp=.23,.15,.10.
%     - mnYincon -- [1x3] mean output, incongruent stimuli, 3 contrasts
%     - sdY0con -- [1x3] the external-noise variability of congr. stimuli
%     - sdYincon -- [1x3] the external-noise var. of incongruent stimuli
%     - sdY -- total std.dev. of the output activation, incl. out_noise
%     - dprime -- [1x3] = (mnYcon - mnYincon) / sdY
%     - bias --   [1x3] = (mnYcon + mnYincon) / 2
%   - Wother: a structure array with the same fields as Wsame, showing
%       the performance in the "other" background (i.e., after a switch).
%   - Wmix:   ditto but based on averaged weights
%   - W0step: ditto but for step-like initial weights (see params.W_init)
%   - W0lin:  ditto but for linear initial weights (see params.W_init)
%
%  Example:
%   Idef=make_PLFrontend(PLM_input_params);PLM_CACHE=cache_inputs(Idef,N);
%   Mparams = PLM_params ; Hstats = Hebb_stats(Mparams);
%   print_Hebb_stats(Hstats) ; plot(Hstats.Wsame.weights,'.-') ;
%
%  See also PRINT_HEBB_STATS, PLM_HEBB1, CACHE_INPUTS, PLM_PARAMS.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.0 $  $Date: 2003/05/27 18:34 $

global PLM_CACHE OPTIMAL_W

if (nargin<2)  N = PLM_CACHE(1).N ; end    % use all reps by default

%-- Decipher and store the relevant model parameters
Hstats.N = N ;
rep_size = params.rep_size ; N_rep = prod(rep_size) ;
N_orient = rep_size(1) ; N_spfreq = rep_size(2) ;
max_act = params.max_act ;        % min_act assumed = -max_act
minus_rep_gain = -params.rep_gain ; Hstats.rep_gain = params.rep_gain ;
rep_noise = params.rep_noise ; Hstats.rep_noise = rep_noise ;
out_noise = params.out_noise ; Hstats.out_noise = out_noise ;
W_init = params.W_init ; Hstats.W_init = W_init ;

%-- Identify the conditions with "right" background orientation
arg = cat(1,PLM_CACHE(:).arg) ;   % [Gabor,contrast,background], [12x3]
con_idx = find((arg(:,1)==2)&(arg(:,3)==2)) ;    % R gabor, R backgr
incon_idx = find((arg(:,1)==1)&(arg(:,3)==2)) ;  % L gabor, R backgr
N_contrasts = length(con_idx) ;                  % should be 3
if (length(incon_idx)~=N_contrasts)
   error('The congruent and incongruent conditions mismatch.') ;
end
Hstats.con_idx = con_idx' ;                      % [8,10,12]
Hstats.incon_idx = incon_idx' ;                  % [7,09,11]

%-- Retrieve a sample of representation vectors from CACHE
Xcon   = zeros(N,N_rep,N_contrasts) ;   % [1000,35,3]
Xincon = zeros(N,N_rep,N_contrasts) ;
for k=1:N_contrasts
   Xcon(:,:,k)   = PLM_CACHE(con_idx(k)).rep(1:N,:) ;
   Xincon(:,:,k) = PLM_CACHE(incon_idx(k)).rep(1:N,:) ;
end

%-- Add representation noise, if specified.
%   Note that the same noise is used for all contrasts.
if (rep_noise>0)
   if (nargin<3) 
      Xcon = Xcon + repmat(randn(N,N_rep).*rep_noise,[1 1 N_contrasts]) ;
      Xincon = Xincon+repmat(randn(N,N_rep).*rep_noise,[1 1 N_contrasts]);
   else
      Xcon = Xcon + repmat(chev{1}.*rep_noise,[1 1 N_contrasts]) ;
      Xincon = Xincon + repmat(chev{2}.*rep_noise,[1 1 N_contrasts]) ;
   end
end

%-- Apply the sigmoidal activation function of the representation layer
Xcon   = max(0,max_act.*(2./(1+exp(minus_rep_gain.*Xcon))-1)) ;
Xincon = max(0,max_act.*(2./(1+exp(minus_rep_gain.*Xincon))-1)) ;

%-- Descriptive statistics of the representations
mnXcon = squeeze(mean(Xcon,1))' ;       % [3x35] <-- [:,35,3]
mnXincon = squeeze(mean(Xincon,1))' ;   % [3x35] <-- [:,35,3]
Hstats.mnXcon = mnXcon ;
Hstats.mnXincon = mnXincon ;

sdXcon = squeeze(mean(Xcon.^2,1))' ;
sdXcon = sqrt(sdXcon - mnXcon.^2) ;     % sqrt(mean(X^2)-mean(X)^2)
sdXincon = squeeze(mean(Xincon.^2,1))' ;
sdXincon = sqrt(sdXincon - mnXincon.^2) ;
Hstats.sdXcon = sdXcon ;
Hstats.sdXincon = sdXincon ;

%-- Estimate the Hebbian weights within the same (i.e. "right") context
%   W = Y.*X   averaged across contrasts, where Ycon=+1 and Yincon=-1
W = mean(mnXcon-mnXincon) ;   % [1x35]
W = reshape(W,rep_size) ;     % [7x5]
Hstats.Wsame = do_one_W(W,Xcon,Xincon,out_noise) ;
%Hstats.Wsame = do_one_W(OPTIMAL_W,Xcon,Xincon,out_noise) ;

%-- The weights for the "other" context are symmetrical w.r.t. vertical
Wo = -flipud(W) ;
Hstats.Wother = do_one_W(Wo,Xcon,Xincon,out_noise) ;

%-- Average weights (expected in a mixed environment)
Wm = (W + Wo) ./ 2 ;
Hstats.Wmix = do_one_W(Wm,Xcon,Xincon,out_noise) ;

%-- Two kinds of initial weights
if (N_orient~=7)
   warning('W0 is not implemented for N_orient~=7. Skipped.') ;
else
   W0s = W_init .* [-1 -1 -1 0 +1 +1 +1]' ;
   W0s = repmat(W0s,1,N_spfreq) ;
   Hstats.W0step = do_one_W(W0s,Xcon,Xincon,out_noise) ;
   W0l = W_init .* [-1.5 -1.0 -0.5  0  +0.5 +1.0 +1.5]' ;
   W0l = repmat(W0l,1,N_spfreq) ;
   Hstats.W0lin = do_one_W(W0l,Xcon,Xincon,out_noise) ;
end

%--- Return HSTATS
%%%% End of the main function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  Wstruct = do_one_W(W,Xcon,Xincon,out_noise)

w = W(:) ;    % [35x1]
N_contrasts = size(Xcon,3) ;
mnYcon = zeros(1,N_contrasts) ; sdY0con = zeros(1,N_contrasts) ;
mnYincon = zeros(1,N_contrasts) ; sdY0incon = zeros(1,N_contrasts) ;
for k=1:N_contrasts
   Y = Xcon(:,:,k) * w ;     % [Nx1] = [Nx35] * [35x1]
   mnYcon(k) = mean(Y) ;
   sdY0con(k) = std(Y,1) ;   % [1x1]
   Y = Xincon(:,:,k) * w ;
   mnYincon(k) = mean(Y) ;
   sdY0incon(k) = std(Y,1) ;
end
sdY = sqrt((sdY0con.^2+out_noise^2 + sdY0incon.^2+out_noise^2)./2) ;
dprime = (mnYcon - mnYincon) ./ sdY ;
bias   = (mnYcon + mnYincon) ./ 2 ;

Wstruct.weights = W ;
Wstruct.mnYcon = mnYcon ;
Wstruct.mnYincon = mnYincon ;
Wstruct.sdY0con = sdY0con ;
Wstruct.sdY0incon = sdY0incon ;
Wstruct.sdY = sdY ;
Wstruct.dprime = dprime ;
Wstruct.bias = bias ;

%-- Return WSTRUCT
%%% End of subfunction DO_ONE_W

%%%%% End of file HEBB_STATS.M
