TSP problem. Mixed algorithm

% Функция перебора всех пар ребер для задачи коммивояжера
% модернизированная версия метода отжига (лучшая замена)

% dist - матрица дистанций
% route - маршрут
% openclose - считать как закрытую или открытую задачу

%------------------------------------------------------------------------

function [tour, distance] = all_brute_best(dist,route,n)

    global_diff = -0.1;
    best_dist = inf;
    
    % если перестановки ребер не улучшают результат, выходим из цикла
    while global_diff < 0
        
        global_diff = 0;
        p1 = route(n);
        % Цикл по всем ребрам
        for i = 1:n-2
 
            t1 = p1;
            p1 = route(i);
            t2 = route(i+1);
            spd_var = dist(t1,p1);
            
            for j = i+2:n
                
                p2 = t2;
                t2 = route(j);
                diff = (dist(t1,p2) - dist(p2,t2))...
                    + (dist(p1,t2) - spd_var);
                
                % индексы на инвертацию пути с наибольшим эффектом
                if global_diff > diff
                    global_diff = diff;
                    % "лучшие" индексы для замены ребер
                    imin = i; jmin = j;
                end
                
            end
        end
        
        % Инвертируем путь
        if global_diff < 0
            route(imin:jmin-1) = route(jmin-1:-1:imin);
        end
        
    end

   %--------------------------------------------------------------------

    % считаем общее расстояние
    cur_dist = 0;
    for i = 1:n-1
        cur_dist = cur_dist + dist(route(i),route(i+1));
    end
    % считаем возврат в первый город
    cur_dist = cur_dist + dist(route(1),route(end));
    
   %--------------------------------------------------------------------

    % выбираем лучший маршрут
    if cur_dist < best_dist
        best_route = route;
        best_dist = cur_dist;
    end
    
    % переборный маршрут - его дистанция
    distance = best_dist;
    % сам маршрут
    tour = best_route;

    end

% Функция перебора всех пар ребер для задачи коммивояжера
% модернизированная версия метода отжига (каждая замена)

% dist - матрица дистанций
% route - маршрут
% openclose - считать как закрытую или открытую задачу

%-------------------------------------------------------------------------

function [tour, distance] = all_brute_first(dist,route,n)

    global_diff = -0.1;
    best_dist = inf;
    
    % если перестановки ребер не улучшают результат, выходим из цикла
    while global_diff < 0
        
        global_diff = 0;
        p1 = route(n);
        % Цикл по всем ребрам
        for i = 1:n-2
            
            t1 = p1;
            p1 = route(i);
            t2 = route(i+1);
            spd_var = dist(t1,p1);
            
            for j = i+2:n
                
                p2 = t2;
                t2 = route(j);
                diff = (dist(t1,p2) - dist(p2,t2))...
                    + (dist(p1,t2) - spd_var);
                
                % индексы на инвертацию пути с наибольшим эффектом
                if diff < 0
                    global_diff = diff;
                    % "лучшие" индексы для замены ребер
                    imin = i; jmin = j;
                    break
                end
            end
            
            % выход с цикла
            if diff < 0
                break
            end
            
        end
        
        % Инвертируем путь
        if global_diff < 0
            route(imin:jmin-1) = route(jmin-1:-1:imin);
        end
        
    end

   %--------------------------------------------------------------------

    % считаем общее расстояние
    cur_dist = 0;
    for i = 1:n-1
        cur_dist = cur_dist + dist(route(i),route(i+1));
    end
    % считаем возврат в первый город
    cur_dist = cur_dist + dist(route(1),route(end));
    
   %--------------------------------------------------------------------

    % выбираем лучший маршрут
    if cur_dist < best_dist
        best_route = route;
        best_dist = cur_dist;
    end
    
    % переборный маршрут - его дистанция
    distance = best_dist;
    % сам маршрут
    tour = best_route;

    end

parfor k = 1:n
    
    % временная матрица (выгорает при посещении вершины)
    dist_temp = dist;
    % маршрут по ближайшим городам
    route = zeros(1,n);
    % первая стартовая позиция
    next_point = start_route(k);
    % заполнение маршрута
    route(1) = next_point;
    
    % посещаем каждый ближайший город
    for e = 2:n
        % выгорание матрицы
        dist_temp(next_point,:) = nan;
        % поиск ближайшего оставшегося города
        [~,next_point] = min(dist_temp(:,next_point));
        % заполнение маршрута
        route(e) = next_point;
    end
    
    % првоеряем на уникальность и перестраиваем маршруты с 1-ого города
    if check_routes == true
        route = check_tour(route);
    end
    
    % массив машрутов
    arr_route_in(k,:) = route;
    
end

% (probability  several route)

% Функция, возвращаемая массив маршрутов методом рулетки среди топ вершин

% dist - матрица расстояний
% tour_sev - кол-во создаваемых маршрутов
% n - кол-во вершин
% tour_select - кол-во потенциальных вершин
% check_routes - проверка на уникальность и перестроение вершин
% p_change - вероятность перехода на альтернативный путь

%-------------------------------------------------------------------------

function [arr_several_route] = psr(dist,tour_sev,n,tour_select,cr,p_change)

    % возвращаемый массив (память)
    arr_several_route = zeros(tour_sev,n);
    
   parfor k = 1:tour_sev

        % временная матрица (выгорает при посещении вершины)
        dist_temp = dist;
        % маршрут по ближайшим городам
        route = zeros(1,n);
        % первая стартовая позиция (в данном случае случайная)
        next_point = randi([1,n]);
        % заполнение первой вершины маршрута
        route(1) = next_point;
        % счетчик
        count_var = 0;
        
        % цикл по посещениям вершин
        for e = 2:n
            % выгорание матрицы посещенной вершины
            dist_temp(next_point,:) = nan;
             % счетчик
            count_var = count_var + 1;           
            % вероятность выбора альтернативной вершины
            if rand(1) > p_change
                % выбираем безусловно кратчайший путь к k-ому городу
                % поиск ближайшего оставшегося города
                [~,next_point] = min(dist_temp(:,next_point));
            else
                % выбираем альтернативную вершину
                % массив дистанций перехода к следующей вершине (обратные)
                arr_next = 1./(dist_temp(:,next_point));
                % сайт хабрахабр транспонирование матрицы в массиве не 
                % воспринимает, поэтому переменная "habrahabr_vis" лишняя
                habrahabr_vis = (1:1:n)';
                % приписываем индексы к массиву
                arr_next = [arr_next, habrahabr_vis];
                % сортируем массив вместе с индексами
                arr_next_sort = sortrows(arr_next,1);
                % чтобы диапозон rand был корректен
                if (n - count_var) < tour_select
                    tour_select_fun = (n - count_var);
                else
                    tour_select_fun = tour_select;
                end
                % считаем вероятность перехода среди топ вершин
                P = arr_next_sort(1:tour_select_fun,:);
                P(:,1) = P(:,1)./sum(P(:,1));
                randNumb = rand;
                next_point = P(find(cumsum(P) >= randNumb, 1, 'first'),2);
            end
            % заполнение маршрута
            route(e) = next_point;
        end
        % првоеряем и перестраиваем маршруты с 1-ого города
        if cr == true
            route = check_tour(route);
        end
        % заполнение массива маршрутов (возвращаемый результат)
        arr_several_route(k,:) = route;

    end
    

% TSP() (на примере задачи коммивояжера)

%------------------------------------------------------------------------%

if exist('cities','var') == 1
    clearvars -except cities
    n = length(cities(:,1));
else
    clearvars
    intvar = inputdlg({'Кол-во вершин'}...
        ,'TSP',1);
    
    n = str2double(intvar(1));
    
    if isnan(n) || n < 0
        msgbox 'Введите корректные значения'
        return
    else
        cities = rand(n,2)*100;
    end
    
end

% начало счетчика времени алгоритма
tic

%-------------------------------Параметры--------------------------------%

% старт с определенного города (номер строки x и y координат)
start_tour = 1;

% проверка на уникальность маршрутов, перестроение маршрутов с первых
% городов (для проверки корректности алгоритма)
check_routes = false;

% использование метода рулетки при выбора определенного города
roulete_method = false;

% метод рулетки (кол-во маршрутов)
if roulete_method == true
    crrm = 1000;
end

% метод рулетки для топ "ts" кратчайших путей
top_roulete_method = true;

if top_roulete_method == true
    % кол-во создаваемых маршрутов методом "топ рулетки"
    ctrrm = 1000;
    % выбор топ самых привлекательных городов для перехода
    ts = 10;
    % вероятность принятия альтернативного решения
    p_change = 0.05;
end

%--------------------------------Память----------------------------------%

% матрица расстояний между городами
dist = zeros(n,n); 

% лучший маршрут (расстояние)
best_dist = inf;

% массив маршрутов на вход функции brute_edge
arr_route_in = zeros(n,n);

% бета - коэффициент (на доработку)
% b = 1;

%------------------------------------------------------------------------%

% матрица расстояний
for i = 1:n
    for j = 1:n
        dist(i,j) = sqrt((cities(i,1) - cities(j,1))^2 + ...
           (cities(i,2) - cities(j,2))^2);
    end
end

% округление матрицы дистанций
dist = round(dist);

% матрица обратных расстояний
returndist = 1./dist;

% генерируем стартовый маршрут (алгоритм ближайшего соседа)
start_route = linspace(1,n,n);

% стартовый маршрут с определенного города
if start_tour > n || start_tour < 1
    start_tour = 1;
end

% перестраиваем стартовый маршрут с определнного города "start_tour"
start_route(start_tour) = 1; start_route(1) = start_tour;

% создаем массив маршрутов методом ближайшего соседа
% стартуем с каждого последующего города
% parfor - параллельный цикл

parfor k = 1:n
    
    % временная матрица (выгорает при посещении вершины)
    dist_temp = dist;
    % маршрут по ближайшим городам
    route = zeros(1,n);
    % первая стартовая позиция
    next_point = start_route(k);
    % заполнение маршрута
    route(1) = next_point;
    
    % посещаем каждый ближайший город
    for e = 2:n
        % выгорание матрицы
        dist_temp(next_point,:) = nan;
        % поиск ближайшего оставшегося города
        [~,next_point] = min(dist_temp(:,next_point));
        % заполнение маршрута
        route(e) = next_point;
    end
    
    % првоеряем на уникальность и перестраиваем маршруты с 1-ого города
    if check_routes == true
        route = check_tour(route);
    end
    
    % массив машрутов
    arr_route_in(k,:) = route;
    
end

% создаем массив маршрутов методом топ ближайших соседей (методом рулетки)
% вероятность посещения вершины обратно пропорциональна длине ребра
if top_roulete_method == true
    arr_several = psr(dist,ctrrm,n,ts,check_routes,p_change);
    arr_route_in = [arr_route_in; arr_several];
end

% создаем массив маршрутов методом рулетки
% вероятность посещения вершины обратно пропорциональна длине ребра
if roulete_method == true
    % получаем массив маршрутов, выбором рулетки
    [randRoute] = prr(returndist,crrm,n,check_routes);
    % совмещаем два массива путей    
    arr_route_in = [arr_route_in;randRoute];
end

% % перестраиваем маршруты с первого города
% parfor r = 1:size(arr_route_in, 1);
%     temp_route = arr_route_in(r,:);
%     [~,numb_ix] = find(temp_route == 1);
%     if numb_ix ~= 1
%         temp_route = [temp_route(numb_ix:end),temp_route(1:numb_ix-1)];
%     elseif numb_ix == n
%         temp_route = fliplr(temp_route);
%     end
%     arr_route_in(r,:) = temp_route;
% end

% если есть, удаляем дубликаты маршрутов
arr_route_in = unique(arr_route_in,'rows');

%--------------------------------Память----------------------------------%
% массив дистанций с учетом оптимизаций (выход)
arr_tour_out = zeros(size(arr_route_in,1),n);
% массив результатов с учетом оптимизаций (выход)
arr_result = zeros(size(arr_route_in,1),1);
%------------------------------------------------------------------------%

% параллельное вычисление маршрутов (модернизированная оптимизация)

spdvar = size(arr_route_in, 1);

parfor i = 1:spdvar
%     [tour, distance] = all_brute_first(dist,arr_route_in(i,:),n);
    [tour, distance] = all_brute_best(dist,arr_route_in(i,:),n);
    arr_tour_out(i,:) = tour;
    arr_result(i,1) = distance;
end

% конец счетчика времени алгоритма
toc

% лучшая найденная дистанция
[best_length,best_indx] = min(arr_result);

% перестраиваем лучший маршрут с 1-ого города
best_route = arr_tour_out(best_indx,:);

% % среднее значение результатов
% mean(arr_result)

%-------------------------------Графика----------------------------------%

tsp_plot = figure(1);
set(tsp_plot,'name','TSP','numbertitle','off'...
    ,'color','w','menubar','none')
set(gcf,'position',[2391 410 476 409]);

opt_tour(:,[1,2]) = cities(best_route,[1,2]);

plot([opt_tour(:,1);opt_tour(1,1)],[opt_tour(:,2);opt_tour(1,2)],'-g.', ...
    'linewidth',1.5)
set(gca,'color','k')

title(['TSP best distance = ',num2str(round(best_length,2))])
xlabel('cities');
ylabel('cities')

%------------------------------------------------------------------------%

fprintf('%s %dn','Лучшая дистанция равна',round(best_length));

clearvars -except cities arr_result arr_tour_out arr_route_in ...
    n best_length opt_tour dist openclose