%  [preproc, netarch] = findarch(train_data, train_label, options)

%

%    Find optimal architecture for given training data, using given options.

%

%    Input:   train_data  =  training set data

%             train_label =  training set target labels

%             options     =  parameters for architecture seek. The following struct:

%                            .preproc_options.method    = 'minimal entropy' OR

%                                                         'k-means clustering'

%                                            .inputdim  = desired input layer

%                                                         dimension (after reducing)

%                            .net_options.method   = 'kurtosis' OR 'bcm'

%

%    Output:  preproc     =  found preprocessing parameters

%             netarch     =  found network architecture parameters

%



function [preproc, netarch] = findarch(train_data, train_label, options)



global theta mom3 mom4



disp('findarch / Yuval Lifshitz, Liron Weissman, Gilad Arnold');



% check arguments

if (nargin < 3)

   error('findarch: invalid num of arguments');

end



% decode architecture and preprocessing options

procoptions = options.preproc_options;



if (strcmp(options.net_options.method, 'kurtosis'))

   method = 1;

elseif (strcmp(options.net_options.method, 'bcm'))

   method = 2;

else

   error('findpreproc: training constraint options not provided or invalid');

end



% definition of the arcitecture structure

netarch = struct('net', network, 'mse_test_error',0,'moment_test_error',0);



% Create and train the ensemble of networks:

% - We create an ensamble of 15 networks, each net is trained on 80% samples out of the 

%   possible samples. 

% - The choosing mechanism is random and uniform and enables duplicate rows.

dataSize = size(train_data, 1);

subset_size = ceil(dataSize*0.8);



% creating the error matrix for finding the best archtecture

% if the best error is bigger than 100 than probably our net is useless

avg_mse_test_error = zeros(ceil(subset_size / 2));

avg_mse_test_error = avg_mse_test_error + 100;



NNTWARN OFF;



% iterate over input dimension to find optimal dimension

for inputdim = ceil(subset_size / 2) : -1 : ceil(sqrt(subset_size)) 

   

   % use auxiliary function to find preprocessing down to set dimension

   preproc = [];

   procoptions.inputdim = inputdim;

   [processed, preproc] = findpreproc(train_data, train_label, procoptions, preproc);

      

   % iterate over hidden layer dimension to find optimal dimention

   for layerdim = inputdim : -1 : ceil(inputdim / 2)

      

      % train ensemble networks

      for i = 1 : 15

         % choose random vectors and their matching targets (uniform)

         randIndx = ceil(rand(subset_size, 1) * dataSize);

         data_subset   = processed(randIndx, :)';

         target_subset = train_label(randIndx)';

         

         k = 0;

         testIndx = [];

         % Find the complement of these vectors for testing  

         for j = 1 : dataSize

            if (sum(j == randIndx) == 0)

               k = k + 1;

               testIndx(k) = j;

            end

         end

         

         data_test   = processed(testIndx, :)';
         target_test = train_label(testIndx)';

         

         min_max_values = minmax(data_subset);

         netarch(i).net = newff(min_max_values,[layerdim,1], ...

            {'tansig', 'tansig'} ,'trainlm', 'learngd', 'mse');

         

         w1 = netarch(i).net.IW{1};

         b1 = netarch(i).net.b{1};

         w2 = netarch(i).net.LW{2};

         b2 = netarch(i).net.b{2};

         tp(1) = 25;

         tp(2) = 1;

         tp(3) = 500;

         tp(4) = 0.02;

         tp(5) = 0.1;

         tp(6) = 0.1;

         tp(7) = 0.02;

         tp(8) = 0.9;

         theta = [];

         mom3  = [];

         mom4  = [];

         if method == 1

            % traing back-prop with Kurtosis

            [w1,b1,w2,b2,epoch,tr,risk_rec] = ...

               tbpkurtosis2(w1,b1,'tansig',w2,b2,'tansig',data_subset, ...

               target_subset,tp,data_test,target_test);

         else

            % train back-prop with BCM

            [w1,b1,w2,b2,epoch,tr,risk_rec] = ...

               tbpbcm2(w1,b1,'tansig',w2,b2,'tansig',data_subset,target_subset, ...

               tp,data_test,target_test);

         end

         

         % store resulting architecture for later use

         netarch(i).net.IW{1} = w1;

         netarch(i).net.b{1}  = b1;

         netarch(i).net.LW{2} = w2;

         netarch(i).net.b{2}  = b2;

         netarch(i).mse_test_error    = tr;

         netarch(i).moment_test_error = risk_rec;   

      end

            

      % calculate the average test error of the ensamble

      tmp_mse_test_error = 0;

      for i = 1 : 15

         tmp_mse_test_error = tmp_mse_test_error + mean(netarch(i).mse_test_error);

      end

      avg_mse_test_error(inputdim,layerdim) = tmp_mse_test_error / 15;      

   end

      

end



% find the best input dimention and hidden layer dimention from the error matrix

[err_vec inpudim_vec_idx] = min(avg_mse_test_error);

[min_err layerdim] = min(err_vec);

inputdim = inpudim_vec_idx(layerdim);

   

% rebuild the ensamble for the chosen architecture

% train ensemble networks

for i = 1 : 15

   % choose random vectors and their matching targets (uniform)

   randIndx = ceil(rand(subset_size, 1) * dataSize);

   data_subset   = processed(randIndx, :)';

   target_subset = train_label(randIndx)';

   

   k = 0;

   testIndx = [];

   % Find the complement of these vectors for testing  

   for j = 1 : dataSize

      if (sum(j == randIndx) == 0)

         k = k + 1;

         testIndx(k) = j;

      end

   end

   

   data_test   = processed(testIndx, :)';
   target_test = train_label(testIndx)';

   

   min_max_values = minmax(data_subset);

   netarch(i).net = newff(min_max_values,[layerdim,1], ...

      {'tansig', 'tansig'} ,'trainlm', 'learngd', 'mse');

   

   w1 = netarch(i).net.IW{1};

   b1 = netarch(i).net.b{1};

   w2 = netarch(i).net.LW{2};

   b2 = netarch(i).net.b{2};

   tp(1) = 25;

   tp(2) = 1;

   tp(3) = 500;

   tp(4) = 0.02;

   tp(5) = 0.1;

   tp(6) = 0.1;

   tp(7) = 0.02;

   tp(8) = 0.9;

   theta = [];

   mom3  = [];

   mom4  = [];

   if method == 1

      % traing back-prop with Kurtosis

      [w1,b1,w2,b2,epoch,tr,risk_rec] = ...

         tbpkurtosis2(w1,b1,'tansig',w2,b2,'tansig',data_subset, ...

         target_subset,tp,data_test,target_test);

   else

      % train back-prop with BCM

      [w1,b1,w2,b2,epoch,tr,risk_rec] = ...

         tbpbcm2(w1,b1,'tansig',w2,b2,'tansig',data_subset,target_subset, ...

         tp,data_test,target_test);

   end

   

   % store resulting architecture for later use

   netarch(i).net.IW{1} = w1;

   netarch(i).net.b{1}  = b1;

   netarch(i).net.LW{2} = w2;

   netarch(i).net.b{2}  = b2;

   netarch(i).mse_test_error    = tr;

   netarch(i).moment_test_error = risk_rec;   

end



NNTWARN ON;

