function [C, R, w, info] = rbf_bcm2(X, y, conf)
%
% Perform training of RBF network using BCM 
% regularization constraint by gradient descent. The  
% constraint is applied to the weights of RBF
%
% Input: 
%	X    - Training data set, where each column
%	       describes sample vector
%	y    - Training targets
%	conf - Configuration structure with the following
%	       relevant parameters:
%
%	         conf.cen     - RBF centers
%	         conf.rad     - RBF radii
%	         conf.u       - BCM regularization parameter
%	         conf.verb    - verbose output
%		 conf.scales  - List of RBF radii scales to try
%		 conf.lambdas - List of regularizaion scales to try

  % Standard parameters validity checks. Taken from Orr's code.

  % Program name (for error messages).
  prog = 'rbf_bcm2';

  % Configuration specification.
  spec(1) = struct( ...
    'comment', 'Centres', ...
    'name', 'cen', ...
    'type', 'matrix', ...
    'options', [], ...
    'default', []);
  spec(2) = struct( ...
    'comment', 'Radii', ...
    'name', 'rad', ...
    'type', 'matrix', ...
    'options', [], ...
    'default', []);
  spec(3) = struct( ...
    'comment', 'Model selection criterion', ...
    'name', 'msc', ...
    'type', 'string', ...
    'options', {{'uev', 'fpe', 'gcv', 'bic', 'mml'}}, ...
    'default', 'bic');
  spec(4) = struct( ...
    'comment', 'Termination threshold', ...
    'name', 'thresh', ...
    'type', {{'number', 'positive', 'integer'}}, ...
    'options', [], ...
    'default', 10000);
  spec(5) = struct( ...
    'comment', 'Hard limit on re-estimation iterations', ...
    'name', 'hard', ...
    'type', {{'number', 'positive', 'integer'}}, ...
    'options', [], ...
    'default', 100);
  spec(6) = struct( ...
    'comment', 'Basis function type', ...
    'name', 'type', ...
    'type', 'string', ...
    'options', [], ...
    'default', 'g');
  spec(7) = struct( ...
    'comment', 'Basis function sharpness', ...
    'name', 'exp', ...
    'type', 'number', ...
    'options', [], ...
    'default', 2);
  spec(8) = struct( ...
    'comment', 'Verbose output', ...
    'name', 'verb', ...
    'type', 'number', ...
    'options', {{0 1}}, ...
    'default', 0);
  spec(9) = struct( ...
    'comment', 'Display graphical timer', ...
    'name', 'timer', ...
    'type', 'string', ...
    'options', [], ...
    'default', '');
  spec(10) = struct( ...
    'comment', 'BCM learning rate', ...
    'name', 'u', ...
    'type', 'number', ...
    'options', [], ...
    'default', 0.1);
  spec(11) = struct( ...
    'comment', 'Basis function scales', ...
    'name', 'scales', ...
    'type', {{'vector', 'positive'}}, ...
    'options', [], ...
    'default', 1);
  spec(12) = struct( ...
    'comment', 'Initial regularisation parameter(s)', ...
    'name', 'lambdas', ...
    'type', {{'vector', 'positive'}}, ...
    'options', [], ...
    'default', 1);
  spec(13) = struct( ...
    'comment', 'BCM inhibition parameter', ...
    'name', 'l', ...
    'type', 'number', ...
    'options', [], ...
    'default', 0.1);

  % Check number of arguments.
  switch nargin
  case 3
  otherwise
    error([prog ': illegal number of arguments'])
  end

  % Check type of input arguments.
  if ~isnumeric(X) | ndims(X) ~= 2
    error([prog ': first argument (X) should be a matrix'])
  end
  if ~isnumeric(y) | ndims(y) ~= 2 | size(y,2) ~= 1
    error([prog ': second argument (y) should be a column vector'])
  end

  % Check for consistent size between X and y.
  [d,p] = size(X);
  if size(y,1) ~= p
    error([prog ': number of samples inconsistent between X and y'])
  end

  conf=conf_check(conf, spec, prog);

  % Check or set defaults for centres.
  if isempty(conf.cen)
    % The inputs are the centres, if the user doesn't say otherwise.
    conf.cen = X;
  else
    % The centre dimension must be equal to the input dimension.
    if size(conf.cen,1) ~= d
      error([prog ': inconsistent dimensions between X and conf.cen'])
    end
  end

  C = conf.cen;
  m = size(C,2);

  % Check or set defaults for radii.
  if isempty(conf.rad)
    % Default: scaled in each dimension, same for each centre.
    conf.rad = (max(X,[],2) - min(X,[],2));
    conf.rad = conf.rad(:,ones(m,1));
  else
    % Expand conf.rad using the same rules as for rbf_dm.
    if max(size(conf.rad)) == 1
      % It's a scalar - same size for each centre and in each dimension.
      conf.rad = conf.rad(ones(d,1),ones(m,1));
    elseif size(conf.rad,1) == 1 & size(conf.rad,2) == m
      % A row vector - same size in each dimension.
      conf.rad = conf.rad(ones(d,1),:);
    elseif size(conf.rad,1) == d & size(conf.rad,2) == 1
      % A column vector - same size for each centre.
      conf.rad = conf.rad(:,ones(m,1));
    elseif size(conf.rad,1) ~= d
      error([prog ': conf.rad has wrong input dimension'])
    elseif size(conf.rad,2) ~= m
      error([prog ': conf.cen and conf.rad are inconsistent'])
    end
  end

  %%%%%%%%%%%%%%% End of parameters checking %%%%%%%%%%%%%%%%%

  R = conf.rad;
  info.dmc = struct('type', conf.type, 'exp', conf.exp, 'bias', 0);
  info.verb = conf.verb;
  info.thresh = conf.thresh;
  info.l = conf.l;

  info.err = Inf;

  % Perform optimization of radiuses scale and the regularization
  % parameter

  for scale = conf.scales
    for lambda = conf.lambdas
    
      info.u = conf.u * lambda;

      % Train the network and get optimal weights for current choise
      % of scales and regularization parameter
      [w, err] = bcm_train(X, y, C, R .* scale, info);
      if err < info.err
	best_w = w;
        info.err = err;
	info.lam = lambda;
	info.scale = scale;        
      end
    end
  end
    
  w = best_w;
  R = info.scale .* R;


function [w, err] = bcm_train(X, y, C, R, info)
%
% Train the network with given centers and radii using BCM regularization
% constraint on the weights to prevent overfitting. info parameter
% includes the regularization parameter u. This function uses gradient
% descent method to reach minimal total error in weights space. The error
% is a sum of SSE with BCM projection index
%

  m = size(C,2);
  [d,p] = size(X);

  % Get starting weights
  H = rbf_dm(X, C, R, info.dmc);
  w = (pinv(H' * H) * H') * y + normrnd(0, .01, m, 1);

  % w = ones(m, 1);
  w1 = w;
  C1 = C;
  R1 = R;

  prev_err = inf;

  calcinfo = info.dmc;
  calcinfo.u = info.u;
  calcinfo.l = info.l;

  out = H*w;
  sse = sum((out-y).^2);
  risk1 = bcm_risk(out, info.u);
  [dc, dr, risk2] = CalcBCM(X, C, R, calcinfo);
  risk = risk1 + risk2;
  err = sse + risk;

  if info.verb
    disp(['sse = ',num2str(sse),', risk1 = ',num2str(risk1),', risk2 = ',num2str(risk2),', err = ',num2str(err)]);
  end

  % Gradient descent rate
  rate = 0.1;

  for i=1:info.thresh

    if rate <= eps
      break;
    end

    if prev_err >= err
      w = w1;
      C = C1;
      R = R1;

      if (prev_err - err) < eps
        break;
      end
      prev_err = err;

      dw1 = -2 .* ((out - y)' * H);
      dw2 = bcm_grad(H, out, info.u);
      dw2 = dw2 * norm(dw1);
      grad = (dw1 + dw2)';

      [dc1, dr1] = rbf_grad(X, H, out - y, C, R, w);     
%      [dc2, dr2, risk2] = CalcBCM(X, C, R, calcinfo);
%      dc = dc1 + dc2 * norm(dc1) / norm(dc2);
%      dr = dr1 + dr2 * norm(dr1) / norm(dr2);
      dc = dc1;
      dr = dr1;
   
      rate = 1.1 * rate;
    else
      rate = 0.5 * rate;
    end
    w1 = (1 - 0.1 * rate) * w + rate .* grad;
    C1 = (1 - 0.1 * rate) .* C + rate .* dc;
    R1 = (1 - 0.1 * rate) .* R + rate .* dr;

    H = rbf_dm(X, C, R, info.dmc);
    out = H*w1;
    sse = sum((out-y).^2);
    risk1 = bcm_risk(out, info.u);
    risk = risk1 + risk2;
    err = sse + risk;

    if info.verb
      disp(['sse = ',num2str(sse),', risk1 = ',num2str(risk1),', risk2 = ',num2str(risk2),', err = ',num2str(err)]);
    end
  end

function dw = bcm_grad(H, out, u)
%
% Gradient of risk function update according to BCM projection index
%

  p = size(out,1);

  % the design matrix is actualy the cell activation of all cells over all input
  c = out';

  theta = mean(c' .^ 2); 
  phi = (c - (ones(1,p)' * theta)') .* c;

  dw = phi * H;
  dw_norm = norm(dw);
  dw = u .* dw ./ p;

function risk = bcm_risk(out, u)
%
% BCM Risk function
%

  risk = -u * (mean(out.^3, 1)/3 - mean(out.^2, 1).^2/4);
  risk = sum(risk);

function [dm, dr] = rbf_grad(X, H, outdiff, C, R, w)

  [d,p] = size(X);
  m = size(C,2);

  dm = zeros(d,m);
  dr = zeros(d,m);

  for n=1:p
    for j=1:m
      diffvec = X(:,n) - C(:,j);
      hj = H(n,j);
      rj = R(:,j);
      dfedm = diffvec * hj ./ (rj.*rj); 
      rifac = (diffvec .^ 2) ./ (rj.*rj.*rj);
      dfedr = hj*rifac;
      del_mul_w = outdiff(n)' * w(j);
      dedmj = del_mul_w * dfedm;
      dedrj = del_mul_w * dfedr;
      dm(:,j) = dm(:,j) - dedmj;
      dr(:,j) = dr(:,j) - dedrj;
    end
  end

function [dm, dr] = rbf_grad2(X, H, outdiff, C, R, w)

  [d,p] = size(X);
  m = size(C,2);

  dm = zeros(d,m);
  dr = zeros(d,m);

  for j=1:m
    wj = w(j);
    mj = C(:,j);
    rj = R(:,j);
    rj2 = rj .* rj;
    rj3 = rj2 .* rj;
    hj = H(:,j)';
    hj = hj(ones(1,d),:);
  
    diff = X - mj(:,ones(1,p));
    s = sum(outdiff(:,ones(1,d))' .* hj .* diff, 2);
    dm(:,j) = 2 .* wj .* s ./ rj2;
    s = sum(outdiff(:,ones(1,d))' .* hj .* diff.^2, 2);
    dr(:,j) = 2 .* wj .* s ./ rj3;
  end
