function info(layer)
    % Single and multiple cell information analysis for any layer 1-4
    % NBESTCELLS cells is the number of cells for each stimulus (default 5) 
    % copyright E T Rolls, Oxford Centre for Computational Neuroscience
    
    % Create infos_matlab.mat for infors_visnet and inform_visnet
    % infos_matlab.mat contains:
    % num_c the number of neurons
    % num_s the number of stimuli 
    % infom: each row is for one trial and has the stimulus number followed by the rates for each neuron
    % Each rate is for a different transform (view)
    % the stimulus numbers start at 1. Each trial (row) can be in any order.
    
global EFR EFRDATA NLAYERS NET_SIZE SQNET_SIZE N_FILES N_GROUPS N_VIEWS N_LOCS BESTCELLS NBESTCELLS ...
       NBESTCELLSINFORM NTRANSFORMS PARAMETER_LIST SYNMAT SYNMAT_CHANGE ...
       SYNMAT_CHANGE_GROUP LOGNUM FIXEDTRACE SPARSENESS BETA ETA LRNRATE NSYN SYNMATL1 SYNMATL2 SYNMATL3 ...
       SYNMATL4 VNLOGINDEX IMLIST IMLIST_TEST N_VIEWS_TEST FILTIM NSYNL1 Images base_dir image_dir ...
       FiltImSize SqFiltImSize Display ...
       filter_layer_1 filter_layer_2  filter_layer_3 filter_layer_4 SIGMOID LATERALINHIB TraceDisp DisplaySYNMATprevious
    
    if layer < 1
        layer = 1
    elseif layer > 4
        layer = 4
    end
    disp(['Single Cell Information Analysis using infors_visnet for layer (1-4) ', num2str(layer)]);
    
    infom = zeros(N_GROUPS*N_VIEWS_TEST, SQNET_SIZE + 1); % for a Matlab file. Each row is a trial with as many columns as there are cells +1
    for group = 1 : N_GROUPS
        for view = 1 : N_VIEWS_TEST
            infom(view + (group - 1) * N_VIEWS_TEST, 1) = group; % the stimulus number starting at 1
            Rate = EFRDATA(group, view, layer, :) * 100; % In VisNet the rates are 0-1, so this is more like spikes/s
            infom(view + (group - 1) * N_VIEWS_TEST, 2 : SQNET_SIZE + 1) = Rate;
        end
    end
    num_c = SQNET_SIZE; % number of cells
    num_s = N_GROUPS; % number of stimuli
    save('../VNanalysis/infos_matlab.mat', 'num_c', 'num_s', 'infom');
    infors_visnet(); % infors_visnet.m reads in data from ../VNanalysis/infos_matlab.mat

    % Read in the file plt produced by infors_visnet.m which has the stim-specific info for each cell
    StimSpecificInfo = dlmread('plt');                  % cell numbers 1 to SQNET_SIZE. col1=cell; col2=group; col3=info; col4=rate; col5=RateDifference
    StimSpecificInfoS = sortrows(StimSpecificInfo, -3); % sort by col 3, stim-specific info, descending order
    n = size(StimSpecificInfoS, 1);
    StimSpecificInfoS;
    % Comment in the following line to sort by the RateDifference
    %StimSpecificInfoS = sortrows(StimSpecificInfo, -5); % sort by col 5, RateDifference= 2*rate(s) - Sumrates, descending order
    
    counter = zeros(N_GROUPS, 1); % this will count for each group the number of cells already selected until 5 is reached
    for group = 1 : N_GROUPS
        i = 1;
        while counter(group) < NBESTCELLS && i <= n; % get 5 cells for each stimulus
            if (StimSpecificInfoS(i,2) == group) && (StimSpecificInfoS(i,5) > 5.0) % the latter checks that the rate is high.
                % StimSpecificInfoS(i,1) % the neuron number
                % StimSpecificInfoS(i,3) % the stimulus-specific info
                BESTCELLS(counter(group) +1 + (group-1)*NBESTCELLS, 1) = StimSpecificInfoS(i,1); % the neuron number, starts at 1
                BESTCELLS(counter(group) +1 + (group-1)*NBESTCELLS, 2) = StimSpecificInfoS(i,2); % the stimulus (i.e. group) number
                BESTCELLS(counter(group) +1 + (group-1)*NBESTCELLS, 3) = StimSpecificInfoS(i,3); % stimulus-specific info
                counter(group) = counter(group) + 1;
            end
            i = i + 1;
        end
    end
    format longG
    %BESTCELLS
    % Do a check that we have enough cells for the info analysis, useful for NET_SIZE=16
    for i = 1 : NBESTCELLS * N_GROUPS
        if BESTCELLS(i, 1) == 0
            BESTCELLS(i, 1) = 1; % set an artifical cell number
            BESTCELLS(i, 2) = floor((i-1) / NBESTCELLS) + 1; % set up the correct object / group number
            BESTCELLS(i, 3) = 0; % No information
        end
    end
    %BESTCELLS
    
    %BestCellsFilename = ['../VNanalysis/BestCells', int2str(LOGNUM), '.mat'];
    %save(BestCellsFilename, 'NBESTCELLS', 'BESTCELLS', '-v7.3'); % version 7.3 is good for matfile operations
    %save(BestCellsFilename, 'NBESTCELLS', 'BESTCELLS'); % For Octave
    
    format short
    
    % Get avg and max of best single cell info across groups
    AvgSingleCellInfo = 0.0;
    MaxSingleCellInfo = 0.0;
    for group = 1 : N_GROUPS
        AvgSingleCellInfo = AvgSingleCellInfo + BESTCELLS((group-1)*NBESTCELLS+1, 3);
        MaxSingleCellInfo = max(MaxSingleCellInfo, BESTCELLS((group-1)*NBESTCELLS+1, 3));
    end
    AvgSingleCellInfo = AvgSingleCellInfo / N_GROUPS;
    MaxSingleCellInfo

    % Now build infom.dat with the NBESTCELLS (5) best cells for each stimulus
    % infom.dat has col 1 for the stimulus number, then one column for each cell;
    % and one row for each transform of a stimulus.
    disp('Multiple Cell Information Analysis using inform_visnet');
    infom = zeros(N_GROUPS * N_VIEWS_TEST, N_GROUPS * NBESTCELLS + 1); % infom(stimulus_number, cell number)
    for i = 1 : NBESTCELLS * N_GROUPS    % the total number of cells in infom.dat
        cell = BESTCELLS(i, 1);          % the first cell has number 1 in this Matlab version
        infomRow = 1;
        for group = 1 : N_GROUPS
            for view = 1 : N_VIEWS_TEST
                infom(infomRow, 1) = group; % the first column is the stimulus number
                infom(infomRow, i + 1 ) = EFRDATA(group, view, layer, cell) * 100; % In VisNet the rates are 0-1, so this is more like spikes/s
                infomRow = infomRow + 1;
            end
        end
    end
  
    num_c = NBESTCELLS * N_GROUPS; % number of cells
    num_s = N_GROUPS; % number of stimuli
    save('../VNanalysis/infom_matlab.mat', 'num_c', 'num_s', 'infom');
    
    inform_visnet(); % this reads in the data from infom_matlab.dat
    
    fprintf('AvgSingleCellInfo= %5.4f bits\n', AvgSingleCellInfo);
    
    % correlation matrix between different test images in the firing of the specified layer
    N_GROUPS_TEST = N_GROUPS; % Set this up thus for now
    Rates = zeros(SQNET_SIZE, N_GROUPS_TEST * N_VIEWS_TEST * N_LOCS);
    counter = 1;
    for group = 1 : N_GROUPS_TEST
        for view = 1 : N_VIEWS_TEST
            EFRtmp = EFRDATA(group,view,layer,:,:); %object,view,loc,layer,R, C
            EFRtmp = reshape(EFRtmp, SQNET_SIZE, 1);
            Rates(:, counter) = EFRtmp;
            counter = counter + 1;
        end
    end 
    CorMat = zeros(N_GROUPS_TEST * N_VIEWS_TEST);
    for i = 1 : size(Rates, 2)
        for j = 1 : size(Rates, 2)
    %            distance(j,i) = norm(rates(i,:) - rates(j,:));
            %CorMat(j,i) = dot(Rates(:, i), Rates(:, j));
            R = corrcoef(Rates(:, i), Rates(:, j));
            CorMat(j,i) = R(1,2);
        end
    end
    figure(100)
    imagesc(CorMat);
    colormap(gray(256));
    title([' Correlation of firing for layer ', num2str(layer), ' between the transforms']);
    %title([' Correlation of firing for layer ', num2str(layer), ' between the transforms  min=', num2str(min(min(CorMat)))]);
    xlabel('Transform Number');
    ylabel('Transform Number');
    c = colorbar;
    ylabel(c, 'Correlation');
    colormap(gray(256));
    % min(min(CorMat))
    % max(max(CorMat))
    
    % Now calculate an object selectivity measure
    % i.e. how much the firing of the layer responds to all views of a given object
    % and to no views of other objects
    % A value of 1.0 indicates perfect object selectivity with complete transform invariance
    RespToObj = 0;
    RespToOtherObjs = 0;
    L = length(CorMat);
    NTransforms = N_VIEWS_TEST * N_LOCS; % number of transforms of each object
    for r = 1 : L % rows
        groupr = floor((r - 1) / NTransforms) + 1; % the group or object for this row
        for c = 1 : L % columns
            CorClipped = CorMat(r, c);
            if CorClipped < 0
                CorClipped = 0; % some correlations may be just < 0
            end
            groupc = floor((c - 1) / NTransforms) + 1; % the group or object for this column
            if groupr == groupc % if this is the same object
                RespToObj = RespToObj + CorClipped;
            else
                RespToOtherObjs = RespToOtherObjs + CorClipped;
            end
        end
    end
    RespToObj;
    RespToOtherObjs;
    ExpectedRespToObj = NTransforms * NTransforms * N_GROUPS; % This is what is expected if for each object the neurons respond to all views / tramsformsns of the object
    ObjectSelectivityExpected = RespToObj / (ExpectedRespToObj + RespToOtherObjs); % 1.0 would be perfect performance
    fprintf('ObjectSelectivity = %4.2f\n',    ObjectSelectivityExpected);
end

