function [dRdm, dRdr, risk] = CalcBCM(X, C, R, conf)
% Learning and neural computation project by - Igor Bogodlov and Ron Sidi
% Function: CalcBCM(X, C, R, conf)
%
% Calculate the BCMRisk derivative according to each of the RBFs
% center and radius, from this it is possible to calculate the 
% unsupervised learning rule for both the centers and the radiuses
%
% Output: dRdm - the derivative of the BCM risk according to the
%                RBFs centers (an array of m center vectors)
%         dRdr - the derivative of the BCM risk according to the
%                RBFs radiuses (an array of m radius vectors)
%
% Input: X - the training set
%        C - the centers of the RBFs [m,d]
%        R - the radiuses of the RBFs:
%            [m,d] - the radiuse of all RBFs in all dimensions
%            [1,d] - all RBFs have same radius, radius is given for each dimension
%            [m,1] - every RBF has different radius but the same for all dimensions
%            [1,1] - all RBFs has the same radius in all dimensions
%        conf - a structure for configuring rbf_dm with the additional fields:
%            l - Between cell threshold (for calculating the inhibited activation)
%            u - Single cell learning rate
%

% Program name (for error messages) and configuration spec.
prog = 'CalcBCMRiskDerivatives';
spec(1) = struct( ...
  'comment', 'Type of RBF function', ...
  'name',    'type', ...
  'type',    'string', ....
  'options', {{{'g', 'gaussian'}, {'c', 'cauchy'}, {'m', 'multiquadric'}, {'i', 'inverse', 'inverse multiquadric'}}}, ...
  'default', 'g');
spec(2) = struct( ...
  'comment', 'Sharpness of RBF function', ...
  'name',    'exp', ...
  'type',    {{'number', 'positive'}}, ...
  'options', [], ...
  'default', 2);
spec(3) = struct( ...
  'comment', 'Location of bias unit', ...
  'name',    'bias', ...
  'type',    {{'number', 'nonnegative', 'integer'}}, ...
  'options', [], ...
  'default', 0);
spec(4) = struct( ...
  'comment', 'Single cell learning rate', ...
  'name',    'u', ...
  'type',    {{'number', 'positive'}}, ...
  'options', [], ...
  'default', 0.1);
spec(5) = struct( ...
  'comment', 'Between cell threshold', ...
  'name',    'l', ...
  'type',    {{'number', 'positive'}}, ...
  'options', [], ...
  'default', 0.01);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% code taken from Orrs rbf_dm to check and set the arguments to the correct dimension

% Check number of parameters. Take special action if only one or two.
switch nargin
case 1
  if isstring(X)
    switch X
    case 'conf'
      conf_print(prog, spec)
      return
    otherwise
      error([prog ': ''' X ''' unrecognised for single string argument'])
    end
  else
    error([prog ': unrecognised type for single argument'])
  end
case 2
  if isstring(X) & isstring(C)
    switch X
    case 'conf'
      conf_print(prog, spec, C)
      return
    otherwise
      error([prog ': ''' X ''' unrecognised for two string arguments'])
    end
  else
    error([prog ': unrecognised type(s) for two arguments'])
  end
case 3
  conf = [];
case 4
otherwise
  error([prog ': wrong number of arguments'])
end

% Check type of input arguments.
if ~isnumeric(X) | ndims(X) ~= 2
  error([prog ': argument X should be a matrix'])
end
if ~isnumeric(C) | ndims(C) ~= 2
  error([prog ': argument C should be a matrix'])
end
if ~isnumeric(R) | ndims(R) ~= 2
  error([prog ': argument R should be a matrix, vector or scalar'])
end

% Check for consistent dimensions amongst X, C and R. Reshape R if necessary.
[d,p] = size(X);
if size(C,1) ~= d
  error([prog ': input dimension inconsistent between X and C'])
end
m = size(C,2);
if max(size(R)) == 1
  R = R(ones(d,1),ones(m,1));
elseif size(R,1) == 1 & size(R,2) == m
  R = R(ones(d,1),:);
elseif size(R,1) == d & size(R,2) == 1
  R = R(:,ones(m,1));
elseif size(R,1) ~= d
  error([prog ': argument R has wrong input dimension'])
elseif size(R,2) ~= m
  error([prog ': number of basis functions inconsistent between C and R'])
end

% Check the configuration is okay and set defaults (if required).
conf = conf_check(conf, spec, prog);

%%%%%%%%% end of code taken from Orrs rbf_dm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the design matrix is actualy the cell activation of all cells over all input
dmc = struct('type', conf.type, 'exp', conf.exp, 'bias', 0);
c = rbf_dm(X, C, R, dmc);

% calculate c' (using conf.l as a threshold)
% c'k = ck - l +(j!=k) cj
ctag = zeros(p, m);
for i = 1:m
   ctag(:,i) = c(:,i) - conf.l * (sum([c(:,1:(i-1)) c(:,(i+1):m)], 2));
end

% calculate the risk for every neoron
risk = mean(ctag.^3, 1)/3 - mean(ctag.^2, 1).^2/4;

% the global risk is sum of all the local risks
risk = sum(risk);

theta = mean(ctag.^2, 1);

dctagdm = zeros(d, p, m);
dctagdr = zeros(d, p, m);

% calculate dc'k/dmk and dc'/drk
% there is no need to calculate dc'k/dmj since it is equal to -l * dc'j/dmj
for i = 1:m
   curcenter = C(:, i);

   % dc'k/dmk = (x - mk)/rk^2 * ck
   rvec = 2.*(R(:,i).^(-2));
   dctagdm(:,:,i) = rvec(:,ones(1,p)) .* (X - curcenter(:, ones(1,p))) .* c(:,repmat(i,1,d))';
   
   %dc'k/drk = 2 * (x - mk)^2 / rk^3 * ck
   rvec = 2.*(R(:,i).^(-3));
   dctagdr(:,:,i) = rvec(:,ones(1,p)) .* ((X - curcenter(:, ones(1,p))).^2) .* c(:,repmat(i,1,d))';
end

% calculate dR/dmk
% calculate f(c'k, theta_k) (there is no need for f(c'k, theta_j) )
% f(c,t) = c^2 - t*c,
f = ctag.^2 - ctag.*theta(ones(1,p), :);
mean_f_dctagdm = zeros(d,m);
mean_f_dctagdr = zeros(d,m);

% calculate E[f(c'k, theta_k) * dc'k/dmk] and E[f(c'k, theta_k) * dc'k/drk]
for i = 1:m
   mean_f_dctagdm(:,i) = mean(f(:,repmat(i,1,d)) .* dctagdm(:,:,i)', 1)';
   mean_f_dctagdr(:,i) = mean(f(:,repmat(i,1,d)) .* dctagdr(:,:,i)', 1)';
end
   
dRdm = zeros(d, m);
dRdr = zeros(d, m);

% now calculate the risk derivatives
for i = 1:m
%dR/dmk = -u ( E[f(c'k, theta_k) dc'k/dmk] - l +(j!=k) E[f(c'j, theta_j) dc'j/dmj] )   
   dRdm(:,i) = -conf.u * (mean_f_dctagdm(:,i) - conf.l * sum(mean_f_dctagdm(:,[1:(i-1),(i+1):m]), 2));
%dR/drk = -u ( E[f(c'k, theta_k) dc'k/drk] - l +(j!=k) E[f(c'j, theta_j) dc'j/drj] )
   dRdr(:,i) = -conf.u * (mean_f_dctagdr(:,i) - conf.l * sum(mean_f_dctagdr(:,[1:(i-1),(i+1):m]), 2));
end