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

%PLM_HEBB1_ALT -- Alternative version of PLM_HEBB1.
%
%  [output,W_hist] = PLM_Hebb1_alt(Mdata)
%  [output,W_hist] = PLM_Hebb1_alt(repres,params,fdbk)
%  [output,W_hist] = PLM_Hebb1_alt(repres,params,fdbk,chev)
%
%  The canonical PLM_HEBB1 takes a STIM descriptor matrix as an input
%  argument and reads the representations from a global PLM_CACHE. This
%  alternative version takes the preprocessed representations instead.
%  REPRES is a NxK matrix of stimulus representations, where N is the
%  number of trials and K=prod(Mparams.rep_size)=35 is the number of
%  representational units under the current parametrization of the model.
%  Each row in REPRES must be generated by PLFRONTEND (or REPRES_PLE?_SEQ)
%  PLM_HEBB1_ALT does NOT depend on the global variable PLM_CACHE.
%  PARAMS is a structure produced by PLM_PARAMS. [Note: do not confuse
%  with the IDEF struct produced by MAKE_PLFRONTEND(PLM_INPUT_PARAMS).]
%  FDBK is a Nx1 vector of 1s, 2s, and possibly 0s. FDBK is a required
%  argument; use [] as shorthand for FDBK=zeros(size(REPRES,1),1).
%  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 three (or four) input arguments, the same
%  information can be passed in a single structure MDATA with fields
%  'rep', 'params', 'fdbk', and optionally 'chev'. If the first input
%  argument is a structure, no other arguments can be passed.
%
%  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 PLM_HEBB1, PLFRONTEND, REPRES_PLE4_SEQ, PLM_PARAMS.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.3 $  $Date: 2005/05/16 17:00 $

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

N_trials = size(repres,1) ;
if (isempty(fdbk))  fdbk = zeros(N_trials,1) ;  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 current representn, add noise, and pass through sigmoid
   %  NOTE: This line is the principal difference from PLM_HEBB1.M.
   rep = repres(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_ALT.M
