function [arch, preproc]=findarch(X, y, conf, preproc)
% Learning and neural computation project by - Igor Bogodlov and Ron Sidi
% Function: findarch(X, y, conf, preproc)
%
% Finds the best architecture and creates a comitee of experts using the
% architecture, if preproc is not defined the best preprocessing will be
% searched for (otherwised the defined preproc will be done).
% Each expert is an ensemble of several networks, created by bootstrapping
% technique.
%
% Output:
%    arch - A structure representing the experts consisting of the following parameters:
%       nexperts - number of experts in the structure (the size of nclusters)
%       nclusters- a vector of the number of clusters (RBFs/hidden units) in each experts ensemble
%       experts  - an array of network ensembles (each entry is returned by train_expert)
%
% Input:
%    X - training set
%    y - output
%    conf - a configuration structure for rbf_dm with the following parameters added:
%       units - the number of neuron in each expert (the size of the vector determines the 
%               number of experts in the system)

disp('Learning and neural computation project by - Igor Bogodlov and Ron Sidi');
disp('Function: findarch(X, y, conf, preproc)');

[p,n] = size(X);
if size(y, 1) ~= p
   error(['Inconsistent dimensions between X and y']);
end

if nargin < 3
   conf = [];
end

if ~isfield(conf, 'nclusters')
   % arch.nclusters = [3, 5, 7, 9, 11]; 
   arch.nclusters = [4,5,6,7,8];
else
   arch.nclusters = conf.nclusters;
   rmfield(conf, 'nclusters');
end

% if preprocessing is not given then find the best preproc
if nargin < 4
   [processed, preproc] = findpreproc(X, y, arch);
else
   processed = findpreproc(X, y, arch, preproc);
end

disp(['Creating ', num2str(size(arch.nclusters, 2)), ' experts with the following number of hidden units: ']);
disp(['   ', num2str(arch.nclusters)]);

% the number of experts in the system is the number of different clusters in the preproc
arch.nexperts = size(preproc.clusterings, 2);

arch.experts = [];

% train the experts using the different number of clusters
for i = 1:arch.nexperts
   conf.cen = preproc.clusterings(i).cen';
   conf.rad = preproc.clusterings(i).rad';

   disp(['Training ensemble #',num2str(i),' out of ',num2str(arch.nexperts),'...']);
   disp(' ');

   E = train_expert(processed, y, conf);
   
   if isempty(arch.experts)
      arch.experts = [E];
   else
      arch.experts = [arch.experts, E];
   end
end

function E=train_expert(X, y, conf)
% Learning and neural computation project by - Igor Bogodlov and Ron Sidi
% Function: train_expert(X, y, conf)
%
% Train expert - ensemble of networks with same
% architecture. The function performs preprocessing
% on the training set, or, if argument conf.clusters is a matrix,
% uses it as clusters matrix.
%
%
% Output:
%
%    E - structure representing the ensemble with the following parameters:
%       W - weight vector for calculating prediction of the entire ensemble
%       s - number of network in ensemble
%       ens - a vector of networks, each networ is a structure with the following parameters:
%         C1, R1, w1 - centers and radiuses of network that predicts the correct output
%         C2, R2, w2 - centers and radiuses of network that predicts the inverted output
%         dmc1, dmc2 - the design matrix configuration structure for both of the netowrks
%         W          - weight matrix for reconstruction (used for confidence)
%         err        - the error of the network on the test set   
%
% Input:
%    X       - training set 
%    y       - training targets
%    conf    - options structure with following parameters:
%      func       - training function to use (defaults to rbf_bcm)
%      clusters   - number of clusters for clustering algorithm to find
%      cerr       - allowed error for clustering algorithm
%      mincluster - minimal number of points in cluster 
%      ensize     - ensemble size
%      cen	       - cluster centers
%      rad	       - cluster radii
%      u	       - regularization parameter
%      scales	  - a list of radii scales to try
%

[p,n] = size(X);
if size(y, 1) ~= p
   error(['Inconsistent dimensions between X and y']);
end

if ~isfield(conf, 'func')
   conf.func = 'rbf_bcm';
end
if ~isfield(conf, 'cerr')
   conf.cerr = eps;
end
if ~isfield(conf, 'ensize')
   conf.ensize = floor(2 * p / log(p));
end

% if clustering is requested perform it according to the given parameters
if isfield(conf, 'clusters')
   
   k = conf.clusters;
   
   disp('Performing K-Means clustering...');
   if isfield(conf, 'mincluster')
      [cen, rad] = mykmeans(k, X, conf.cerr, conf.mincluster);
   else
      [cen, rad] = mykmeans(k, X, conf.cerr);
   end
   
   conf.cen = cen';
   conf.rad = rad';
else
   % otherwise just add the parameter so that rmfield could remove it
   conf.clusters = p;
end

% add the parameter so that rmfield could remove it
if ~isfield(conf, 'mincluster')
   conf.mincluster = 1;
end

% Set up semi-logarithmic scale
% range = 1:9;
% conf.scales = [.1 .* range, range, 10 .* range];

conf.scales = 1:5;

% Number of ensembles
s = conf.ensize;

% Set up configuration structure for training routine
conf2 = rmfield(conf, {'func','clusters','mincluster','cerr','ensize'});

% Train the ensemble
E = train_ensemble(conf.func, X, y, conf2, s);

% Now we calculate weights for each network in ensemble
% based on their predictions

E.W = rank_ensemble(E, X, y);

function E=train_ensemble(func, X, y, conf, s)
%
% Train ensemble of networks. Each network is trained on the
% replicated fraction of the training set using bootstrapping
% technique with noise injection. Signals, which don't participate
% in training are used for cross-validation accross the ensemble. Each
% signal is thus used to validate different subset of ensemble members,
% and the scores are combined to give overall mean square error.
% As soon as the cross-validation error stops decreasing, the construction
% of ensemble stops. Thus best number of networks in the ensemble is
% determined on-line, from training results.
% All networks in the ensemble use the same architecture, defined in
% conf structure.
%

% number of ensembles
E.s = s;
ens = [];
[m,n] = size(X);

% Lists of network where given signal doesn't participate in
% the training - those sets are used for cross-validation over
% the ensemble
cvsets = cell(m, 1);
cvresults = cell(m, 1);

prev_cverr = Inf;

for j=1:s
   
   disp(['Training network #',num2str(j),' out of ',num2str(s),' ...']);
   
   % Create a dataset for the network by replication - bootstraping
   idxs = unidrnd(m,1,m);
   XX = X(idxs, :);
   yy = y(idxs);

   % and add noise to make the dataset independent
   XX = XX + normrnd(0,.001,m,n);
   
   % Train two networks one with the given output value and one with the inverted output value (this can be
   % done since the output is binary, it makes prediction better)
   [F.C1, F.R1, F.w1, info1] = feval(func, XX', yy, conf);
   [F.C2, F.R2, F.w2, info2] = feval(func, XX', ~yy, conf);
   
   % save the network design matrix configuration structure
   F.dmc1 = info1.dmc;
   F.dmc2 = info2.dmc;
   
   % Get weight matrix for reconstruction
   H = rbf_dm(XX', F.C1, F.R1, F.dmc1);
   F.W = pinv(H)*XX;
   
   % Use test set to estimate the error of each netwrok, this will be used to calculate the weighting of the networks later
   
   % create the test set, all of the samples that were not in the training set are the test set
   idxs_test = setdiff(1:m, unique(idxs));
   
   % find the error between the prediction of the trained network and the actual output
   Ht1 = rbf_dm(XX(idxs_test,:)', F.C1, F.R1);
   Ht2 = rbf_dm(XX(idxs_test,:)', F.C2, F.R2);
   testerr = (Ht1*F.w1 + 1 - Ht2*F.w2)./2 - y(idxs_test);
   F.err = mean(testerr, 1);
   
   % Try to estimate current cross-validation error among the ensemble.
   % In order to obtain realistic estimate of the ensemble, we calculate
   % cross validation of the signals only in networks where they don't
   % participate in training
 
   % Add new cv-set and cross-validation results
   for i=1:size(idxs_test,2)
     t = idxs_test(i);
     cvsets{t} = [cvsets{t}, j];
     cvresults{t} = [cvresults{t}, testerr(i,:)];
   end

   % Estimate current ensemble by giving cross-validation score for
   % sets among the ensemble, which don't participate in training.
   cverr = 0; evaluated = 0;
   for i=1:m
     if ~isempty(cvresults{i})
       evaluated = evaluated + 1;
       cverr = cverr + mean(cvresults{i},2).^2;
     end
   end
   cverr = cverr / evaluated;

   disp(['cverr = ',num2str(cverr), ', evaluated = ', num2str(evaluated)]);

   % If all points has been included into the estimate, and the cross
   % validation error has stopped decreasing, stop ensemble construction
   if (evaluated == m) & (cverr > prev_cverr)
     j = j - 1;
     break;
   end
   prev_cverr = cverr;

   % Add network to the ensemble
   if isempty(ens)
      ens = [F];
   else
      ens = [ens, F];
   end
end

E.ens = ens;
E.s = j;
E.cvresults = cvresults;
E.cverr = prev_cverr;
% E.stat_errs = errs;
% E.stat_evaln = evaln;


function W=rank_ensemble(E, X, y)
%
% Get weights corresponding to 'goodness' of each network in
% the ensemble
%

% create a vector of errors of each network in the enseble
Y = [];
for i=1:E.s
   Y = [Y, E.ens(i).err];
end

% calculate the covariance matrix of errors in the ensemble
C = Y' * Y;
% Normally, Cij would have to be E[ei*ej], covariance
% matrix of errors in the ensemble. But since each error
% is estimated by different test set, and each test set is 
% added with a random noise, this would no longer be correct.
% Instead, we make somewhat inaccurate assumption that the
% test sets are independent, and then Cij is simply E[ei]E[ej].

invC = pinv(C);

% Normalized weights are Bayesian estimation of posterior predictions 
W = sum(invC) / sum(sum(invC));


function [C, R, w, info] = rbf_bcm(X, y, conf)
%
% Perform training of RBF network using BCM 
% regularization constraint by gradient descent. The  
% constraint is applied to the weights of RBF
%
% Output:
%    C - the centers of the RBFs
%    R - the radiuses of the RBFs
%    w - the weights of the RBFs used to calculate the output
%    info - a structure consisting of the following relevant parameters:
%        lam   - the lambda used in the best network
%        scale - the scale used in the best network
%        err   - the estimated error of the best network
%        dmc   - a structure that describes the used design matrix (used by rbf_dm)
%
% Input: 
%    X    - Training data set, where each column
%           describes sample vector
%    y    - Training targets
%    conf - Configuration structure for rbf_dm with the following
%           parameters added:
%
%           u      - BCM regularization parameter
%           thresh - a limit on the number of iterations done by the gradient descent algorithm

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

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

% 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', '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);

% 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.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 using the analitical solution with no weights regularization, and add a little noise
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;

prev_err = inf;

out = H*w;
sse = sum((out-y).^2);
risk = bcm_risk(out, info.u);
err = sse + risk;

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

% Gradient descent rate (the rate is dynamic and will change throughout the execution of the alg)
rate = 0.1;

% don't do too many iterations
for i=1:info.thresh
   
   if rate <= eps
      break;
   end
   
   % if the last fix is good (reduced the error)
   if prev_err >= err
      % update the weights
      w = w1;
      
      % if the fix was not very relevant stop)
      if (prev_err - err) < eps
         break;
      end
      prev_err = err;
      
      % calculate the derivative of the sse according to the weights
      dw1 = -2 .* ((out - y)' * H);
      
      % calculate the derivative of the bcm risk according to the weights
      dw2 = bcm_grad(H, out, info.u);
      % make sure that the bcm derivatives and the sse derivatives have the same scale
      dw2 = dw2 * norm(dw1);
      
      % combime the two gradients
      grad = (dw1 + dw2)';
      
      % increase the learning rate in order to converge faster
      rate = 1.1 * rate;
   else
      % if the last fix wasn't good (the error was bigger then the previous iteration)
      % then reduce the learning rate and try again (without updating the weights)
      rate = 0.5 * rate;
   end
   
   % update the weights using the new gradient
   w1 = (1 - 0.1 * rate) * w + rate .* grad;
   
   % calculate new sse and bcmrisk
   out = H*w1;
   sse = sum((out-y).^2);
   risk = bcm_risk(out, info.u);
   err = sse + risk;
   
   if info.verb
      disp(['sse = ',num2str(sse),', risk = ',num2str(risk),', 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);
