function CreateNet(receptiveFieldWidth, fieldWidth0, fieldWidth, Nfreq, Nang, Nsign) 
% Creates the connectivity of the synaptic matrices for VisNet, SYNMATL1-SYNMATL4
% receptiveFieldWidth, 
% fieldWidth0: the receptive field sizes for neurons in layer 1
% fieldWidth: the receptive field sizes for neurons in layers 2-4
% Nfreq, Nang, Nsign: the number of frequencies (4), angles (4) and signs (2) for the Gabor filtered images read into FILTIM
% If ScaleFactor = ScaleFactor / 2 below, use only the central 256x256 of the 512x512 filtered images,
% which is where the object is.
% This works well as translation invariance is not implemented by an image shift

global EFR EFRDATA NLAYERS NET_SIZE SQNET_SIZE N_FILES N_GROUPS N_VIEWS N_LOCS BESTCELLS NBESTCELLS ...
       NBESTCELLSINFORM NTRANSFORMS XY_VIEW XY_VIEW_TEST 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;


% Do layer 1 connections to the input filters of the images
    % There is a single .filt file for each image.
    % Each filtered image is 512x512 pixels in unsigned byte format,
    % stored by row, that is with columns varying most rapidly.
    % The order is highest frequencies first [1 2 4 8].
    % Within each freq the order is:
    % ang0_pos ang0_neg ang45_pos ang45_neg ang90_pos ang90_neg ang135_pos ang135_neg
    % For now, the byte pixels are held in a Matlab double array, to facilitate viewing + checking.
    % We set up the pointers to this image array, into which each image is loaded.
    % To implement translation invariance training, the aim is to slide the
    % pointer along this filter array just as in VisNetC, 
    % because using circshift would need a separate array for every filter,
    % and this would complicate the connectivity to SYNMATL1.
    
    tic;
    a1 = clock;
    layer0 = [FiltImSize FiltImSize];
    layer1 = [NET_SIZE NET_SIZE];
    fprintf('Connections for layer  1  ');
        
    ScaleFactor = FiltImSize / NET_SIZE;  % 512 / 32 = 16
    ScaleFactor = ScaleFactor / 2; % ETR for CreateNet2
    idx = 1; % This is the higher level neuron for which connections are being set up
    
    % Create a sparse matrix from a list of nonzero elements using the sparse function with five arguments:
    % S = sparse(i,j,v,m,n,nz)
    % i and j are vectors of row and column indices, respectively, for the nonzero elements of the matrix. 
    % v is a vector of nonzero values whose indices are specified by the corresponding (i,j) pairs. 
    % m is the row dimension of the resulting matrix, and n is the column dimension.
    % nz allocates space for nz nonzero elements. This syntax to allocates space for nonzero values to be filled in after construction.
    % SYNMATL1 = sparse(i, j, v, SqFiltImSize * Nfreq * Nang * Nsign, SQNET_SIZE, NSYN(1) * SQNET_SIZE); % as created below
    
    i = ones(NET_SIZE * NET_SIZE * NSYN(1), 1);  % row indices
    j= ones(NET_SIZE * NET_SIZE * NSYN(1), 1);   % column indices
    v = zeros(NET_SIZE * NET_SIZE * NSYN(1), 1); % the synaptic weights
    iter = 1;
    
    for col = 1 : NET_SIZE
        for row = 1 : NET_SIZE % build up the connectivity row by row - i.e each receptive field
            tmpRow = zeros([FiltImSize FiltImSize], 'double'); % care here
            connection = 1;
            randWeights = rand(1, NSYN(1)); % define the initial weights
            %randWeights = ones(1, NSYN(1)); % just to check the architecture
            randWeights = randWeights / norm(randWeights); % normalize the synaptic weights on each neuron
            while connection <= NSYN(1)
                randPos = mod(floor(randn(1,2) * receptiveFieldWidth(1)) + fieldWidth0, layer0) + 1;
                if tmpRow(randPos(1), randPos(2)) == 0
                    tmpRow(randPos(1), randPos(2)) = randWeights(connection);
                    connection = connection + 1;
                end
            end
            %tmpRow = circshift(tmpRow, [(row * ScaleFactor - FiltImSize/2 -1) (col * ScaleFactor - FiltImSize/2 -1)]);
            tmpRow = circshift(tmpRow, [(row * ScaleFactor - FiltImSize/4 -1) (col * ScaleFactor - FiltImSize/4 -1)]); % ETR: use this if ScaleFactor / 2
            % tmpRow now contains the col row indices for a single 512* 512 filtered image
            tmpRow = reshape(tmpRow, 1, FiltImSize * FiltImSize);
            synapseIndex = find(tmpRow);
            
            % The connections are now divided up across all the filters
            % We need to offset to the correct freq, angle and sign filtered component
            synindex = 1;
            for s = 1 : NSYNL1(1) % the highest spatial frequency
                filter = randi([1 8]); % one of the 8 filters, 4 ang * 2 signs
                index = synapseIndex(synindex) + ((filter - 1) * SqFiltImSize);
                i(iter) = index;
                j(iter) = idx;
                v(iter) = tmpRow(synapseIndex(synindex));
                synindex = synindex + 1;
                iter = iter + 1;
            end
            
            for s = 1 : NSYNL1(2) % the second highest spatial frequency
                filter = randi([1 8]);% one of the 8 filters, 4 ang * 2 signs
                index = synapseIndex(synindex) + ((filter - 1) * SqFiltImSize) + (SqFiltImSize * Nang * Nsign);
                i(iter) = index;
                j(iter) = idx;
                v(iter) = tmpRow(synapseIndex(synindex));
                synindex = synindex + 1;
                iter = iter + 1;
            end
            
            for s = 1 : NSYNL1(3) % the third highest spatial frequency
                filter = randi([1 8]);% one of the 8 filters, 4 ang * 2 signs
                index = synapseIndex(synindex) + ((filter - 1) * SqFiltImSize) + 2 * (SqFiltImSize * Nang * Nsign);
                i(iter) = index;
                j(iter) = idx;
                v(iter) = tmpRow(synapseIndex(synindex));
                synindex = synindex + 1;
                iter = iter + 1;
            end
            
            for s = 1 : NSYNL1(4) % the lowest spatial frequency
                filter = randi([1 8]);% one of the 8 filters, 4 ang * 2 signs
                index = synapseIndex(synindex) + ((filter - 1) * SqFiltImSize) + 3 * (SqFiltImSize * Nang * Nsign);
                i(iter) = index;
                j(iter) = idx;
                v(iter) = tmpRow(synapseIndex(synindex));
                synindex = synindex + 1;
                iter = iter + 1;
            end
            idx = idx + 1;
            
        end
    end
    SYNMATL1 = sparse(i, j, v, SqFiltImSize * Nfreq * Nang * Nsign, SQNET_SIZE, NSYN(1) * SQNET_SIZE);
    
    a2 = clock;
    disp(['Layer time: ', num2str(etime(a2, a1))]);
    
    % Do connections for layers 2-4
    for layer = 2 : 4
        
        a1 = clock;
        
        fprintf('Connections for layer %2d  ', layer);
        idx = 1; %  the neuron number
        for col = 1 : NET_SIZE
            for row = 1 : NET_SIZE
            % build up the connectivity row by row - i.e each receptive field
                tmpRow = zeros([NET_SIZE NET_SIZE], 'double'); % care here
                connection = 1;
                randWeights = rand(1, NSYN(layer)); % define the initial weights with a random distribution
                %randWeights = ones(1, NSYN(1)); % just to check the architecture
                randWeights = randWeights / norm(randWeights); % normalize the synaptic weights on each neuron
                while connection <= NSYN(layer)
                    randPos = mod(floor(randn(1,2) * receptiveFieldWidth(layer)) + fieldWidth, layer1) + 1;
                    if (tmpRow(randPos(1), randPos(2)) == 0)
                        tmpRow(randPos(1), randPos(2)) = randWeights(connection);
                        connection = connection + 1;
                    end
                end
                tmpRow = circshift(tmpRow, [(row - NET_SIZE/2 -1) (col - NET_SIZE/2 -1)]);
                if Display > 1  % Show the connections of a typical single neuron
                    if layer == 2 && col == NET_SIZE / 2 + 1 && row == NET_SIZE / 2 + 1
                    %if layer == 4 && col == 2 && row == NET_SIZE / 2 + 1
                        figure(99)
                        imagesc(tmpRow);
                        xlabel('X');
                        ylabel('Y');
                        title(['Connections of layer ', num2str(layer), ' neuron= ', num2str(idx)]);
                        colormap(gray);
                        drawnow;
                    end
                end    
                tmpRow = sparse(reshape(tmpRow, 1, NET_SIZE * NET_SIZE));
                if layer == 2
                    SYNMATL2(:, idx) = tmpRow;
                elseif layer == 3
                    SYNMATL3(:, idx) = tmpRow;
                elseif layer == 4
                    SYNMATL4(:, idx) = tmpRow;
                end
                idx = idx + 1;
            end
        end
        
        a2 = clock;
        disp(['Layer time: ', num2str(etime(a2, a1))]);
    
    end
    if NET_SIZE == 32
        save('../VNanalysis/netUntrained32.mat', 'NSYN', 'SQNET_SIZE', 'FiltImSize', 'SYNMATL1', 'SYNMATL2', 'SYNMATL3', 'SYNMATL4');
    elseif NET_SIZE == 64
        save('../VNanalysis/netUntrained64.mat', 'NSYN', 'SQNET_SIZE', 'FiltImSize', 'SYNMATL1', 'SYNMATL2', 'SYNMATL3', 'SYNMATL4');
    elseif NET_SIZE == 128
        save('../VNanalysis/netUntrained128.mat', 'NSYN', 'SQNET_SIZE', 'FiltImSize', 'SYNMATL1', 'SYNMATL2', 'SYNMATL3', 'SYNMATL4');
    elseif NET_SIZE == 256
        save('../VNanalysis/netUntrained256.mat', 'NSYN', 'SQNET_SIZE', 'FiltImSize', 'SYNMATL1', 'SYNMATL2', 'SYNMATL3', 'SYNMATL4');
    else
        save('../VNanalysis/netUntrained.mat', 'NSYN', 'SQNET_SIZE', 'FiltImSize', 'SYNMATL1', 'SYNMATL2', 'SYNMATL3', 'SYNMATL4');
    end
    
    fprintf('\nConnectivity creation time: %f sec\n', toc);
end
    

