%% Clean up
clear;
close all;
clc;

%% Setup Parameters
ntrials = 1; % number of trials
dataset = 'MNIST'; % MNIST or YaleB (for YaleB, there's max 64 points for each person)
num_array = [0:1];  % which MNIST numbers to use? or, which YaleB people to use?
labels    = [1:2]; % labels for numbers
plist = [5000]; % number of data points in training set for EACH number (total of p*length(num_array) training points)
n = 800; % number of data points in testing set for EACH number (total of n*length(num_array) testing points)
mlist = [100 200]; % number of measurements
preprocess_data = 1;  % preprocess data? (0 = no, 1 = yes)
plot_data       = 1;
mu_train_test   = 0;  % include both training and testing data in mu/mean when preprocessing
nlayerslist     = 3;  % number of layers in neural net
RFlist         = [3];    % different options for redness factors
RFTYPES        = 6;    % number of RF types defined
saveFigs = 0;

% Parameter quality check
if ~( min(mlist) > max(nlayerslist) )
    error('Please ensure smallest m (min(mlist)) is larger than the largest number of layers (max(nlayerslist)))!');
end

G = length(num_array); % number of groups

%% Loop through training points, measurements, and trials

% Loop on number of layers
for iL = 1:length(nlayerslist)
    nlayers = nlayerslist(iL);
    % Loop on number of training points
    for ip = 1:length(plist)
        p = plist(ip);

        % Load training and testing data
        switch dataset

            case 'MNIST'
                [X_train b_train X_test b_test] = Load_MNIST(num_array, labels, p, n);

            case 'YaleB'
                [X_train b_train X_test b_test] = Load_YaleB(num_array, labels, p, n);
        end
        d = size(X_train,1);

        if (plot_data && ip == 1)
            % Plot Training and Testing Data
            faceW = sqrt(d); 
            faceH = sqrt(d); 
            numPerLine = p; 
            ShowLine = length(num_array); 

            Y = zeros(faceH*ShowLine,faceW*numPerLine); 
            for i=0:ShowLine-1 
                for j=0:numPerLine-1 
                    switch dataset
                        case 'MNIST'
                            Y(i*faceH+1:(i+1)*faceH,j*faceW+1:(j+1)*faceW) = reshape(X_train(:,i*numPerLine+j+1),[faceH,faceW])';
                            
                        case 'YaleB'
                            Y(i*faceH+1:(i+1)*faceH,j*faceW+1:(j+1)*faceW) = reshape(X_train(:,i*numPerLine+j+1),[faceH,faceW]); 
                    end
                end 
            end 

            figure; set(gcf,'WindowStyle','normal'); set(gcf,'PaperPositionMode','Auto');
            imagesc(Y);colormap(gray);
            %title(['Training p/' num2str(G) ' = ' num2str(p)],'fontsize',16);
            title('Training Data','fontsize',16);
            axis off
            set(gcf,'Position',[366 523 644 146]);
            filename = sprintf('Classify%sRF_TrainingData_%dlayers_%dp_%dn_%dd',dataset,nlayers,p,n,d);
            if (saveFigs)  
                print('-depsc','-r0',[filename '.eps'])
                saveas(gcf,[filename '.png'],'png');  
                saveas(gcf,[filename '.fig'],'fig');
            end

            numPerLine = n; 

            Y = zeros(faceH*ShowLine,faceW*numPerLine); 
            for i=0:ShowLine-1 
                for j=0:numPerLine-1 
                    switch dataset
                        case 'MNIST'
                            Y(i*faceH+1:(i+1)*faceH,j*faceW+1:(j+1)*faceW) = reshape(X_test(:,i*numPerLine+j+1),[faceH,faceW])'; 
                        case 'YaleB'
                            Y(i*faceH+1:(i+1)*faceH,j*faceW+1:(j+1)*faceW) = reshape(X_test(:,i*numPerLine+j+1),[faceH,faceW]); 
                    end
                end 
            end 

            figure; set(gcf,'WindowStyle','normal'); set(gcf,'PaperPositionMode','Auto');
            imagesc(Y);colormap(gray);
            %title(['Testing n/' num2str(G) ' = ' num2str(n)],'fontsize',16);
            title('Testing Data','fontsize',16);
            axis off
            set(gcf,'Position',[127 681 1086 116]);
            filename = sprintf('Classify%sRF_TestingData_%dlayers_%dp_%dn_%dd',dataset,nlayers, p,n,d);
            if (saveFigs)  
                print('-depsc','-r0',[filename '.eps'])
                saveas(gcf,[filename '.png'],'png');  
                saveas(gcf,[filename '.fig'],'fig');
            end
        end

        if (preprocess_data)
            % Rescale training data
            X_train_orig = X_train;
            if (mu_train_test)
                mu = sum([X_train X_test],2)/size([X_train X_test],2);
            else
                mu = sum(X_train,2)/size(X_train,2);
            end
            X_train = X_train - repmat(mu,1,p*length(num_array));

            % Rescale testing data
            X_test_orig = X_test;
            %mu = sum(X_test,2)/size(X_test,2);
            X_test = X_test - repmat(mu,1,n*length(num_array));
        end

        % Loop on number of measurements
        for im = 1:length(mlist)
            m = mlist(im);

            % Loop on trials
            for in = 1:ntrials
                tic;

                % Generate Layers and Redness Factors (RF) from Training Data      
                A = randn(m,d);
                Y_train = sign(A*X_train); 
                
                % Classify Testing Data
                Y_test = sign(A*X_test); 
                TestRF = zeros(RFTYPES,G,size(X_test,2)); % total redness factors for testing data

                level = [];
                for L = 1:nlayers
                   tic;
                   display(['L = ' num2str(L)]);
                   mm = 1;
                   otherIdx = zeros(1,L);
                   while (mm <= m)
                      if (L == 1)
                         level(L,mm).hyperplane_idx = mm; 
                      else
                         % choose L hyperplanes at random (no duplicates)
                         tmp = randperm(m);
                         idx_candidate = tmp(1:L);
                         idx_candidate = sort(idx_candidate);
                         if (sum(all(bsxfun(@eq,otherIdx, repmat(idx_candidate,size(otherIdx,1),1)),2)) > 0)
                             continue;
                         else
                            otherIdx(mm,1:L) = idx_candidate;
                            level(L,mm).hyperplane_idx = idx_candidate;
                         end
                      end

                      %level(L,mm).signPatterns = Y_train(level(L,mm).hyperplane_idx,:); % L x p
                      tmpSignPatterns = Y_train(level(L,mm).hyperplane_idx,:); % L x p
                      [tmpSignPatternsUnique, ia, ib] = unique(tmpSignPatterns','rows'); % return unique columns of level(L,mm).signPatterns (the unique sign patterns from the training data)
                      tmpSignPatternsUnique = tmpSignPatternsUnique';
                      
                      %level(L,mm).RF = nan(RFTYPES,G,size(tmpSignPatternsUnique,2));
                      tmpRF = nan(RFTYPES,G,size(tmpSignPatternsUnique,2));
                      
                      for pattern = 1:size(tmpSignPatternsUnique,2)
                          %locsOfThisPattern = all(bsxfun(@eq,tmpSignPatterns, ...
                          %                      repmat(tmpSignPatternsUnique(:,pattern),1,size(tmpSignPatterns,2))),1);
                          locsOfThisPattern = (pattern == ib);
                          
                          for iG = 1:G
                              % Group iG: count for pattern
                              %level(L,mm).groupCount(iG,pattern) = sum(locsOfThisPattern((iG-1)*p+1:iG*p));
                              tmpGroupCount(iG,pattern) = sum(locsOfThisPattern((iG-1)*p+1:iG*p));
                          end
                                            
                          % Group 1: count for pattern
                          %level(L,mm).groupCount(1,pattern) = sum(locsOfThisPattern(1:p));

                          % Group 2: count for pattern
                          %level(L,mm).groupCount(2,pattern) = sum(locsOfThisPattern(p+1:end));
                          
                          for irf = 1:length(RFlist)
                              
                              RF_type = RFlist(irf);
                          
                              switch RF_type

                                  case 1 % # of group g with pattern / total # with pattern
                                      level(L,mm).RF(RF_type,:,pattern) = level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern));           

                                  case 2 % # of group g with pattern
                                      level(L,mm).RF(RF_type,:,pattern) = level(L,mm).groupCount(:,pattern);

                                  case 3 % multiply by |RednessFactor - BluenessFactor|
                                      %level(L,mm).RF(RF_type,:,pattern) = (level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern))) * (abs(level(L,mm).groupCount(1,pattern)-level(L,mm).groupCount(2,pattern))./sum(level(L,mm).groupCount(:,pattern)));
                                      %level(L,mm).RF(RF_type,:,pattern) = (level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern))) .* (sum(abs(repmat(level(L,mm).groupCount(:,pattern),1,G)-repmat(level(L,mm).groupCount(:,pattern)',G,1)),2)./sum(level(L,mm).groupCount(:,pattern)));
                                      tmpRF(RF_type,:,pattern) = (tmpGroupCount(:,pattern)./sum(tmpGroupCount(:,pattern))) .* (sum(abs(repmat(tmpGroupCount(:,pattern),1,G)-repmat(tmpGroupCount(:,pattern)',G,1)),2)./sum(tmpGroupCount(:,pattern)));
                                      
                                  case 4 % multiply by (|RednessFactor - BluenessFactor| * # points in bin )/# training points total
                                      %level(L,mm).RF(RF_type,:,pattern) = (level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern))) * ((abs(level(L,mm).groupCount(1,pattern)-level(L,mm).groupCount(2,pattern)))/size(X_train,2));
                                      level(L,mm).RF(RF_type,:,pattern) = (level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern))) .* (sum(abs(repmat(level(L,mm).groupCount(:,pattern),1,G)-repmat(level(L,mm).groupCount(:,pattern)',G,1)),2)./size(X_train,2));

                                  case 5 % multiply by (|RednessFactor - BluenessFactor| * # points in bin * 2^L )/# training points total
                                      %level(L,mm).RF(RF_type,:,pattern) = (level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern))) * ((abs(level(L,mm).groupCount(1,pattern)-level(L,mm).groupCount(2,pattern)) * 2^L)/size(X_train,2));
                                      level(L,mm).RF(RF_type,:,pattern) = (level(L,mm).groupCount(:,pattern)./sum(level(L,mm).groupCount(:,pattern))) .* ((sum(abs(repmat(level(L,mm).groupCount(:,pattern),1,G)-repmat(level(L,mm).groupCount(:,pattern)',G,1)),2) * 2^L)./size(X_train,2));
                                      
                                  case 6 % # of group g with pattern / number of group g * |# of group 1 with pattern / number of group 1 - # of group 2 with pattern / number of group 2 | 
                                      %level(L,mm).RF(RF_type,:,pattern) = level(L,mm).groupCount(:,pattern)./p * abs(level(L,mm).groupCount(1,pattern)/p - level(L,mm).groupCount(2,pattern)/p);
                                      level(L,mm).RF(RF_type,:,pattern) = level(L,mm).groupCount(:,pattern)./p .* (sum(abs(repmat(level(L,mm).groupCount(:,pattern),1,G)/p-repmat(level(L,mm).groupCount(:,pattern)',G,1)/p),2));
                              end
                          
                          end
                      end
                      
                      % Apply to test data
                      for nn = 1:size(Y_test,2)
                            curTestSignPattern = Y_test(level(L,mm).hyperplane_idx,nn);                           
                            %curTestSignPatternLoc = all(bsxfun(@eq,tmpSignPatternsUnique, ...
                            %                            repmat(curTestSignPattern,1,size(tmpSignPatternsUnique,2))),1);
                            pattern = find(tmpSignPatternsUnique'*curTestSignPattern == L);
                            
                            %if (sum(curTestSignPatternLoc)>1)
                            if (length(pattern)>1)
                               error('Sign pattern should only match once!'); 
                            end
                            %pattern = find(curTestSignPatternLoc);
                            if (~isempty(pattern))
                                for irf = 1:length(RFlist)
                                    RF_type = RFlist(irf);
                                    TestRF(RF_type,:,nn) = TestRF(RF_type,:,nn) + tmpRF(RF_type,:,pattern);
                                end
                            end
                      end
                      
                      clear tmpSignPatterns;
                      clear tmpSignPatternsUnique;
                      clear tmpRF;
                      clear tmpGroupCount;
                      
                      mm = mm + 1;

                   end
                   toc
                end
                
                % Plot RF data?
%                 if (plot_data && in == 1 && im == length(mlist) && ip == length(plist))
%                     for irf = 1:length(RFlist)
%                         RF_type = RFlist(irf);
%                         figure; hold on; mymarkers = {'ro','bx','gs','md','k*','c<','y>','bo','rx','ms'};
%                         for L = 1:nlayers
%                            for mm = 1:m
%                                for iG = 1:G
%                                    plot((mm+(L-1)*m)*ones(size(level(L,mm).RF,3)), squeeze(level(L,mm).RF(RF_type,iG,:)),mymarkers{iG});
%                                end
%                            end
%                            plot(L*m+0.5*[1 1], get(gca,'ylim'),'k-');
%                         end
%                         xlabel('Node Index per Layer','fontsize',14);
%                         ylabel('Redness Factor (RF)','fontsize',14);
%                         title({['Redness Factors for Trial ' num2str(in) ' (group 1=red, group 2=blue)'];['nlayers = ' num2str(nlayers) ', m = ' num2str(m) ', p/' num2str(G) ' = ' num2str(p) ', RF-type = ' num2str(RF_type)]},'fontsize',14);
%                         set(gca,'fontsize',14);
%                         %ylim([0 1]);
%                         xlim([0 m*nlayers+1]);
%                         grid on;
%                         filename = sprintf('Classify%sRF_RF_%dlayers_%dm_%dp_%dn_%dd_%dpreprocess_%dmutesttrain_%dRFtype',dataset,nlayers,m,p,n,d,preprocess_data,mu_train_test,RF_type);
%                         if (saveFigs)  
%                             print('-depsc','-r0',[filename '.eps'])
%                             saveas(gcf,[filename '.png'],'png');  
%                             saveas(gcf,[filename '.fig'],'fig');
%                         end
%                     end
%                     
%                     figure; hold on;
%                     maxSignPatterns = 0;
%                     for L = 1:nlayers
%                        for mm = 1:m
%                            plot((mm+(L-1)*m), size(level(L,mm).signPatternsUnique,2),'b.');
%                            maxSignPatterns = max(maxSignPatterns,size(level(L,mm).signPatternsUnique,2));
%                        end
%                     end
%                     for L = 1:nlayers
%                         plot(L*m+0.5*[1 1], [0 plot(L*m+0.5*[1 1], [0 1],'k-');],'k-');
%                     end
%                     xlabel('Node Index per Layer','fontsize',16);
%                     ylabel('Number of Unique Sign Patterns','fontsize',16);
%                     title({['Number Unique Sign Patterns for Trial ' num2str(in) ', nlayers = ' num2str(nlayers)];['m = ' num2str(m) ', p/' num2str(G) ' = ' num2str(p) ', RF-type = ' num2str(RF_type)]},'fontsize',16);
%                     set(gca,'fontsize',14);
%                     xlim([0 m*nlayers+1]);
%                     ylim([0 maxSignPatterns+1]);
%                     grid on;
%                     filename = sprintf('Classify%sRF_numSignPatterns_%dlayers_%dm_%dp_%dn_%dd_%dpreprocess_%dmutesttrain_%dRFtype',dataset,nlayers,m,p,n,d,preprocess_data,mu_train_test,RF_type);
%                     if (saveFigs)  
%                         print('-depsc','-r0',[filename '.eps'])
%                         saveas(gcf,[filename '.png'],'png');  
%                         saveas(gcf,[filename '.fig'],'fig');
%                     end
%                     
%                 end

                TestRF = TestRF / (m*nlayers);
                for irf = 1:length(RFlist)
                    RF_type = RFlist(irf);
                    [junk b_test_guess] = max(squeeze(TestRF(RF_type,:,:)),[],1);
                    correct_label_cnt = sum(b_test_guess == b_test);

                    classificationRate(iL,ip,im,in,RF_type) = correct_label_cnt/size(X_test,2);

                    if (plot_data && in == 1 && im == length(mlist) && ip == length(plist))
                        figure; hold on; mymarkers = {'ro','bx','gs','md','k*','c<','y>','bo','rx','ms'};
                        for iG = 1:G
                            plot(1:size(X_test,2),squeeze(TestRF(RF_type,iG,:)),mymarkers{iG});
                        end
                        xlabel('Testing Data Point','fontsize',14);
                        ylabel('Total Redness Factor (RF)','fontsize',14);
                        title({['Test Data Redness Factors for Trial ' num2str(in) ' (group 1=red, group 2=blue)'];['nlayers = ' num2str(nlayers) ', m = ' num2str(m) ', p/' num2str(G) ' = ' num2str(p) ', RF-type = ' num2str(RF_type)]},'fontsize',14);
                        set(gca,'fontsize',14);
                        %ylim([0 1]);
                        xlim([0 size(X_test,2)]);
                        grid on;
                        filename = sprintf('Classify%sRF_TestRF_%dlayers_%dm_%dp_%dn_%dd_%dpreprocess_%dmutesttrain_%dRFtype',dataset,nlayers,m,p,n,d,preprocess_data,mu_train_test,RF_type);
                        if (saveFigs)  
                            print('-depsc','-r0',[filename '.eps'])
                            saveas(gcf,[filename '.png'],'png');  
                            saveas(gcf,[filename '.fig'],'fig');
                        end
                    end
                end
                
                elapsedTime = toc;
                disp(['Completed: nlayers=' num2str(nlayers) ', p=' num2str(p)  ', m=' num2str(m)  ', trial=' num2str(in) ', elapsed time=' num2str(elapsedTime)]);

            end
        end
    end
end

%% Plot Results

%%% Plot classificationRate %%%

mymarkers = {'o','+','s','d','<','>','^'};

% Classification rate versus number of measurements for a given number of training points p
for irf = 1:length(RFlist)
    RF_type = RFlist(irf);
    for iL = 1:length(nlayerslist)
        nlayers = nlayerslist(iL);
        figure; hold all;
        for ip = 1:length(plist)
            plot(mlist,mean(squeeze(classificationRate(iL,ip,:,:,RF_type)),2),'-','marker',mymarkers{ip},'linewidth',2);
            %legendString1{ip} = sprintf('p/%d = %d',G,plist(ip));
            legendString1{ip} = sprintf('%d Training per Group',plist(ip));
        end
        legend(legendString1,4);
        xlabel('Number of Measurements (m)','fontsize',16);
        ylabel('Average Correct Classification Rate (Total)','fontsize',16);
%         if (preprocess_data)
%             if (mu_train_test)
%                 mu_in_preprocess = 'mu from training and test data';
%             else
%                 mu_in_preprocess = 'mu from training only';
%             end
%             title({['Testing Size Per Group (n/' num2str(G) ') = ' num2str(n) ', Dim (d) = ' num2str(d) ', Trials = ' num2str(ntrials) ', nlayers = ' num2str(nlayers)];['Numbers: ' num2str(num_array) ', Preprocess = ' num2str(preprocess_data) ', ' mu_in_preprocess ', RF-type = ' num2str(RF_type)]},'fontsize',14);
%         else
%             title({['Testing Size Per Group (n/' num2str(G) ') = ' num2str(n) ', Dim (d) = ' num2str(d) ', Trials = ' num2str(ntrials) ', nlayers = ' num2str(nlayers)];['Numbers: ' num2str(num_array) ', Preprocess = ' num2str(preprocess_data) ', RF-type = ' num2str(RF_type)]},'fontsize',14);
%         end
        set(gca,'fontsize',16);
        ylim([0 1]);
        grid on;
        filename = sprintf('Classify%sRF_RateVersusM_%dlayers_%dtrials_%dn_%dgroup1_%dgroup2_%dpreprocess_%dmutesttrain_%dRFtype',dataset,nlayers,ntrials,n,num_array(1),num_array(2),preprocess_data,mu_train_test,RF_type);
        if (saveFigs)  
            print('-depsc','-r0',[filename '.eps'])
            saveas(gcf,[filename '.png'],'png');  
            saveas(gcf,[filename '.fig'],'fig');
        end
    end
end

% % Classification rate versus number of training points for a given number of measurements m
% for irf = 1:length(RFlist)
%     RF_type = RFlist(irf);
%     for iL = 1:length(nlayerslist)
%         nlayers = nlayerslist(iL);
%         figure; hold all;
%         for im = 1:length(mlist)
%             plot(plist,mean(squeeze(classificationRate(iL,:,im,:,RF_type)),2),'-','marker',myMarkers{ceil(ip/7)},'linewidth',2);
%             legendString2{im} = sprintf('m = %d',mlist(im));
%         end
%         legend(legendString2,4);
%         xlabel(['Number of Training Points Per Group (p/' num2str(G) ')'],'fontsize',16);
%         ylabel('Average Correct Classification Rate (Total)','fontsize',16);
%         if (preprocess_data)
%             if (mu_train_test)
%                 mu_in_preprocess = 'mu from training and test data';
%             else
%                 mu_in_preprocess = 'mu from training only';
%             end
%             title({['Testing Size Per Group (n/' num2str(G) ') = ' num2str(n) ', Dim (d) = ' num2str(d) ', Trials = ' num2str(ntrials) ', nlayers = ' num2str(nlayers)];['Numbers: ' num2str(num_array) ', Preprocess = ' num2str(preprocess_data) ', ' mu_in_preprocess ', RF-type = ' num2str(RF_type)]},'fontsize',14);
%         else
%             title({['Testing Size Per Group (n/' num2str(G) ') = ' num2str(n) ', Dim (d) = ' num2str(d) ', Trials = ' num2str(ntrials) ', nlayers = ' num2str(nlayers)];['Numbers: ' num2str(num_array) ', Preprocess = ' num2str(preprocess_data) ', RF-type = ' num2str(RF_type)]},'fontsize',14);
%         end
%         set(gca,'fontsize',14);
%         ylim([0 1]);
%         grid on;
%         filename = sprintf('Classify%sRF_RateVersusP_%dlayers_%dtrials_%dn_%dgroup1_%dgroup2_%dpreprocess_%dmutesttrain_%dRFtype',dataset,nlayers,ntrials,n,num_array(1),num_array(2),preprocess_data,mu_train_test,RF_type);
%         if (saveFigs)  
%             print('-depsc','-r0',[filename '.eps'])
%             saveas(gcf,[filename '.png'],'png');  
%             saveas(gcf,[filename '.fig'],'fig');
%         end
%     end
% end

% Example Image
for iG = 1:G
    figure;
    switch dataset
        case 'MNIST'
            if (preprocess_data)
                subplot(1,2,1);
                imagesc(reshape(X_train_orig(:,(iG-1)*p+1),28,28)'); colormap(gray); colormap;
                title('Example Original','fontsize',16);
                subplot(1,2,2);
                imagesc(reshape(X_train(:,(iG-1)*p+1),28,28)'); colormap(gray); colormap;
            else
                imagesc(reshape(X_train(:,(iG-1)*p-1),28,28)'); colormap(gray); colormap;
            end

        case 'YaleB'
            if (preprocess_data)
                subplot(1,2,1);
                imagesc(reshape(X_train_orig(:,(iG-1)*p+1),32,32)); colormap(gray); colormap;
                title('Example Original','fontsize',16);
                subplot(1,2,2);
                imagesc(reshape(X_train(:,(iG-1)*p+1),32,32)); colormap(gray); colormap;
            else
                imagesc(reshape(X_train(:,(iG-1)*p+1),32,32)); colormap(gray); colormap;
            end
    end
    if (preprocess_data)
        if (mu_train_test)
            mu_in_preprocess = 'mu from training and test data';
        else
            mu_in_preprocess = 'mu from training only';
        end
        title({['Example Training Point for Number ' num2str(num_array(iG))];['Preprocess = ' num2str(preprocess_data) ', ' mu_in_preprocess ', RF-type = ' num2str(RF_type)]},'fontsize',16);
    else
        title({['Example Training Point for Number ' num2str(num_array(iG))];['Preprocess = ' num2str(preprocess_data) ', RF-type = ' num2str(RF_type)]},'fontsize',16);
    end
    filename = sprintf('Classify%sRF_ExampleImage_Group%d_%dlayers_%dtrials_%dn_%dgroup1_%dgroup2_%dpreprocess_%dmutesttrain_%dRFtype',dataset,iG,nlayers,ntrials,n,num_array(1),num_array(2),preprocess_data,mu_train_test,RF_type);
    if (saveFigs)  
        print('-depsc','-r0',[filename '.eps'])
        saveas(gcf,[filename '.png'],'png');  
        saveas(gcf,[filename '.fig'],'fig');
    end
end


