function  [output,W_hist] = PLM_Hebb1(arg1,params,chev)

%PLM_HEBB1 -- Hebbian version of the Perceptual Learning Model.
%
%  [output,W_hist] = PLM_Hebb1(Mdata)
%  [output,W_hist] = PLM_Hebb1(stim,params)
%  [output,W_hist] = PLM_Hebb1(stim,params,chev)
%
%  STIM is is a matrix describing the trial sequence. It is usually
%  generated by BLOCK_FROM_CACHE and has the following five columns:
%   1. Gabor orientation: 1=-10 or 2=+10 deg, see PLEXP1_PARAMS
%   2. Gabor peak contrast: 1=0.10, 2=0.15 or 3=0.23, see PLEXP1_PARAMS
%   3. noise orientation: 1=-15 or 2=+15 deg, see PLEXP1_PARAMS
%   4. cache_idx: b/n 1 and 12  -- see global variable PLM_CACHE below
%   5. rep_idx: b/n 1 and PLM_CACHE(cache_idx).N
%  In addition, STIM may have a sixth column (=STIM(:,1) by default):
%   6. feedback: 0=none, 1="Left", 2="Right"
%
%  PARAMS is a structure produced by PLM_PARAMS. [Note: do not confuse
%  with the IDEF struct produced by MAKE_PLFRONTEND(PLM_INPUT_PARAMS).]
%  If PARAMS contains a field 'W0', it is used as the initial weight
%  vector, [N_rep+1+1 x 1]. Otherwise, PARAMS must contain two fields:
%  a scalar 'W_init' and a vector 'W0_seed'. W0 = W_init .* W0_seed.
%
%  The "chance event descriptor" CHEV can be one of the following:
%   1. A column vector of Gaussian "output noise" for each trial
%   2. A full NOISE matrix with N_trials rows and 1+N_rep+1 columns
%   3. A N_trials-by-2 "entry point" matrix generated by MAKE_CHEV.
%        It is passed to RANDOM_TABLE to generate a full NOISE matrix.
%   3. [] (the default) -- NOISE is generated by a call to RANDN.
%
%  Instead of supplying two (or three) input arguments, the same
%  information can be passed in a single structure MDATA with fields
%  'stim', 'params', and optionally 'chev'. If the first input
%  argument is a structure, no other arguments can be passed.
%
%  PLM_HEBB1 depends on the global variable PLM_CACHE, which must 
%  contain a struct generated by CACHE_INPUTS (see also AGGREG_CACHE).
%
%  OUTPUT is a matrix with N_trials rows and 5 columns:
%   1. response R(t) -- 1="Left" or 2="Right"
%   2. output activation O(t) 
%   3. bottom-up evidence I(t)
%   4. response criterion C(t)
%   5. feedback F(t) -- -1="Left" or +1="Right"
%  W_HIST(N+1,:)=running average of own responses,
%  W_HIST(N+2,:)=running average of teaching signal.
%
%  See also PLMODEL1, PLFRONTEND, PLZCONGR_HEBB1, PLM_PARAMS.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.3 $  $Date: 2005/05/03 $
%
% Part of the Hebbian Perceptual Learning Model version 1.1 for MATLAB.
% http://alexpetrov.com/proj/plearn/
% Copyright (c) Alexander Petrov 1999-2008, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in PercLearn_license.m

global PLM_CACHE

%-- Handle input arguments
if (isstruct(arg1))   % a single structure with subfields
   if (nargin~=1) warning('ARG1 is a struct; I ignore the rest.') ; end
   stim = arg1.stim ;
   params = arg1.params ;
   if (isfield(arg1,'chev'))  chev = arg1.chev ;  else chev = [] ; end
else  % STIM, PARAMS, and perhaps CHEV explicitly supplied
   stim = arg1 ;
   if (nargin<3) chev = [] ; end
end

%-- Decipher STIM
[N_trials,N_cols] = size(stim) ;
stim_type = stim(:,4) ;            % CACHE index, from 1 to 12
stim_idx = stim(:,5) ;             % index in CACHE(stim_type).rep
if (N_cols < 6)     % default FDBK = GABOR_ORIENT
   fdbk = stim(:,1) ;    % 1="Left", 2="Right"; recoded below
else    % explicit feedback (or lack thereof)
   fdbk = stim(:,6) ;    % 0=none, 1="Left", 2="Right"; recoded below
end

%-- Decipher the "chance event" descriptor CHEV
N_rep = prod(params.rep_size) ;    % number of representational units
rep_noise_std = params.rep_noise ; % s.d. of additive Gaussian noise
out_noise_std = params.out_noise ;
fdbk_fract = params.fdbk_fract ;
if (isempty(chev))      
   out_noise = randn(N_trials,1) .* out_noise_std ;
   if (rep_noise_std > 0)
      rep_noise = randn(N_trials,N_rep) .* rep_noise_std ;
   else
      rep_noise = zeros(N_trials,N_rep) ;
   end
   if (fdbk_fract>0 & fdbk_fract<1)
      z = randn(N_trials,N_rep) ;     % used with FDBK_FRACT later
   end
elseif (size(chev,2)==2)      % entry points
   if (rep_noise_std > 0)
      z = random_table(chev,1+N_rep+1) ;
      out_noise = z(:,1) .* out_noise_std ;
      rep_noise = z(:,2:N_rep+1) .* rep_noise_std ;
      z = z(:,end) ;          % used with FDBK_FRACT later
   else
      z = random_table(chev,1+1) ;
      out_noise = z(:,1) .* out_noise_std ;
      rep_noise = zeros(N_trials,N_rep) ;
      z = z(:,end) ;          % used with FDBK_FRACT later
   end
else    % explicit noise matrix supplied
   out_noise = chev(:,1) .* out_noise_std ;
   if (rep_noise_std > 0)
      rep_noise = chev(:,2:N_rep+1) .* rep_noise_std ;
   else
      rep_noise = zeros(N_trials,N_rep) ;
   end
   z = chev(:,end) ;          % used with FDBK_FRACT later
end

%-- Attend to the feedback on only a fraction of the trials
if (fdbk_fract==0)         % pure self-supervized mode
   fdbk = zeros(N_trials,1) ;
elseif (fdbk_fract < 1)    % feedback on some trials, own resp. on others
   % Clobber 1-FDBK_FRACT fraction entries in FDBK
   % (RAND would do just as well, but Z is used from CHEV...)
   fdbk_Z_cutoff = norminv(fdbk_fract) ;
   idx = find(z > fdbk_Z_cutoff) ;    % length(idx)/N_trials=1-fdbk_fract
   fdbk(idx) = 0 ;
%else (fdbk_fract==1)      % FDBK always attended
   % do nothing: FDBK is already set
end

%-- Recode FDBK from {1,2} format to {-1,+1} format
idx = find(fdbk~=0) ;           % recode 0-->0, 1-->-1, 2-->+1
fdbk(idx) =sign(fdbk(idx)-1.5); % 0=none, -1="Left",+1="Right"

%-- Decipher PARAMS
max_act = params.max_act ;        % min_act assumed = -max_act
minus_rep_gain = -params.rep_gain ; % gain of the repres-layer act.fun
Wmin = params.W_minmax(1) ; Wmax = params.W_minmax(2) ;
learn_rate = params.learn_rate * max_act ;
runav_rate = params.runav_rate ;
criterion = params.criterion ;
%trigger = params.trigger ;

%-- Initial weights, biases, and running averages
if (isfield(params,'W0'))    % initial weights explicitly supplied
   weights = params.W0 ;     % [N_rep,1]
else    % make W0 a multiple of the generic seed
   weights = params.W_init .* params.W0_seed(1:N_rep,1) ;   % [N_rep,1]
end

%-- Main loop
output = zeros(N_trials,5) ;           % [resp,act,bottom_up,crit,fdbk]
W_hist = zeros(N_rep+1+1,N_trials+1) ; % [weights;runav;teach_runav]
runav = 0 ;
teach_runav = 0 ;
W_hist(:,1) = [weights ; runav ; teach_runav] ;
for  tr = 1:N_trials
   %- Retrieve REP from PLM_CACHE, add noise, and pass through sigmoid
   rep = PLM_CACHE(stim_type(tr)).rep(stim_idx(tr),:) ;
   rep = rep + rep_noise(tr,:) ;   % scaled by REP_NOISE_STD, possibly 0
   rep = max(0,max_act.*(2./(1+exp(minus_rep_gain.*rep))-1)) ;
   %- Multiply by weights to get "bottom-up evidence"
   bottom_up = rep * weights ;      % [1x1] = [1xN]*[Nx1]
   %- Activation of the output unit; linear activation function
   act = bottom_up + out_noise(tr) ;
   %- Determine the current boundary between "Left" and "Right"
   %if (abs(runav)>trigger) crit=runav*criterion;else crit=0;end
   crit = runav * criterion ;
   %- Give out the overt response of the the model
   resp = 1 + (act>crit) ;      % 1="Left"=(act<=0), 2="Right"=(act>0)
   %- Hebbian weight update, hard feedback: "Left"=-1, "Right"=+1.
   %  The activation of the output unit is incorporated in LEARN_RATE.
   curr_fdbk = fdbk(tr) ;           % 0=none, -1="Left", +1="Right"
   if (curr_fdbk==0)    % no feedback available; use own RESP instead
      curr_fdbk = sign(resp-1.5) ;          % -1="Left", +1="Right"
   end
 %%dW = (learn_rate*curr_fdbk) .* rep' ;    % [Nx1], feedback version
 %%dW = (learn_rate*(curr_fdbk-runav)) .* rep'; % no-fdbk ver.1.2 Apr2004
   dW = (learn_rate*(curr_fdbk-teach_runav)) .* rep' ;  % ver.1.3 May2005
   %- Update the weights, soft weight bounding
   idx = find(dW<0) ; dW(idx) = dW(idx).*(weights(idx)-Wmin) ;
   idx = find(dW>0) ; dW(idx) = dW(idx).*(Wmax-weights(idx)) ;
   weights = weights + dW ;         % [Nx1]
   %- Running average of model's own responses
   runav = runav + runav_rate*(sign(resp-1.5)-runav) ;
   %- Running average of the teaching signal
   teach_runav = teach_runav + runav_rate*(curr_fdbk-teach_runav) ;
   %- Record output and weight history and move to the next trial
   output(tr,:) = [resp,act,bottom_up,crit,curr_fdbk] ;
   W_hist(:,tr+1) = [weights ; runav ; teach_runav] ;
end

%-- Return OUTPUT and W_HIST
%%%%% End of file PLM_HEBB1.M
