Прежде чем перейти к статье, хочу вам представить, экономическую онлайн игру Brave Knights, в которой вы можете играть и зарабатывать. Регистируйтесь, играйте и зарабатывайте!
Вот мы и добрались до написания программ.
В этой статье напишем скрипты для расчётов резистивного расстояния и для моделирования случайных блужданий. В качестве ЯП был выбран Octave (всё-таки математикой занимаемся).
И так, в прошлых двух частях мы познакомились с теорией данной темы. Настало время проверить теорию на практике.
Резистивное расстояние
Определим задачи для нашей программы:
Обеспечение универсальности (выбор файла с матрицей смежности графа, а также начального и конечного узлов).
Расчет резистивного расстояния между двумя произвольными вершинами графа.
Измерение времени выполнения программы.
Ввод и предварительная проверка данных
Начальный фрагмент определяет исходные данные задачи. Одним из ключевых аспектов являлась универсальность использования, что требовало реализации ввода имени файла для чтения и номеров узлов.
Здесь проверяется количество переданных параметров. В случае если их недостаточно, программа завершает работу с сообщением об ошибке. Далее происходит парсинг параметров: первый — имя файла, второй — начальный узел, третий — конечный узел. После происходит проверка введенных параметров на корректность. В случае неправильных данных вызывается ошибка с пояснением. Запускается таймер выполнения программы. В случае, если переданный файл уже обрабатывался, данные загрузятся из сохранённого файла.
args = argv();
if numel(args) < 3
error('Необходимо указать имя файла с матрицей смежности и',
'номера двух узлов.');
end
filename = args{1};
a = str2double(args{2});
b = str2double(args{3});
if isnan(a) || isnan(b) || a <= 0 || b <= 0
error('Номера узлов должны быть положительными целыми ',
'числами.');
endif
tic;
[~, name, ~] = fileparts(filename);
pseudo_inverse_filename = [name, '_L_plus.mat'];
if exist(pseudo_inverse_filename, 'file') == 2
load(pseudo_inverse_filename, 'L_plus');
disp(['Псевдообратная матрица Лапласа загружена из файла ',
pseudo_inverse_filename, '.']);
n = size(L_plus, 1);
Основные расчеты
Этот блок отвечает за тяжеловесные расчеты программы.
Считывается матрица смежности из указанного файла. Подсчитывается количество вершин. Вычисляются степени вершин, и создается диагональная матрица степеней. Затем вычисляется матрица Лапласиана. Происходит проверка корректности вычисления псевдообратной матрицы. Вычисляется псевдообратная матрица Лапласиана. Результат вычислений записывается в специальный файл. Считывание псевдообратной матрицы из этого файла позволит сократить время расчёта при очередном запуске программы с той же матрицей смежности.
else
A = dlmread(filename);
n = size(A, 1);
D = diag(sum(A, 2));
L = D - A;
if rank(L) < n-1
error('Матрица Лапласиана не имеет полного ранга.');
end
L_plus = pinv(L);
save(pseudo_inverse_filename, 'L_plus');
disp(['Псевдообратная матрица Лапласа вычислена и ',
'сохранена в файл ', pseudo_inverse_filename, '.']);
end
Расчет разницы потенциалов
Проверяется, что указанные узлы существуют в матрице. Создается результирующий вектор и производится расчет резистивного расстояния (эффективного сопротивления). Для оптимизации программы эти действия записаны в одну строчку.
if a > n || b > n
error('Номера узлов выходят за пределы размерности ',
'матрицы.');
endif
answer = (L_plus(a, a) - L_plus(b, a)) - (L_plus(a, b) -
L_plus(b, b));
Вывод результата
Остановка таймера выполнения и вывод результата:
elapsed_time = toc;
disp("===========================================================");
disp(['Резистивное расстояние между узлами ', num2str(a), ' и ',
num2str(b), ': ', num2str(answer)]);
disp(['Время выполнения программы: ', num2str(elapsed_time),
' секунд']);
Тестирование вычисления резистивного расстояния
Программа была протестирована на различных графах в три этапа:
Проверка правильности работы алгоритма на малых графах (до 5 вершин, заполненность 100%).
Тестирование на графах средних размеров (100-600 вершин, заполненность 90%).
Тестирование на больших графах (1000 и более вершин, заполненность 25-90%).
Данные для тестирования использовались без подгрузки предрасчитанной псевдообратной матрицы.
Результаты тестирования
Малые графы: алгоритм выдает решение менее чем за сотую долю секунды.
Средние графы: время выполнения менее полусекунды для графов с количеством вершин менее 500. При 500-600 вершинах расчет занимает 3.5 секунды.
Большие графы: для графа с 1000 вершинами и заполненностью 90% расчет занимает 6.5 секунд. Граф с 3000 вершинами и заполненностью 25% рассчитывается за 500 секунд.
![Зависимость времени выполнения программы от количества вершин в графе Зависимость времени выполнения программы от количества вершин в графе](https://habrastorage.org/r/w1560/getpro/habr/upload_files/078/398/61b/07839861bfd059b88df829781b0cd00b.png)
Основное время выполнения программы приходится на вычисление псевдообратной матрицы (сложность ). Поэтому, время выполнения увеличивается на три порядка при увеличении размера графа в 10 раз.
Реализация программы случайного блуждания
Аналогично использовался Octave. Но т.к. в этом скрипте достаточно большое количество переменных, думаю не лишним будет дать описание каждой.
Входные переменные
- матрица смежности графа.
- индекс начального узла.
- индекс конечного узла.
- количество тестов.
Выходные переменные
- вектор, содержащий время первого попадания в конечный узел для каждого теста.
- вектор, содержащий время полного обхода графа для каждого теста.
- вектор, содержащий время перемещения от начального узла к конечному и обратно для каждого теста.
- отсортированный вектор времени первого попадания и их количества.
- отсортированный вектор времени полного обхода графа и их количества.
- отсортированный вектор времени перемещения и их количества.
- среднее время первого попадания в конечный узел.
- среднее время полного обхода графа.
- среднее время перемещения от начального узла к конечному и обратно.
- эффективное сопротивление графа, рассчитанное как среднее время перемещения, деленное на удвоенное количество ребер.
Промежуточные переменные
- количество узлов в графе.
- количество ребер в графе, рассчитанное как половина суммы всех элементов матрицы смежности.
- текущий узел в процессе блуждания.
- булевый вектор, указывающий, были ли посещены узлы.
- количество шагов, сделанных в текущем тесте.
- булева переменная, указывающая, было ли достигнуто первое попадание в конечный узел.
- количество посещенных узлов.
- вектор соседних узлов для текущего узла.
- следующий узел для перехода.
- количество шагов для перемещения от начального узла к конечному и обратно в текущем тесте.
- путь к файлу с матрицей смежности.
- время, затраченное на выполнение тестов.
Скрипт случайных блужданий
В этом блоке инициализируются переменные для количества узлов и ребер в графе, векторов для хранения времени первого попадания, времени полного обхода, и времени перемещения.
function [fht, ct, mfht, mct, eff_res, mcmt] = random_walk(adj, start, end_, num_sim)
n = size(adj, 1);
m = sum(adj(:)) / 2;
fht = zeros(num_sim, 1);
ct = zeros(num_sim, 1);
cmt = zeros(num_sim, 1);
Следующией код выполняет цикл по числу тестов для выполнения случайных блужданий.
Этот блок начинает наш основной цикл, инициализируя текущий узел как начальный узел , а также массивы для отслеживания посещенных узлов и количества шагов. Также инициализируются переменные для отслеживания первого попадания в конечный узел и количества посещенных узлов.
for sim = 1:num_sim
curr = start;
visited = false(n, 1);
visited(start) = true;
steps = 0;
first_hit = false;
visit_count = 1;
Этот цикл продолжается до тех пор, пока не будут посещены все узлы графа. Внутри цикла определяется следующий узел для перехода случайным образом из соседей текущего узла.
while visit_count < n
neighbors = find(adj(curr, :));
next = neighbors(randi(length(neighbors)));
steps = steps + 1;
curr = next;
Если это первое попадание в конечный узел, сохраняется количество шагов до первого попадания.
if ~first_hit && curr == end_
fht(sim) = steps;
first_hit = true;
end
Если узел еще не был посещен, он помечается как посещенный, и увеличивается счетчик посещений.
if ~visited(curr)
visited(curr) = true;
visit_count = visit_count + 1;
end
end
После завершения обхода всех узлов записывается количество шагов.
ct(sim) = steps;
Затем рассчитывается время перемещения от начального узла к конечному и обратно. Это выполняется двумя циклами: один до конечного узла, и один обратно к начальному узлу.
cms = 0;
curr = start;
while curr ~= end_
neighbors = find(adj(curr, :));
next = neighbors(randi(length(neighbors)));
cms = cms + 1;
curr = next;
end
while curr ~= start
neighbors = find(adj(curr, :));
next = neighbors(randi(length(neighbors)));
cms = cms + 1;
curr = next;
end
cmt(sim) = cms;
end
Здесь рассчитывается среднее время первого попадания, среднее время обхода, среднее время перемещения и эффективное сопротивление.
mfht = mean(fht);
mct = mean(ct);
mcmt = mean(cmt);
eff_res = mcmt / (2 * m);
end
Для запуска функции random walk и получения результатов был разработан скрипт run random walk. Этот скрипт задаёт параметры графа, запускает функцию random walk и выводит результаты на экран.
Указание пути к файлу с матрицей смежности, начального узла, конечного узла и числа тестов. Затем чтение матрицы смежности из указанного файла с помощью функции и сохранение в переменную матрицы смежности.
args = argv();
if numel(args) < 3
error('Необходимо указать имя файла с матрицей смежности и',
'номера двух узлов.');
end
file_path = args{1};
start = str2double(args{2});
end_ = str2double(args{3});
num_sim = 1000;
adj = dlmread(file_path);
Выполнение функции случайного блуждания с заданными параметрами и измерение затраченного времени с помощью и
. Результаты сохраняются в соответствующих переменных, а время выполнения - в переменной
tic;
[fht, ct, mfht, mct, eff_res, mcmt] = random_walk(adj, start, end_, num_sim);
elapsed_time = toc;
Этот блок кода выводит основные результаты вычислений, такие как: среднее время первого попадания, среднее время прохода, среднее время обхода всего графа, эффективное сопротивление, время выполнения программы
fprintf('Среднее время первого попадания в вершину %d из' ,
'вершины %d: %f шага.\n', end_, start, mfht);
fprintf('Среднее время прохода из вершины %d в вершину %d и' ,
'обратно: %f шага.\n', start, end_, mcmt);
fprintf('Среднее время обхода всего графа: %f шага.\n', mct);
fprintf('Эффективное сопротивление: %f.\n', eff_res);
fprintf('Время выполнения программы: %f секунд.\n', elapsed_time);
Проведение тестов
В качестве примеров мы рассмотрим три графа. Первый граф будет тестовым для понимания методологии решения, вторым графом является Московское метро, третий - линейный граф.
Пример 1
Начальная вершина была выбрана как вершина 1, а конечная вершина как вершина 4. Количество тестов было установлено на 1000. Результаты тестов показали следующее:
Среднее время первого попадания в вершину 4 из вершины 1: 5.047200 шага.
Среднее время прохода из вершины 1 в вершину 4 и обратно: 10.004400 шага.
Среднее время обхода всего графа: 5.934400 шага.
Сопротивление: 1.000440.
Время выполнения программы: 7.940940 секунд.
Стоит заметить, что при стократном увеличении тестов, результаты изменятся незначительно.
Пример 2
В качестве второго примера была взята матрица Московского метро (442 вершины, 630 ребер), так как она позволяет проверить случайное блуждание в реальной транспортной сети. Начальная вершина была выбрана как вершина 1 (станция Новокосино), а конечная вершина как вершина 442 (станция Апрелевка). Количество тестов было установлено на 100000. Результаты тестов показали следующее:
Среднее время первого попадания в вершину 442 из вершины 1: 17159.75 шага.
Среднее время прохода из вершины 1 в вершину 442 и обратно: 20142.68 шага.
Среднее время обхода всего графа: 40777.67 шага.
Сопротивление: 15.99.
Время выполнения программы: 3194.03 секунд.
Резистивное расстояние по результатам выполнения программы для расчёта сопротивления с помощью законов Киргхофа равно 16.
Тогда из формулы (\ref{C_xy}) следует, что среднее время коммутирования равно:
Среднее время коммутирования по результатам случайных блужданий:
Результаты вычислений показывают, что среднее время прохождения из вершины 1 в вершину 442 и обратно почти точно удовлетворяет соотношению формулы (\ref{C_xy}).
Ошибка составлятет 0,01%.
Пример 3
Третьим примером стал линейный граф с количеством вершин 442 и количеством ребер 441. Количество тестов было установлено на 10000.
Среднее время первого попадания в вершину 442 из вершины 1: 191980.58 шага.
Среднее время прохода из вершины 1 в вершину 442 и обратно: 389058.45 шага.
Среднее время обхода всего графа: 191980.58.
Сопротивление: 441.11.
Время выполнения программы: 32729.89 секунд.
Резистивное расстояние по результатам выполнения программы для расчёта сопротивления с помощью законов Киргхофа равно 441.
Тогда из формулы (\ref{C_xy}) следует, что среднее время коммутирования равно:
Среднее время коммутирования по результатам случайных блужданий:
Результаты вычислений показывают, что среднее время прохождения из вершины 1 в вершину 442 и обратно достаточно точно удовлетворяет соотношению формулы (\ref{C_xy}).
Ошибка составлятет 0,03%.
Теоретические расчеты
Пример 1
Теперь с помощью теории рассчитаем количество шагов до попадания в вершину 4 из вершины 1. Составим матрицу степеней из матрицы смежности, используемой выше:
Из этих двух матриц составим матрицу Лапласа (Лапласиан), используя данную формулу:
Далее мы получаем собственные значения и собственные вектора матрицы Лапласа, чтобы их получить я воспользуюсь Octave:
laplacian_matrix = [2 -1 -1 0;
-1 3 -1 -1;
-1 -1 3 -1;
0 -1 -1 2];
[eigenvectors, eigenvalues_matrix] = eig(laplacian_matrix);
eigenvalues = diag(eigenvalues_matrix);
disp('Собственные значения:');
disp(eigenvalues);
disp('Собственные вектора:');
disp(eigenvectors);
Где eigenvalues – это , а eigenvectors – это
У собственных значений убираем нулевое значения и берем нулевую и третью строки из матрицы для расчета
. Получаем следующее:
Теперь нам нужны значения и
. У нас
уже известна, она равна 4, поэтому осталось найти
. Для этого находим сумму всех элементов матрицы смежности и делим их на 2.
В итоге получаем, что = 5 и
= 4. Получив эти значения, мы можем посчитать сопротивление
Теперь нужно посчитать commute time, для этого воспользуемся формулой из предыдущей статьи (см часть 2):
Нам осталось посчитать только hitting time, но так как наш граф является маленьким и симметричным, то у нас есть возможность воспользоваться данным вариантом формулы:
Для нашего маленького и симметричного графа, это равносильно что:
Следовательно равен:
Пример 2
Пример 1 был приведен для демонстрации методологии. Поскольку нахождение вручную данных графа Московского метро практически невозможно, воспользуемся программным методом. Для этого была написана программа, которая воспроизводит вычисления формул: ,
,
.
В этой части кода, мы загружаем матрицу смежности графа из файла. Затем определяем количество вершин и ребер.
args = argv();
if numel(args) < 3
error('Необходимо указать имя файла с матрицей смежности и',
'номера двух узлов.');
end
file_path = args{1};
adj_matrix = dlmread(file_path);
n = size(adj_matrix, 1);
m = sum(adj_matrix(:)) / 2;
Здесь мы вычисляем по формулам собственные значения и собственные векторы. Далее извлекаем собственные значения в вектор. И отбрасываем первое нулевое собственное значение и соответствующий ему собственный вектор (они не нужны для дальнейших вычислений).
degree_matrix = sum(adj_matrix, 2);
laplacian_matrix = diag(degree_matrix) - adj_matrix;
[eigenvectors, eigenvalues_matrix] = eig(laplacian_matrix);
eigenvalues = diag(eigenvalues_matrix);
eigenvalues_nonzero = eigenvalues(2:end);
eigenvectors_nonzero = eigenvectors(:, 2:end);
В этой части код инициализируем матрицу H нулями. Задаем максимальное количество итераций и допустимую ошибку.
H = zeros(n, n);
max_iter = 200000;
tol = 1e-16;
Здесь начинаем подсчет времени выполнения. Затем создаем список соседей для каждой вершины. В каждом элементе этого массива neighbors будут храниться массивы индексов соседей соответствующей вершины.
tic;
neighbors = cell(n, 1);
for i = 1:n
neighbors{i} = find(adj_matrix(i, :) == 1);
end
В этом блоке в каждой итерации обновляем значения H для всех вершин, используя их соседей. Затем каждые 10000 итераций выводим текущую ошибку. Если ошибка становится меньше заданного порога (tol), процесс останавливается. И в конце записываем время выполнения итерационного процесса.
for iteration = 1:max_iter
H_prev = H;
for u = 1:n
if degree_matrix(u) > 0
H(u, :) = 1 + (1 / degree_matrix(u)) * sum(H(neighbors{u}, :), 1);
end
H(u, u) = 0; % Условие H_ii = 0
end
if mod(iteration, 10000) == 0
disp(['Итерация: ', num2str(iteration), ' H - ', num2str(max(max(abs(H - H_prev))))]);
end
if max(max(abs(H - H_prev))) < tol
disp(['Сошлось после ', num2str(iteration), ' итераций']);
break;
end
end
elapsed_time = toc;
Здесь мы рассчитаем среднее время прохождения, время возвращения и сопротивление между заданными вершинами на основе полученных данных.
i = str2double(args{2});
j = str2double(args{3});
H_ij = H(i, j);
R_ij = sum((eigenvectors_nonzero(i, :) - eigenvectors_nonzero(j, :)).^2 ./ eigenvalues_nonzero');
C_ij = 2 * m * R_ij;
Выводим вычисленные результаты.
disp(['Среднее время прохода из вершины ', num2str(i), ' в вершину ', num2str(j), ': ', num2str(H_ij)]);
disp(['Среднее время прохода из вершины ', num2str(i), ' в вершину ', num2str(j), ' и обратно: ', num2str(C_ij)]);
disp(['Сопротивление: ', num2str(R_ij)]);
disp(['Время выполнения программы: ', num2str(elapsed_time), ' секунд']);
Результаты данных вычислений:
Сошлось после 293492 итераций
Среднее время прохождения из вершин 1 в вершину 442: 17085.752085552034
Среднее время прохождения из вершины 1 в вершину 442 и обратно: 20199.377213668093
Сопротивление: 16.03125175687944
Время выполнения программы: 2841.8106 секунд
Пример 3
Результаты вычислений для линейного графа:
Сошлось после 2016601 итераций
Среднее время прохождения из вершин 1 в вершину 442: 194480.99999830942
Среднее время прохождения из вершины 1 в вершину 442 и обратно: 388961.9999977657
Сопротивление: 440.9999999974668
Время выполнения программы: 8452.6199 секунд
Заключение
Результаты вычислительных экспериментов и теоретических расчетов для графа демонстрируют высокую степень согласованности, что свидетельствует о правильности применяемых моделей и алгоритмов. Тестирование проводилось на графе с четырьмя вершинами, на котором было установлено, что теоретические предсказания совпадают с практическими результатами. Кроме того, применение данных методов на лиенйном графе и на графе Московского метро показало, что алгоритмы могут работать не только в теоретических условиях, но и с реальной транспортной сетью. Среднее время первого попадания, прохождения и обхода всего графа, а также сопротивление, полученные в результате тестов, находятся в хорошем согласии с теоретическими расчетами. Эти три случая подтверждают, что предложенные модели и алгоритмы могут быть успешно применены для анализа и прогнозирования поведения более сложных и масштабных транспортных сетей.
Все скрипты доступны на GitHub.
Вот и написаны программы. Теперь можно применить статистический аппарат для анализа полученных данных. Но эта уже тема следующей главы.