function som_iris_ex4()
    % SOM on IRIS for different final neighborhood sizes
    % (port of your Python code to MATLAB)

    rng(7);  % reproducible random generator
    
    R = 12; 
    C = 12;          % bigger lattice
    M = R * C;
    iters = 4000;    % more iterations for bigger maps
    
    eta0 = 0.12;
    sigma0 = max(R, C) / 2;       % initial neighborhood (global ordering)
    sigma_final_list = [1 2 3];   % final sigmas to test
    
    % --- load and scale iris data (150x4, each feature to [0,1]) ---
    load fisheriris   % variables: meas (150x4), species
    X = meas;         % 150x4
    S = size(X, 1);
    
    Xmin = min(X, [], 1);
    Xmax = max(X, [], 1);
    X    = (X - Xmin) ./ (Xmax - Xmin);   % min-max scaling to [0,1]
    
    fprintf('loaded iris: %d samples, %d features\n', S, size(X, 2));
    fprintf('grid: %dx%d = %d neurons\n', R, C, M);
    fprintf('testing sigma_final: [');
    fprintf('%d ', sigma_final_list);
    fprintf(']\n');
    
    % grid coordinates (same idea as Python meshgrid + ravel)
    [gx, gy] = meshgrid(1:C, 1:R);   % gx,gy: R x C
    grid2d   = [gx(:), gy(:)];       % M x 2
    
    % --- figure for topographic error ---
    fig1 = figure('Color', 'w', 'Position', [100 100 800 500]);
    ax1 = axes('Parent', fig1);
    hold(ax1, 'on');
    xlabel(ax1, 'Iteration');
    ylabel(ax1, 'Topographic Error');
    title(ax1, sprintf('R=%d, C=%d — Effect of final \\sigma in {1,2,3}', R, C));
    grid(ax1, 'on');
    box(ax1, 'on');
    
    results = struct('sigma_final', {}, 'W', {}, 'topoErr', {}, 'topoErr_s', {});
    
    tau = iters;   % single time constant for whole run
    
    for k = 1:numel(sigma_final_list)
        sigma_final = sigma_final_list(k);
        fprintf('\n=== Testing sigma_final = %d ===\n', sigma_final);
        
        % weights in [0,1], same scale as X
        W = rand(M, size(X, 2));
        
        topoErr = zeros(iters, 1);
        
        for t = 1:iters
            % --- time-dependent sigma and eta ---
            sigma_t = sigma_final + (sigma0 - sigma_final) * exp(-t / tau);
            eta_t   = eta0 * exp(-t / tau);
            
            % --- random sample from X ---
            idx_sample = randi(S);
            x = X(idx_sample, :);   % 1 x D
            
            % --- BMU1, BMU2 ---
            dists = sum((W - x).^2, 2);  % M x 1
            [~, idx_sorted] = sort(dists, 'ascend');
            bmu1 = idx_sorted(1);
            bmu2 = idx_sorted(2);
            
            % --- Gaussian neighborhood around BMU1 in grid space ---
            dgrid2 = sum((grid2d - grid2d(bmu1, :)).^2, 2);  % M x 1
            h = exp(-dgrid2 / (2 * sigma_t^2));
            h = h / max(h);    % normalize
            
            % --- update weights ---
            % implicit expansion: h (M x 1) with (x - W) (M x D)
            W = W + eta_t * (h .* (x - W));
            
            % --- topographic error: 0 if BMU1,BMU2 are neighbors, else 1 ---
            topoErr(t) = ~are_neighbors(grid2d, bmu1, bmu2);
            
            if mod(t, 1000) == 0
                fprintf('  iteration %d/%d, sigma=%.3f, eta=%.5f\n', ...
                        t, iters, sigma_t, eta_t);
            end
        end
        
        % smoothing TE with moving average (okno = 50)
        win = 50;
        topoErr_s = movmean(topoErr, win);
        
        results(k).sigma_final = sigma_final;
        results(k).W           = W;
        results(k).topoErr     = topoErr;
        results(k).topoErr_s   = topoErr_s;
        
        % plot TE curve
        plot(ax1, 1:iters, topoErr_s, 'LineWidth', 1.5, ...
             'DisplayName', sprintf('\\sigma_{final}=%d', sigma_final));
    end
    
    legend(ax1, 'Location', 'northeast');
    
    % --- PCA projection and final lattices ---
    fig2 = figure('Color', 'w', 'Position', [100 100 1200 350]);
    
    for k = 1:numel(results)
        W = results(k).W;
        
        % PCA to 2D
        [coeff, score] = pca(W); %#ok<ASGLU>
        W2 = score(:, 1:2);   % M x 2
        
        WX = reshape(W2(:, 1), R, C);
        WY = reshape(W2(:, 2), R, C);
        
        ax = subplot(1, numel(results), k, 'Parent', fig2);
        hold(ax, 'on');
        axis(ax, 'equal');
        box(ax, 'on');
        title(ax, sprintf('\\sigma_{final}=%d', results(k).sigma_final));
        
        % horizontal lines (rows)
        for r = 1:R
            plot(ax, WX(r, :), WY(r, :), '-', 'LineWidth', 1.0);
        end
        
        % vertical lines (columns)
        for c = 1:C
            plot(ax, WX(:, c), WY(:, c), '-', 'LineWidth', 1.0);
        end
        
        % neuron markers
        plot(ax, W2(:, 1), W2(:, 2), 'ko', 'MarkerSize', 2, 'MarkerFaceColor', 'k');
    end
    
    sgtitle('Final Lattices Projected by PCA');
    
    % --- summary ---
    fprintf('\n== Summary (lower TE is better) ==\n');
    for k = 1:numel(results)
        topoErr = results(k).topoErr;
        last500 = topoErr(end-499:end);
        te_mean = mean(last500);
        fprintf('sigma_final=%d -> mean TE (last 500 iters): %.4f\n', ...
                results(k).sigma_final, te_mean);
    end
end

% =======================================================================
% local helpers
% =======================================================================

function tf = are_neighbors(grid2d, i, j)
    % odleglosc Chebyshev = 1 (8-neighborhood)
    % grid2d: M x 2, i,j: indices of neurons
    di = abs(grid2d(i,1) - grid2d(j,1));
    dj = abs(grid2d(i,2) - grid2d(j,2));
    tf = max(di, dj) == 1;   % true if neighbors
end