function train(layer, locs, nepochs, TrainMode) 
    % This works for view+translation invariance training if there is > 1 view of each object in VNFileList
    % and for translation only training if there is only 1 view of each object in VNFileList
    % TrainMode=1 for translation invariance;   Not implemented
    % TrainMode=2 for view invariance or view and translation invariance
    % layer must be in range 1-4.
    % The code implements the synaptic update for layers 2-4 based on the
    % firing at t-1, so even if ETA is zero, some invariance learning
    % will occur. To check this, there is a commented trace update at time
    % t, and this results in no invariance learning.
    % N_LOCS and locs will always be 1 for this Matlab version
    % copyright E T Rolls, Oxford Centre for Computational Neuroscience
    
    % disp(['Matlab training layer ', int2str(layer)]);
    fprintf('Matlab training layer %1d    ', layer);
    
    tic

    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 ...
       NORM MAX_WEIGHTS

    % BestCells indices start at 1 in this Matlab version
    EFRsize = [NET_SIZE NET_SIZE];
    TraceDisp = zeros(nepochs * N_GROUPS * N_VIEWS, 3); % use this for theTrace, Rate and Trace-Rate for layer 4 display
    
    if layer == 1               % use Hebbian learning, no trace
        tic
        SYNMATL1_copy = SYNMATL1;              % a local copy instead of a global makes the Matlab code run faster
        lrnrate1 = LRNRATE(1);                 % similar
        for epoch = 1 : nepochs
            for group = 1 : N_GROUPS % the number of objects
                %Trace = zeros(SQNET_SIZE,1); % just in case we want to try a trace for layer 1, which is not the default
                %TransformCounter=0; % needed only for trace learning
                r = randperm(NTRANSFORMS);      % with permutation.
                for transform = r
                    view = transform;
                    im = view + (group - 1) * N_VIEWS;   % the image numbers start at 1, with N_VIEWS images or files per group
                    %fprintf('group=%d  view=%d im=%d IMLIST(im)=%s\n', group, view, im, char(IMLIST(im,1)));
                    image_name = char(IMLIST(im, 1));
                    image_dir = strcat(base_dir, image_name, '.filtered/');
                    %disp(strcat(image_dir, image_name,'.filt'));
                    %filename = strcat(image_dir, image_name,'.filt')
                    f_filt = fopen(strcat(image_dir, image_name,'.filt'), 'r');
                    FILTIM_copy = fread(f_filt); % this is to speed up Matlab
                    fclose(f_filt);
                    FILTIM_copy = FILTIM_copy';
                    if Display >= 2
                        i = 5; % should be the vertical pos filter
                        figure(10);
                        tmpFiltim = FILTIM_copy;
                        tmpFiltim = reshape(tmpFiltim, SqFiltImSize, 32);
                        tmp = tmpFiltim(:, i);
                        tmp =reshape(tmp, FiltImSize, FiltImSize);
                        imagesc(tmp');
                        colormap(gray);
                        title([' Filter ', num2str(i)]);
                        max(max(tmp));
                    end
                    EFR_copy = FILTIM_copy * SYNMATL1_copy; % activation

                    if  LATERALINHIB == 1
                        tmpD=EFR_copy;
                        tmpD=reshape(EFR_copy,NET_SIZE,NET_SIZE);
                        tmpD = conv2(tmpD, filter_layer_1, 'same');  % like Tinhibition
                        EFR_copy=reshape(tmpD,1,NET_SIZE*NET_SIZE);
                    end

                    mn = min(EFR_copy);        % convert to Rate using sparseness and threshold linear activation function
                    mx = max(EFR_copy);

                    EFR_copy = (EFR_copy - mn);
                    EFR_copy = EFR_copy / (mx - mn);
                    if Display >= 2
                        figure(1);
                        tmp = EFR_copy;
                        tmp =reshape(tmp, NET_SIZE, NET_SIZE);
                        imagesc(tmp');
                        colormap(gray);
                        title([' Actvn ', ' layer=', num2str(layer), '  group=', num2str(group), '  view=', num2str(view)]);
                        %disp(['layer=', num2str(layer), '  group=', num2str(group), '  view=', num2str(view)]);
                        drawnow
                    end
                    tmp = sort(EFR_copy, 'descend');
                    Threshold = tmp(floor(SQNET_SIZE * SPARSENESS(layer)));

                    if SIGMOID == 1 % following is for sigmoid
                        EFR_copy = EFR_copy - Threshold;% shift the array so that the value is 0 at p,
                        beta = -2.0 * BETA(layer);   % beta = -2¦Â(r-¦Á) 
                        for neuron = 1 : SQNET_SIZE
                            d = beta * EFR_copy(neuron);
                            if d > 70.0
                                EFR_copy(neuron) = 0.0;
                            elseif d < -50.0
                                EFR_copy(neuron) = 1.0;
                            else
                               EFR_copy(neuron) = 1.0 / (1.0 + exp(d));
                            end
                        end
                    else            %following is for linear threshold 
                        EFR_copy(EFR_copy < Threshold) = 0; % vectorize
                        % OR
%                         for neuron = 1 : SQNET_SIZE % apply a threshold linear activation function
%                             if EFR_copy(neuron) < Threshold
%                                 EFR_copy(neuron) = 0;
%                             end
%                         end
                    end

                    mn = min(EFR_copy);      
                    mx = max(EFR_copy);
                    EFR_copy = (EFR_copy - mn);
                    EFR_copy = EFR_copy / (mx - mn); % double check that rates are 0-1
                    
                    if Display >= 2
                        figure(2);
                        tmp = EFR_copy;
                        tmp =reshape(tmp, NET_SIZE, NET_SIZE);
                        imagesc(tmp'); 
                        colormap(gray);
                        title([' Rate ', ' layer=', num2str(layer), '  group=', num2str(group), '  view=', num2str(view)]);
                        disp(['layer=', num2str(layer), '  group=', num2str(group), '  view=', num2str(view)]);
                        drawnow
                        %pause
                    end


                    for i = 1 : SQNET_SIZE       % do Hebbian learning
                        SYNMATL1_col = SYNMATL1_copy(:, i); % column-wise operations speed up the operations on the sparse matrix
                        postRate = EFR_copy(i); % default Hebbian learning for layer 1
%                         postRate = Trace(i); % trace learning rule
%                         postRate = Trace(i) - (0.5 * EFR_copy(i)); % error correction rule
%                         postRate(postRate < 0) = 0;  % clip at 0. Necessary only for error correction. This will prevent negative weights, and run problems
                        synapseIndex = find(SYNMATL1_col);
                        preRates = FILTIM_copy(synapseIndex); % These rates will be in the range 0-255.
                        preRates = preRates ./ 255.0; % so scale so that they are now in the range 0-1. Alter the LRNRATE. VisNet2
                        if NORM == 1 %|| NORM == 2 % || NORM == 3 % disable Oja and MAX_WEIGHTS for layer 1 temporarily
                            dW = (preRates .* postRate) * lrnrate1;
                            %dW = sparse(dW); % does not speed matters for 32x32
                            temp1 = SYNMATL1_col(synapseIndex) + dW';
                            temp2 = temp1 ./ norm(temp1);
                        elseif NORM == 2 % Oja rule: disable for layer 1 ?
                            dW = (preRates  - (postRate .* SYNMATL1_col(synapseIndex)')) .* postRate .* lrnrate1; % Oja
                            temp1 = SYNMATL1_col(synapseIndex) + dW'; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                            if MAX_WEIGHTS(1) > 0 && SIGMOID == 1
                                temp1(temp1>MAX_WEIGHTS(1)) = MAX_WEIGHTS(1);
                            end
                            temp2 = temp1;
                        elseif NORM == 3 || NORM == 4 % Hebbian learning with Long-Term Depression (LTD) Eqn B.32 of Rolls 2023 Brain Computations and Connectivity
                            dW = (preRates  - SYNMATL1_col(synapseIndex)') .* postRate .* LRNRATE(layer); 
                            temp1 = SYNMATL1_col(synapseIndex) + dW'; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
%                             for j = 1 : length(synapseIndex)
%                                 if temp1(j) < 0
%                                     fprintf('Layer 1 weight is < 0');
%                                 end
%                             end
                            temp1(temp1<0) = 0; % make sure that no weights go negative
                            if MAX_WEIGHTS(1) > 0
                                temp1(temp1>MAX_WEIGHTS(1)) = MAX_WEIGHTS(1);
                            end
                            temp2 = temp1;
                         end
                        SYNMATL1_copy(:, i) = sparse(synapseIndex, ones(NSYN(layer), 1), temp2, length(FILTIM), 1);
                    end

%                     Trace = (ETA(layer) .* Trace) + ((1 - ETA(layer)) .* EFR_copy); % only now update the trace, so that we learn with the trace at t-1
                end
                if Display > 0
                    SYNMATL1 = SYNMATL1_copy; % copy back the local variable only necessary here for SisplaySynWeights
                    DisplaySynWeights(layer, epoch, group);
                    %pause
                end
            end % of group
        end % of epochs
        SYNMATL1 = SYNMATL1_copy; % copy back the local variable
        EFR(:, layer) = EFR_copy;
        fprintf('first layer runtime=%4.1f\n', toc);
        return; % end of layer 1
    else % so must be a higher layer
        FinalLayer = layer; % save the final layer
        for epoch = 1 : nepochs
            for group = 1 : N_GROUPS                 % the number of objects
                Trace = zeros(SQNET_SIZE,1);
                TransformCounter=0;
                r = randperm(NTRANSFORMS);      % with permutation. 
                %r = (1 : NTRANSFORMS);          % uncomment this ONLY IF a regular sequence is wanted for analyses
                r1 = randperm(NTRANSFORMS);      % get 5 random transforms to start the list while the trace builds up 
                r2 = r1(1 : 5);
                r = [r2 r]; % This is a way to learn with all N_VIEWS transforms, as the first 5 are not used for trace learning
                for transform = r
                    view = transform;
                    im = view + (group - 1) * N_VIEWS;    % the image numbers start at 1, with n_VIEWS images or files per group
                    image_name = char(IMLIST(im, 1));
                    image_dir = strcat(base_dir, image_name, '.filtered/');
                    %disp([strcat(image_dir, image_name,'.filt')]);
                    f_filt = fopen(strcat(image_dir, image_name,'.filt'), 'r');
                    FILTIM_copy = fread(f_filt);
                    fclose(f_filt);
                    %size(FILTIM)
                    FILTIM_copy = FILTIM_copy';
                    SYNMATL1_copy = SYNMATL1; % Local copies are faster than globals in Matlab
                    SYNMATL2_copy = SYNMATL2;
                    SYNMATL3_copy = SYNMATL3;
                    SYNMATL4_copy = SYNMATL4;

                    for layer = 1 : FinalLayer % get the firing rates first
                        switch layer
                            case 1
                            EFR(:, layer) = FILTIM_copy * SYNMATL1_copy; % activation
                            if  LATERALINHIB == 1 
                                tmpD=EFR(:,layer);
                                tmpD = reshape(EFR(:,layer), NET_SIZE, NET_SIZE);
                                tmpE = conv2(tmpD, filter_layer_1, 'same');  % like Tinhibition
                                EFR(:, layer) = reshape(tmpE, 1, NET_SIZE * NET_SIZE);
                            end
                            case 2
                            EFR(:, layer) = EFR(:, layer - 1)' * SYNMATL2_copy; % activation
                            if LATERALINHIB == 1
                                tmpD=EFR(:,layer);
                                tmpD = reshape(EFR(:,layer), NET_SIZE, NET_SIZE);
                                tmpE = conv2(tmpD, filter_layer_2, 'same');  % like Tinhibition
                                EFR(:, layer) = reshape(tmpE, 1, NET_SIZE * NET_SIZE);
                            end
                            case 3
                            EFR(:, layer) = EFR(:, layer - 1)' * SYNMATL3_copy; % activation
                            if LATERALINHIB == 1
                                tmpD=EFR(:,layer);
                                tmpD = reshape(EFR(:,layer), NET_SIZE, NET_SIZE);
                                tmpE = conv2(tmpD, filter_layer_3, 'same');  % like Tinhibition
                                EFR(:, layer) = reshape(tmpE, 1, NET_SIZE * NET_SIZE);
                            end
                            case 4
                            EFR(:, layer) = EFR(:, layer - 1)' * SYNMATL4_copy; % activation
                            if LATERALINHIB == 1
                                tmpD=EFR(:,layer);
                                tmpD = reshape(EFR(:,layer), NET_SIZE, NET_SIZE);
                                tmpE = conv2(tmpD, filter_layer_4, 'same');  % like Tinhibition
                                EFR(:, layer) = reshape(tmpE, 1, NET_SIZE * NET_SIZE);
                            end
                        end        

                        mn = min(EFR(:, layer));        % convert to Rate using sparseness and sigmoid or threshold linear activation function
                        mx = max(EFR(:, layer));
                        EFR(:, layer) = (EFR(:, layer) - mn);
                        EFR(:, layer) = EFR(:, layer) / (mx - mn); 
                        tmp = sort(EFR(:, layer), 'descend');
                        Threshold = tmp(floor(SQNET_SIZE * SPARSENESS(layer)));

                        if SIGMOID == 1  % following is for sigmoid activation function
                            EFR(:, layer) = EFR(:, layer) - Threshold;% shift the array so that the value is 0 at p,
                            beta = -2.0 * BETA(layer);   % beta = -2¦Â(r-¦Á)

                            for neuron = 1 : SQNET_SIZE
                                d = beta * EFR(neuron, layer);
                                if d > 70.0
                                    EFR(neuron, layer) = 0.0;
                                elseif d < -50.0;
                                    EFR(neuron, layer) = 1.0;
                                else
                                    EFR(neuron, layer) = 1.0 / (1.0 + exp(d));
                                end
                            end

                        else % threshold linear
                            for neuron = 1 : SQNET_SIZE
                                if EFR(neuron, layer) < Threshold
                                    EFR(neuron, layer) = 0;
                                end
                            end
                        end
                        mn = min(EFR(:, layer));  % set min to 0 and max to 1      
                        mx = max(EFR(:, layer));
                        EFR(:, layer) = (EFR(:, layer) - mn);
                        EFR(:, layer) = EFR(:, layer) / (mx - mn);
                    end % of layer: now we have the rate in the final layer for learning
                    
                    % the following line would be for the trace at time t to be used in the syn update, and should be commented out
                    %Trace = (ETA(layer) .* Trace) + ((1 - ETA(layer)) .* EFR(:, layer)); % only now update the trace, so that we learn with the trace at t-1
 
                    TransformCounter = TransformCounter + 1;
                    if (TransformCounter >= 5) % allow the trace to build up before learning
                        %fprintf('Layer=%d  mean(Trace)=%5.3f  mean(Firing)=%5.3f  mean(Trace-Firing)=%5.3f\n', layer, mean(Trace), mean(EFR(:, layer)), mean(Trace) - mean(EFR(:, layer)));
                        %fprintf('Layer=%d  Epoch=%d  mean(Trace-Firing)=%6.3f\n', layer, epoch, mean(Trace - EFR(:, layer)));
                        for i = 1 : SQNET_SIZE       % do Trace learning
                            postRate = Trace(i); % default trace learning rule
%                             postRate = Trace(i) - (0.5 * EFR(i, layer)); % error correction rule. Alter the 0.5 as useful.
%                             postRate(postRate < 0) = 0;  % clip at 0. Necessary only for error correction. This will prevent negative weights, and run problems
                            if layer == 2
                                SYNMATL2_col = SYNMATL2_copy(:, i); % column-wise operations speed up the operations on the sparse matrix
                                synapseIndex = find(SYNMATL2_col);
                                preRates = EFR(synapseIndex, layer - 1); % the firing rates of the presynaptic neurons connected to neuron i
                                if NORM == 1 % normalisation of the synaptic weight vectors to length=1
                                    dW = (preRates .* postRate) * LRNRATE(layer); % the change of the synaptic weights for neuron 1 to implement learning
                                    temp1 = SYNMATL2_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    temp2 = temp1 ./ norm(temp1); % normalise the synaptic weights for neuron i
                                elseif NORM == 2 % Oja rule
                                    dW = (preRates  - postRate .* SYNMATL2_col(synapseIndex)).* postRate * LRNRATE(layer); % the change of the synaptic weights for neuron 1 to implement learning
                                    temp1 = SYNMATL2_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(2) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(2)) = MAX_WEIGHTS(2);
                                    end
                                    temp2 = temp1;
                                elseif NORM == 3 % Hebbian learning with LTD eq B.32 of Rolls 2023 Brain Computations and Connectivity
                                    dW = (preRates  - SYNMATL2_col(synapseIndex)) .* postRate * LRNRATE(layer); 
                                    temp1 = SYNMATL2_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
%                                     for j = 1 : length(synapseIndex)
%                                         if temp1(j) < 0
%                                             fprintf('Layer %d weight is < 0\n', layer);
%                                         end
%                                     end
                                    temp1(temp1<0) = 0; % make sure that no weights go negative
                                    if MAX_WEIGHTS(2) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(2)) = MAX_WEIGHTS(2);
                                    end
                                    temp2 = temp1;
                                elseif NORM == 4 % MAX_WEIGHTS only with Hebbian learning
                                    dW = (preRates .* postRate) * LRNRATE(layer);
                                    temp1 = SYNMATL2_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(2) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(2)) = MAX_WEIGHTS(2);
                                    end
                                    temp2 = temp1;
                                end
                                SYNMATL2_copy(:, i) = sparse(synapseIndex, ones(NSYN(layer), 1), temp2, length(SYNMATL2), 1); % copy the changed weights back into the matrix of synaptic weights
                            elseif layer == 3
                                SYNMATL3_col = SYNMATL3_copy(:, i); % column-wise operations speed up the operations on the sparse matrix
                                synapseIndex = find(SYNMATL3_col);
                                preRates = EFR(synapseIndex, layer - 1);
                                if NORM == 1
                                    dW = (preRates .* postRate) * LRNRATE(layer);
                                    temp1 = SYNMATL3_col(synapseIndex) + dW;
                                    temp2 = temp1 ./ norm(temp1);
                                elseif NORM == 2 % Oja rule
                                    dW = (preRates  - postRate .* SYNMATL3_col(synapseIndex)).* postRate * LRNRATE(layer); % the change of the synaptic weights for neuron 1 to implement learning
                                    temp1 = SYNMATL3_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(3) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(3)) = MAX_WEIGHTS(3);
                                    end
                                    temp2 = temp1;
                                elseif NORM == 3 % Hebbian learning with LTD eq B.32 of Rolls 2023 Brain Computations and Connectivity
                                    dW = (preRates  - SYNMATL3_col(synapseIndex)) .* postRate * LRNRATE(layer); 
                                    temp1 = SYNMATL3_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(3) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(3)) = MAX_WEIGHTS(3);
                                    end
                                    temp2 = temp1;
                                elseif NORM == 4 % MAX_WEIGHTS only with Hebbian learning
                                    dW = (preRates .* postRate) * LRNRATE(layer);
                                    temp1 = SYNMATL3_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(3) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(3)) = MAX_WEIGHTS(3);
                                    end
                                    temp2 = temp1;
                                end
                                SYNMATL3_copy(:, i) = sparse(synapseIndex, ones(NSYN(layer), 1), temp2, length(SYNMATL3), 1);
                            elseif layer == 4
                                SYNMATL4_col = SYNMATL4_copy(:, i); % column-wise operations speed up the operations on the sparse matrix
                                synapseIndex = find(SYNMATL4_col);
                                preRates = EFR(synapseIndex, layer - 1);
                                if NORM == 1
                                    dW = (preRates .* postRate) * LRNRATE(layer);
                                    temp1 = SYNMATL4_col(synapseIndex) + dW;
                                    temp2 = temp1 ./ norm(temp1);
                                elseif NORM == 2 % Oja rule
                                    dW = (preRates  - postRate .* SYNMATL4_col(synapseIndex)).* postRate * LRNRATE(layer); % the change of the synaptic weights for neuron 1 to implement learning
                                    temp1 = SYNMATL4_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(4) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(4)) = MAX_WEIGHTS(4);
                                    end
                                    temp2 = temp1; 
                                elseif NORM == 3 % Hebbian learning with LTD eq B.32 of Rolls 2023 Brain Computations and Connectivity
                                    dW = (preRates  - SYNMATL4_col(synapseIndex)) .* postRate * LRNRATE(layer); 
                                    temp1 = SYNMATL4_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(4) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(4)) = MAX_WEIGHTS(4);
                                    end
                                    temp2 = temp1;
                                elseif NORM == 4 % MAX_WEIGHTS only with Hebbian learning
                                    dW = (preRates .* postRate) * LRNRATE(layer);
                                    temp1 = SYNMATL4_col(synapseIndex) + dW; % add the synaptic weight change vector to the vector of synaptic weights for neuron i
                                    if MAX_WEIGHTS(4) > 0 && SIGMOID == 1
                                        temp1(temp1>MAX_WEIGHTS(4)) = MAX_WEIGHTS(4);
                                    end
                                    temp2 = temp1; 
                                end
                                SYNMATL4_copy(:, i) = sparse(synapseIndex, ones(NSYN(layer), 1), temp2, length(SYNMATL4), 1);
%                                 if i == 2233 % analyse the trace and rate
%                                 for a neuron that learns to an Object:
%                                 use the gui to select a cell i
%                                     fprintf('Layer=%d  Group=%d View=%d  Epoch=%d  Trace=%6.3f Rate=%6.3f Trace-Rate=%6.3f  Trace-0.5*Rate=%6.3f\n', layer, group, view, epoch, Trace(i), EFR(i, layer), Trace(i) - EFR(i, layer), Trace(i) - 0.5*EFR(i, layer)), 
%                                     TraceDisp( (epoch-1)*N_GROUPS*N_VIEWS + (group-1)*N_VIEWS + view, 1) = Trace(i);
%                                     TraceDisp( (epoch-1)*N_GROUPS*N_VIEWS + (group-1)*N_VIEWS + view, 2) = EFR(i, layer);
%                                     TraceDisp( (epoch-1)*N_GROUPS*N_VIEWS + (group-1)*N_VIEWS + view, 3) = Trace(i) - 0.5*EFR(i, layer);
%                                 end
                            end
                        end % end of for i = 1 : SQNET_SIZE
                        SYNMATL2 = SYNMATL2_copy; % copy back the local variable
                        SYNMATL3 = SYNMATL3_copy; % copy back the local variable
                        SYNMATL4 = SYNMATL4_copy; % copy back the local variable
                    end % of TransformCounter
                    %ETA(layer)
                    Trace = (ETA(layer) .* Trace) + ((1 - ETA(layer)) .* EFR(:, layer)); % only now update the trace, so that we learn with the trace at t-1
                    %{
                    figure(1);
                    tmp = EFR(:, layer);
                    tmp =reshape(tmp, NET_SIZE, NET_SIZE);
                    imagesc(tmp');
                    colormap(gray);
                    title(' Rate ');
                    drawnow
                    figure(2);
                    tmp = Trace(:, layer);
                    tmp =reshape(tmp, NET_SIZE, NET_SIZE);
                    imagesc(tmp');                        
                    colormap(gray);
                    title(' Trace at t-1');
                    disp(['layer=', num2str(layer), '  group=', num2str(group), '  view=', num2str(view),  '  maxTrace=', num2str(max(Trace))]);
                    drawnow
                    %pause
                    %}
                end % of transform
                if Display > 0
                     DisplaySynWeights(layer, epoch, group);
                    %pause
                end % of if Display
            end % of group
        end % of epoch
    end % of layer ==
%     if layer == 4
%         figure(900)
%         plot(TraceDisp(:, 1));
%         hold on;
%         plot(TraceDisp(:, 2));
%         plot(TraceDisp(:, 3));
%         legend('Trace', 'Rate', 'Trace-Rate');
%     end
         
    fprintf('layer runtime=%4.1f\n', toc);
end