Category Archives: Lessons

Проверка попадания точки в симплекс

Предположим, есть d точек \{a_i\}_1^d в пространстве \mathbb R^{d-1}.
Требуется проверить, принадлежит ли произвольное p \in \mathbb R^{d-1} выпуклой оболочке \{a_i\}_1^d.

Для простоты предположим, что объём \mathrm{co}\,\{a_i\}_1^d не равен нулю, тогда это симплекс, вершинами которого выступают \{a_i\}_1^d. Возможное технически простое решение данной задачи подсказывает теорема Каратеодори о выпуклой оболочке. В соответствии с этой теоремой, если p \in \mathrm{co}\,\{a_i\}_1^d, то p можно представить в виде выпуклой линейной комбинации: p = \sum_{i=1}^d{c_i a_i}, где \{c_i\}_1^d есть некоторый набор неотрицательных чисел такой, что \sum_{i=1}^d{c_i} = 1.

Теперь исходную задачу можно представить в виде задачи решения системы линейных уравнений:

\begin{bmatrix}  a_{11} & a_{21} & \cdots & a_{d1} \\  \vdots & \vdots & \cdots & \vdots \\  a_{1,d-1} & a_{2,d-1} & \cdots & a_{d,d-1} \\  1 & 1 & \cdots & 1  \end{bmatrix}  \begin{bmatrix}  c_1 \\ \vdots \\ c_{d-1} \\ c_d  \end{bmatrix} =  \begin{bmatrix}  p_1 \\ \vdots \\ p_{d-1} \\ 1  \end{bmatrix}.

Решив данную систему относительно \{c_i\}_1^d, можно проверить, что все они неотрицательные. Если это так, то точка p лежит в \mathrm{co}\,\{a_i\}_1^d.

Также можно заметить, что если c_i = 0, то p лежит на грани симплекса, противолежащей вершине a_i.
Если c_i < 0, то p лежит в наружном полупространстве, находящемся за этой гранью.
Если же все c_i > 0, то p лежит во внутренности симплекса.

Реклама

Об использовании Visual Studio Code для учебных примеров

Visual Studio Code — редактор на основе платформы Electron, который можно использовать как минималистичную IDE. Мне требовалась поддержка C++17 и переносимость всего «на флешке». Что меня приятно удивило, так это «шустрость» VSCode.

Итак, порядок действий, предпринятых мной был таков:

  1. Скачать VSCode в архиве (в моём случае это была версия 1.27.1 для Windows/x64). Распаковать его, допустим, в папку vscode. Внимание! Путь к vscode не должен содержать не-ASCII символы и пробелы.
  2. Создать пустую папку vscode/data. Там будут плагины и общие настройки. При обновлении VSCode можно просто распаковать новую версию и переместить туда свою папку data из старой.
  3. Распаковать сборку MinGW рядом с vscode, допустим, в папку mingw (mingw/bin содержит g++.exe). В качестве таковой я брал сборку отсюда, там уже есть Boost и ряд других библиотек. Я к ним ещё добавляю Eigen и кое-что своё. Добавляю, просто складывая соответствующие папки в mingw/include.
  4. Создать рядом с vscode пустую папку для примеров. Допустим, cpp.
  5. Запустить VSCode (vscode/code.exe) и установить плагин C/C++ for Visual Studio Code от Microsoft (с Microsoft Intellisense).
  6. Открыть в VSCode папку cpp.
  7. Настройки создаются в виде JSON-файлов в папке cpp/.vscode. Ниже я опишу возможное их содержимое. Их можно создать непосредственно в VSCode, заполнить в нём же, и затем перезапустить VSCode.
  8. После того, как файлы настроек созданы, и VSCode перезапущен, следует протестировать его. Создать cpp/test.cpp с произвольным содержимым на C++ (например, helloworld). Terminal/Run build task… или Ctrl+Shift+B должно запускать сборку и завершаться без ошибок. F5 запускает сборку и результат в отладчике (gdb). Заголовочные файлы должны успешно открываться (правой кнопкой мыши и Go to definition). Intellisense должен работать адекватно.
  9. Полученный «пакет» в составе vscode, mingw, cpp можно упаковать в архив и развернуть на другой системе или просто скопировать на внешний носитель и запускать с него.

Итак, файлы настроек.

Плагин C/C++ использует файл c_cpp_properties.json.
Его содержимое в моём случае:

{
    "configurations": [
        {
            "name": "MinGW",
            "intelliSenseMode": "gcc-x64",
            "compilerPath": "${workspaceFolder}\\..\\mingw\\bin\\g++.exe",
            "cppStandard": "c++17",
            "cStandard": "c11",
            "includePath": [
                "${workspaceFolder}",
                "${workspaceFolder}\\..\\mingw\\include"
            ],
            "defines": [
                "__cplusplus=201703L"
            ]
        }
    ],
    "version": 4
}

Строго говоря, дефайн стандартного макроса __cplusplus, может, и не понадобится, просто у меня в начале без него плагин не видел правильное содержимое стандартных заголовков.

Для сборки требуется создать «задачу сборки». В простейшем случае это компиляция одной единицы трансляции (cpp-файла) сразу в исполняемый файл (в моём случае всегда один и тот же с названием lab.exe).
Описание задач сборки располагается в файле tasks.json. Его содержимое в моём случае:

{
    // See https://go.microsoft.com/fwlink/?LinkId=733558
    // for the documentation about the tasks.json format
    "version": "2.0.0",
    "tasks": [
        {
            "label": "build lab",
            "type": "shell",
            "command": "${workspaceFolder}\\..\\mingw\\bin\\g++",
            "args": [
                "-g", "-Og", // "debug" build
                //"-O3",
                "${file}", 
                "-march=native", "-Wall", "-pedantic",
                "-std=c++17",
                "-olab.exe"
            ],
            "group": {
                "kind": "build",
                "isDefault": true
            }
        }
    ]
}

Наконец, описание запуска отладчика по F5 располагается в launch.json.
Его содержимое в моём случае:

{
    // Use IntelliSense to learn about possible attributes.
    // Hover to view descriptions of existing attributes.
    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
    "version": "0.2.0",
    "configurations": [
        {
            "name": "(gdb) Launch",
            "type": "cppdbg",
            "request": "launch",
            "program": "${workspaceFolder}\\lab.exe",
            "args": [],
            "stopAtEntry": false,
            "cwd": "${workspaceFolder}",
            "environment": [],
            "externalConsole": true,
            "MIMode": "gdb",
            "miDebuggerPath": "${workspaceFolder}\\..\\mingw\\bin\\gdb.exe",
            "setupCommands": [
                {
                    "description": "Enable pretty-printing for gdb",
                    "text": "-enable-pretty-printing",
                    "ignoreFailures": true
                },
                {
                    "description": "Enable break on exception for gdb",
                    "text": "catch throw",
                    "ignoreFailures": true              
                }
            ],
            "preLaunchTask": "build lab"
        }
    ]
}

Вторая команда catch throw активирует перехват исключений отладчиком gdb.
Вообще, к сожалению, я не могу сказать, что с отладкой всё хорошо: есть ряд странностей. Например, точка останова, поставленная в одном месте, может сработать ниже (как будто поставленная в другом месте — может быть, это можно исправить, если заменить -Og на -O0). Точку останова нельзя поставить внутри шаблона (не срабатывает вообще). Не даёт ставить новые точки останова во время отладки. Если во время отладки нажать паузу, то продолжение может не срабатывать (останется только остановить исполнение).

Для меня проблемы с отладчиком не являются критичными, так как я привык к отладке через проверку условий и протоколирование. Но тем, кто привык к отладчику в «нормальной» Visual Studio такое, конечно, будет не по нраву.

UPD. Если всё это запускать на Windows 10, то в качестве терминала может быть по умолчанию выбран PowerShell. Почему-то он «не понимает» относительные пути через .. и не может запустить компилятор. В этом случае рекомендуется заменить его на cmd.exe (в VSCode можно выбрать терминал).

UPD2. Замена флага оптимизации -Og на -O0 всё же помогает отладчику, по крайней мере, точки останова не «убегают» с тех мест, куда они были поставлены.

Граница Парето в 2D: STL-версия

Почти пять лет назад я описал алгоритм вычисления границы Парето множества точек на плоскости.

С тех пор я использовал данную задачу как задачу по программированию на первом курсе (на алгоритмы STL). Всё дело в том, что этот алгоритм очень просто записывается с помощью STL. Требуется применить два алгоритма (sort и unique) и адаптер компаратора (not_fn или самописный, например, переставляющий аргументы местами). Бывают простенькие задачи, к которым невозможно адекватно применить STL без написания своих компонент, а бывает и наоборот. Итак, код.

#include <algorithm>
#include <functional>

template 
  <
    class RanIt, 
    class Comp1, 
    class Comp2
  >
RanIt remove_dominated
  (
    RanIt from,
    RanIt to,
    Comp1 compare1,
    Comp2 compare2
  )
{
  std::sort(from, to, std::not_fn(compare1));
  return std::unique(from, to, std::not_fn(compare2));
}

////////////////////////////
#include <string>
#include <vector>
#include <iostream>
using namespace std;

bool test_remove_dominated()
{
  vector<string> input
  {
    "look",
    "after",
    "the",
    "raven",
    "star"
  };
  
  input.erase(
    remove_dominated(
      input.begin(), input.end(),
      std::less<>{},
      [](auto & a, auto & b)
      {
        return a.size() < b.size();
      }),
    input.end()
    );
    
  sort(input.begin(), input.end());
  return input == vector<string>
    {
      "raven",
      "star",
      "the"
    };
}

int main()
{
  cout << test_remove_dominated() << '\n';
  return 0;
}

В примере два критерия сравнения строк: лексикографический и по длине. Так «raven» доминирует «look» (и лексикографически, и по длине), но не «star».

О вычислении матричной экспоненты

1. Экспонента квадратной матрицы \mathbf A формально определяется как сумма ряда:

e^{\mathbf A} \triangleq \mathbf I + \sum\limits_{i=1}^{\infty} \dfrac{\mathbf A^i}{i!},

где \mathbf I есть единичная матрица соответствующего размера.

Из этого определения несложно получить ряд свойств экспоненты, некоторые из которых перечислены ниже:

  1. Экспонента нулевой матрицы есть единичная матрица.
  2. Транспонирование (а также эрмитово сопряжение) перестановочно с экспоненцированием: e^{\mathbf A^\top} = (e^{\mathbf A})^\top.
  3. Если \mathbf A симметрична (эрмитова), то и e^{\mathbf A} симметрична (эрмитова).
  4. Если \mathbf {A B} = \mathbf {B A}, то e^{\mathbf A} \mathbf B = \mathbf B e^{\mathbf A} и e^{\mathbf A} e^{\mathbf B} = e^{\mathbf B} e^{\mathbf A} = e^{\mathbf A + \mathbf B}. В частности, e^{a \mathbf A} e^{b \mathbf A} = e^{(a + b)\mathbf A}.
  5. Матрица e^{\mathbf A} всегда обратима: e^{\mathbf A} e^{-\mathbf A} = \mathbf I.
  6. Для любой обратимой матрицы \mathbf P верно e^{\mathbf{P A P}^{-1}} = \mathbf P e^{\mathbf A} \mathbf P^{-1}.
  7. Экспонента диагональной матрицы \mathbf A = \mathrm{diag}(a_{ii}) также является диагональной матрицей e^{\mathbf A} = \mathrm{diag}(e^{a_{ii}}).
  8. Последнее свойство обобщается на блочно-диагональные матрицы. Пусть \mathbf A = \mathrm{diag}(\mathbf A_k), где каждая \mathbf A_k есть квадратная матрица. Тогда экспонента также является блочно-диагональной матрицей: e^{\mathbf A} = \mathrm{diag}(e^{\mathbf A_k}).

Таким образом, если матрица \mathbf A диагонализируема, т.е. если найдётся такая матрица \mathbf P, что \mathbf A = \mathbf {P D P}^{-1}, где \mathbf D = \mathrm{diag}(\lambda_i) — диагональная матрица (матрица собственных значений), то после выполнения диагонализации (eigendecomposition) \mathbf A вычислить её экспоненту не составит труда.

К сожалению, не все матрицы диагонализируемы, хотя для матриц над полем комплексных чисел множество недиагонализируемых матриц есть множество меры нуль. Поэтому, кстати, даже экспоненту вещественнозначной матрицы может быть удобнее вычислять в комплексных числах. Те, кто хорошо учил линейную алгебру, вспомнят о «жордановой форме», которая существует для любой матрицы и представляет собой как раз блочно-диагональную матрицу, составленную из блоков, называемых «жордановы клетки». Это означает, что если у нас есть представление \mathbf A = \mathbf{P J P}^{-1}, где \mathbf J = \mathrm{diag}(\mathbf J_k), и если мы умеем вычислять экспоненту жордановой клетки, то мы опять же без проблем вычислим экспоненту исходной матрицы. Ниже я коснусь вопроса вычисления экспоненты жордановой клетки.

Особняком стоит свойство, известное как «формула Якоби». Его можно применять, например, для проверки качества вычисления матричной экспоненты:

\det e^{\mathbf A} = e^{\mathrm{tr}\,\mathbf A}.

2. Вычисление экспоненты прямым суммированием ряда традиционно считается «плохим методом». На самом деле, главный его недостаток заключается в низкой скорости вычислений, но не в точности. Важно организовать этот процесс правильно! В своей практике я применял следующие приёмы:

  • Компенсационное суммирование. Например, алгоритм Ноймайера.
  • Вместо вычисления очередной степени \mathbf A^i домножением, я вычислял степень каждый раз «заново» (до 2\log_2{i} + 1 матричных умножений).
  • Единичная матрица добавляется в конце (опционально: не нужно её добавлять, если вычисляется e^{\mathbf A} - \mathbf I ).
  • Суммирование продолжается до тех пор, пока очередное слагаемое не обнулится из-за слишком большого делителя-факториала (на двойной точности обычно достаточно несколько десятков слагаемых).

Я не претендую на то, что это лучший способ вычислять данную сумму. Но исходя из моего (ограниченного) опыта, этот способ достаточно надёжен и сравнительно прост в реализации. В моей практике собственно скорость вычисления экспоненты не играла большой роли. А точность была важна.

3. Экспонента жордановой клетки. Жорданова клетка \mathbf J (в случае комплексных чисел) есть квадратная матрица, на главной диагонали которой стоит некое число \lambda (любое значение, но одно во всех позициях диагонали), на первой наддиагонали стоят единицы, а все прочие элементы — нули:

\mathbf J = \begin{bmatrix}  \lambda & 1 & 0 & \ldots & 0 \\  0 & \lambda & 1 & \ldots & 0 \\  & & \cdots & & \\  0 & 0 & \ldots & \lambda & 1 \\  0 & 0 & \ldots & 0 & \lambda  \end{bmatrix}.

Заметим, что \mathbf J = \mathbf \Lambda + \mathbf N, где \mathbf \Lambda — шаровая матрица (диагональная с одинаковыми значениями на диагонали), а \mathbf N — нильпотентна, поскольку \mathbf N^n = \mathbf 0 (здесь и ниже n — размер \mathbf J). Более того, легко заметить, что домножение \mathbf N^i на \mathbf N попросту сдвигает единичную наддиагональ на одну позицию вверх. Поэтому имеем

e^{\mathbf N} = \mathbf I + \sum\limits_{i=1}^{n-1}\dfrac{\mathbf N^i}{i!} = \begin{bmatrix}  1 & 1 & 1/2 & 1/6 & \ldots & 1/(n - 1)! \\  0 & 1 & 1 & 1/2 & \ldots & 1/(n-2)! \\  & & & \cdots & & \\  0 & 0 & \ldots & 0 & 1 & 1 \\  0 & 0 & \ldots & 0 & 0 & 1  \end{bmatrix}.

Ну и по свойствам экспоненты и шаровых матриц имеем

e^{\mathbf J} = e^{\mathbf \Lambda + \mathbf N} = e^{\mathbf \Lambda} e^{\mathbf N} = e^\lambda e^{\mathbf N}.

В случае наличия скалярного множителя t:

e^{\mathbf J t} = e^{\mathbf \Lambda t + \mathbf N t} = e^{\lambda t} e^{\mathbf N t},

где

e^{\mathbf N t} = \begin{bmatrix}  1 & t & t^2/2 & t^3/6 & \ldots & t^{n-1}/(n - 1)! \\  0 & 1 & t & t^2/2 & \ldots & t^{n-2}/(n-2)! \\  & & & \cdots & & \\  0 & 0 & \ldots & 0 & 1 & t \\  0 & 0 & \ldots & 0 & 0 & 1  \end{bmatrix}.

К сожалению, вычисление жордановой формы как численный метод обладает врождённым недостатком: чувствительностью к сколь угодно малым изменениям исходной матрицы. Из чего можно сделать вывод, что на практике может быть разумно вносить в недиагонализируемую матрицу малые возмущения с целью получить диагонализируемую матрицу.

4. О вычислении матричной экспоненты (и иных функций, определяемых как сумма ряда) с помощью комбинирования разложения Шура, рекуррентных соотношений Сильвестра и суммировании ряда см. Davies, Higham. A Schur-Parlett Algorithm for Computing Matrix Functions.

Об аналитическом решении линейных ОДУ

1. Пусть дана однородная система (задача Коши):

\dot{\mathbf x} = \mathbf A \mathbf x, \quad \mathbf x(t_0) = \mathbf x_0,

где \mathbf x(t) \in \mathbb R^n — фазовый вектор системы, \mathbf A — постоянная n\times n-матрица над полем \mathbb R. Цель — найти функцию \mathbf x(t), t \geqslant t_0, удовлетворяющую указанной системе. Запись \dot{\mathbf x} обозначает производную \mathbf x'(t).

Как и в случае n = 1 решение можно получить через экспоненту.
В нашем случае — матричную экспоненту, определяемую формально как сумму ряда. Вычисление матричной экспоненты, особенно для больших n, есть отдельная большая тема, которой в данной статье я не намерен касаться подробно. Ниже приведены формулы для прямого вычисления матричной экспоненты в случае n = 2.

Итак, сразу решение:

\mathbf x(t) = e^{\mathbf A(t - t_0)} \mathbf x_0.

Заметим, что, как и в скалярном случае, мы можем продифференцировать функцию e^{\mathbf A t} по t как матрицу покомпонентно, получив:

(e^{\mathbf A t})' = e^{\mathbf A t} \mathbf A = \mathbf A e^{\mathbf A t}.

\mathbf A, \mathbf A^{-1} и e^{\mathbf A t} являются перестановочными матрицами, что легко показать, раскрыв экспоненту в сумму ряда.

2. Пусть теперь дана неоднородная система:

\dot{\mathbf x} = \mathbf A \mathbf x + \mathbf f, \quad \mathbf x(t_0) = \mathbf x_0,

где \mathbf f \colon \mathbb R \to \mathbb R^n. Выполнив замену

\mathbf z(t) = e^{\mathbf A(t_0 - t)} \mathbf x(t),

исключим фазовую переменную из правой части уравнения:

\dot{\mathbf z}(t) = e^{\mathbf A(t_0 - t)} \mathbf f(t), \quad \mathbf z(t_0) = \mathbf x_0,

откуда легко получить окончательную формулу для решения исходной неоднородной системы:

\mathbf x(t) = e^{\mathbf A(t - t_0)} \left( \mathbf x_0 + \int_{t_0}^t {e^{\mathbf A(t_0 - \tau)} \mathbf f(\tau)} \,d\tau\right).

3. Пусть \mathbf f(t) \equiv \mathbf u = const. Пусть также \det\mathbf A \neq 0. Тогда имеем первообразную:

(e^{\mathbf A t} \mathbf A^{-1})' = e^{\mathbf A t}.

После подстановки первообразной в интеграл формула решения из предыдущего пункта обретает вид:

\mathbf x(t) = e^{\mathbf A(t - t_0)} (\mathbf x_0 + \mathbf A^{-1} \mathbf u) - \mathbf A^{-1} \mathbf u.

4. Предположим теперь, что матрица \mathbf A является функцией времени. Оставим также функцию \mathbf f. Можно воспользоваться тем же приёмом, что и в п.2, введя замену:

\mathbf z(t) = e^{\int_t^{t_0} {\mathbf A(\tau)}\,d\tau} \mathbf x(t).

Выполнив подстановку замены в уравнение, обнаружим, что для того, чтобы удалось убрать фазовую переменную из правой части, необходимо выполнение условия:

\forall t_1, t_2 \in [t_0, t] \left(\,\mathbf A(t_1) \mathbf A(t_2) = \mathbf A(t_2) \mathbf A(t_1)\,\right).

Тогда можно записать решение системы в форме:

\mathbf x(t) = e^{\int_{t_0}^t {\mathbf A(\tau)}\,d\tau} \left( \mathbf x_0 + \int_{t_0}^t {e^{\int^{t_0}_\tau {\mathbf A(s)}\,ds} \mathbf f(\tau)} d\tau\right).

5. Формулы для матричной экспоненты размерности 2. Пусть дана матрица

\mathbf A = \begin{bmatrix} a & b \\ c & d  \end{bmatrix}

Положим D = 4 bc + (a - d)^2. Имеем теперь три случая.

Случай D > 0.
Положим q = \sqrt{D}, C = \cosh\frac{q}{2}, S = \sinh\frac{q}{2}. Тогда

e^{\mathbf A} = \frac{1}{q} e^\frac{a+d}{2}  \begin{bmatrix} qC + (a - d) S & 2bS \\                             2cS & qC - (a - d)S  \end{bmatrix}

Случай D < 0.
Формула та же, что и в случае D > 0, если положить q = \sqrt{-D}, C = \cos\frac{q}{2}, S = \sin\frac{q}{2}.

Случай D = 0.

e^{\mathbf A} = e^\frac{a+d}{2}  \begin{bmatrix} 1 + \frac{a - d}{2} & b \\                             c & 1 - \frac{a - d}{2}  \end{bmatrix}

Для множества частных случаев можно вывести более простые формулы. Например, для диагональной матрицы b = c = 0:

e^{\mathbf A} = \begin{bmatrix}  e^a & 0 \\ 0 & e^d \end{bmatrix}

Для верхнетреугольной матрицы c = 0:

e^{\mathbf A} = \begin{bmatrix}  e^a & b\frac{e^a - e^d}{a - d}   \\ 0 & e^d \end{bmatrix}

Один алгоритм раскраски графа II

Недавно я переписал алгоритм субоптимальной правильной вершинной раскраски графа на Python. Решил выложить его сюда.

Вершина графа отождествляется с её индексом. Граф задаётся индексируемой последовательностью (tuple или list) наборов соседей (массив списков соседей). Раскраска — массив (list в терминах Python) цветов.

Итак, функция раскраски у меня имеет следующий вид:

## Построить правильную вершинную раскраску графа.
## \param graph -- граф в виде массива множеств номеров смежных вершин
## \param forbidden -- массив множеств запрещённых цветов для каждой вершины
## \param color -- массив цветов вершин
## \return массив (новых) цветов вершин
def coloring(graph, forbidden=None, color=None):
    if forbidden is None: forbidden = []
    if color is None: color = []
    # Предусловия.
    assert len(forbidden) == 0 or len(forbidden) == len(graph)
    assert len(color) == 0 or len(color) == len(graph)

    if len(forbidden) == 0:
        forbidden = [frozenset()]*len(graph)
    #...

Передав раскраску color, можно задать желаемые цвета, а алгоритм разрешит конфликты, если они есть. Если color уже содержит правильную раскраску, то она и будет возвращена функцией coloring. Массив forbidden позволяет запретить раскрашивать какие-то вершины в какие-то цвета: forbidden[v] есть множество запрещённых цветов (числовых меток) для вершины с индексом v.

Всё, что описано ниже, определено внутри функции coloring.

Если пользователь не передал некую начальную раскраску color, то построим её в два цвета поиском в глубину.

    # Раскраска в два цвета поиском в глубину.
    def coloring2():
        color2 = [0]*len(graph)
        def dfsColoring2(u, c=1):
            if color2[u] == 0:
                for d in range(c, len(graph)+1):
                    if d not in forbidden[u]:
                        color2[u] = d
                        break
                c = 3 - c
                for v in graph[u]:
                    dfsColoring2(v, c)
    
        for start in range(0, len(graph)):
            dfsColoring2(start)
        return color2

    # Построить начальную раскраску, если требуется.
    if len(color) == 0:
        color = coloring2()

Теперь собственно алгоритм. Он заключается в повторении пары действий: выявлении конфликтов цветов (listConflicts) и разрешении конфликтов цветов (resolveConflicts), пока все конфликты не будут разрешены.

    # Перечислить текущие конфликты цветов.
    def listConflicts():
        conflicts = []
        for u in range(0, len(graph)):
            nc = [color[v] for v in graph[u]]
            conflict_rank = nc.count(color[u])
            if conflict_rank > 0:
                nc = frozenset(nc) | forbidden[u]
                new_color = max(nc) + 1
                for c in range(1, new_color-1):
                    if c not in nc:
                        new_color = c
                        break
                conflict = (conflict_rank, new_color, u)
                conflicts.append(conflict)
        conflicts.sort(
            key=lambda c: (c[0], -c[1], c[2]), reverse=True)
        return conflicts
    
    # Разрешить конфликты цветов (потенциально не все).
    def resolveConflicts(conflicts):
        closed = set()
        for _, new_color, u in conflicts:
            if u in closed:
                continue
            color[u] = new_color
            closed |= set(graph[u])
    
    # Собственно алгоритм: начиная с некоторой раскраски,
    # повторять обнаружение и разрешение конфликтов, 
    # пока список конфликтов не опустеет.
    while True:
        conflicts = listConflicts()
        if len(conflicts) == 0:
            return color
        resolveConflicts(conflicts)

Цвет 0 не используется.

Пример:

g = ((1,2), (0,2), (0,1,3,6), (2,4,5), (3,5,6), (3,4), (2,4))
c = coloring(g)
print(c)
> [1, 2, 3, 2, 1, 3, 2]

Раскрашенный граф:

Вычисление массы молекулы по формуле вещества с помощью Python

Я привык писать на C++, и та задача, решение которой я собираюсь здесь продемонстрировать, сначала была мной решена на C++. Однако Python позволяет решить её с помощью намного более короткой программы (думаю, ряд других популярных языков в этом плане не хуже, но речь сейчас не о них). Это первая запись в моём блоге с использованием Python. Шапочное знакомство с этим языком у меня состоялось много лет назад, но пользоваться им не доводилось. Недавно я решил использовать его на занятиях, в первую очередь, из-за пакета модулей SciPy. В качестве первой задачки для себя я взял задачу, для которой достаточно стандартной функциональности.

Задача. Дана строка, содержащая формулу химического вещества. Требуется найти массу молекулы этого вещества в атомных единицах массы. Формула может содержать названия элементов вида K (одна заглавная латинская буква) или Mg (заглавная + строчная буква), круглые скобки (химики также используют и квадратные скобки, но упростим задачу), натуральные числа (количества атомов и групп атомов). Пример: «Ca (OH)2» (пробелы должны игнорироваться), масса равна сумме масс Ca и удвоенной суммы масс O и H.

Исходные данные: словарь (или файл) с массами элементов.

Первая идея была парсить формулу «руками» подряд, пропуская пробельные символы, выявляя названия элементов и числа, а скобочные формы обрабатывая рекурсивно. Подобная программа на C++ (с комментариями, циклом чтения в main и функцией чтения файла с массами элементов) вылилась примерно в 200 строк кода.

Вторая идея была использовать регулярное выражение для выявления частей формулы. Стандартная библиотека C++ включает поддержку регулярных выражений.

/// Break-up formula into tokens.
vector<string> tokenize_formula(const string &formula)
{
  // Regex object used for tokenizing formulae.
  static const regex tokens("[0-9]+|[A-Z][a-z]?|\\(|\\)");
  vector<string> result;
  sregex_iterator 
    from(begin(formula), end(formula), tokens), to;
  for_each(from, to, 
    [&](auto &match) { result.push_back(match.str()); });
  return result;
}

Дальше всё примерно так же, как при выделении лексем «руками». Такой вариант стал короче почти вдвое, но обзавёлся и важным недостатком — ошибки игнорировались, например, формула «N aHCO3» проходила, как если бы было написано «NHCO3», что неприятно.

Python позволяет записать то же самое ещё короче, но в нём есть возможность ещё сверх того ужать программу, не прибегая при этом к сторонним библиотекам. Причём, если в формуле есть ошибки, то мы это увидим. Впрочем, для начала я приведу вариант на основе предыдущего регэкспа, но уже с отловом ошибок (однако без чтения файла и main). Предполагается, что массы элементов заданы в словаре elements, скажем, так:

elements = {
    'H' :  1.008,
    'C' : 12.01,
    'N' : 14.007,
    'O' : 16,
    'Na': 22.99,
    'Mg': 24.306,
    'P' : 30.974,
    'S' : 32.068,
    'Cl': 35.45,
    'K' : 39.098,
    'Ca': 40.078
    }

Итак, код парсера, используещего методы findall и split:

import re
chemass_tokens = re.compile(r'[0-9]+|[A-Z][a-z]?|\(|\)|\[|\]')

def chemass(formula):
    """
    Формула -- это строка, составленная из знаков элементов, 
    количеств и скобок вроде 'H2O' или 'Ca (OH)2'. 
    Пробелы и всё прочее игнорируется (о выброшенном выводится сообщение).
    @return: масса
    """
    # Получить список токенов:
    tokens = chemass_tokens.findall(formula)
    m, l = do_chemass(tokens)
    # То, что не попало в токены, можно получить,
    # разбив формулу на части, разделённые токенами:
    unparsed = chemass_tokens.split(formula)
    # split возвращает кучу пустых строк, а также пробелы -- выкинем их.
    unparsed = [s for s in unparsed if s != '' and not str.isspace(s)]
    if len(unparsed) > 0:
        print('Chemass ignored the following formula parts:')
        print(unparsed)
    
    return m

def do_chemass(l):
    """
    Отдельная функция нужна для рекурсивной обработки вложенных скобочных форм.
    l -- список строк, где каждая строка -- 
    либо имя элемента, либо число, либо скобка.
    """
    mass = 0
    cm = 0
    while l:
        if str.isdecimal(l[0]):
            cm *= int(l[0])
        elif l[0] == '(' or l[0] == '[':
            mass += cm
            cm, l = do_chemass(l[1:])
            continue # чтобы не удалить голову списка l.
        elif l[0] == ')' or l[0] == ']':
            del l[0]
            break
        elif l[0] not in elements: # сообщить о неизвестном элементе.
            print('Chemass: '+l[0]+' is unknown element, ignored.')
        else:
            mass += cm
            cm = elements[l[0]]
        
        del l[0] # удалить голову списка l.

    mass += cm
    return mass, l

Ну, а теперь, «коронный номер». В Python есть встроенная функция eval, вычисляющая значение выражения (программы) на Python… Преобразуем-ка формулу в выражение, в духе ‘Ca(OH)2’ → ‘Ca+(O+H)*2’ и вычислим его! Для этого можно использовать метод регэкспа sub (подстановка по образцу). Далее привожу полученный код:

import re
chemass_names = re.compile('([A-Z][a-z]?)')
chemass_open = re.compile('(\()')
chemass_num = re.compile('([0-9]+)')

def chemass(formula):
    """
    Формула -- это строка, составленная из знаков элементов, 
    количеств и скобок вроде 'H2O' или 'Ca (OH)2'.
    @return: масса
    """
    return eval(chemass_formula(formula))

def chemass_formula(formula):
    """
    Функция, преобразующая строковую формулу 
    к выражению на Python, вычисляющему массу.
    """
    return chemass_num.sub(r'*\1',
        chemass_names.sub(r'+elements["\1"]',
        chemass_open.sub(r'+\1', formula)))

Это всего 20 строк. Если в формуле есть ошибки, о них сообщит сам парсер Python’а. Да, eval нельзя использовать в некоторых случаях, например, из соображений безопасности, но если это скрипт для личного пользования «по-быстрому», то это может быть очень удобно.

Addendum. Я использовал предварительную поддержку C++17 (g++ 7.2), и меня неприятно удивило отсутствие многих очевидных функций для работы со string_view. Например, его нельзя передать в конструктор fstream или в функцию stoi. Да, всегда можно явно сконвертировать в объект string, но это противоречит цели использования string_view: избежать создания лишних строк.

Пример функциональной декомпозиции

Задача вроде тех, что решают студенты на нашем экзамене после первого семестра.

«Дана прямоугольная матрица. Найти в ней первую строку, не содержащую нулей, или определить, что таких строк нет.»

При всей простоте этой задачи она предполагает широкий спектр возможных решений. К сожалению, не смотря на то, что на практике и в самостоятельных работах подобные задачи есть, у многих студентов они всё равно вызывают затруднения на экзамене.

Решение в стиле C подразумевает прототип функции вроде (матрицу можно оформить по-разному, но предпочтём пока простейший вариант):

float* find_row_without_zero
  (float **a, size_t rows, size_t cols);

— вернуть указатель на соответствующую строку, либо nullptr, если такой нет.

Итак, студент видит слово «матрица», из чего делает вывод, что здесь явно двойной цикл, и пишет что-то такое:

float* find_row_without_zero
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
  {
    for (size_t j = 0; j < cols; ++j)
    {
      if (a[i][j] == 0)
        // ???
      else
        // ???

Почему-то многим очень хочется поставить else-секцию с возвратом результата. Т.е. ошибочно подменить квантор «все» (не нули) на «существует». Под if в этом случае идёт break. Положим, стало ясно, что else здесь никак не нужен. Что же делать? Типичное неказистое решение имеет следующий вид:

// Bad structured programming (using a flag variable).
float* find_row_without_zero_1
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
  {
    bool found = true;
    for (size_t j = 0; j < cols; ++j)
    {
      if (a[i][j] == 0)
      {
        found = false;
        break;
      }
    }

    if (found)
      return a[i];
  }

  return nullptr;
}

В принципе, кто-то может решиться на использование goto, что даёт избавление от лишней переменной. Но это маловероятно: в примерах goto почти не встречается, а в книгах последних десятилетий его использование традиционно порицается.

// I'm not afraid of GOTO!
float* find_row_without_zero_2
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
  {
    for (size_t j = 0; j < cols; ++j)
      if (a[i][j] == 0)
        goto next_row;

    return a[i];
  next_row:;
  }

  return nullptr;
}

На самом деле, данная задача подразумевает применение функциональной декомпозиции (основного метода программирования, кстати). Достаточно понять, что решение по сути является композицией двух линейных поисков.

// Good structured programming:
// applying functional decomposition.
float* find_value
  (float *a, size_t size, float val)
{
  for (size_t i = 0; i < size; ++i)
    if (a[i] == val)
      return a + i;
  return nullptr;
}

float* find_row_without_zero_3
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
    if (!find_value(a[i], cols, 0))
      return a[i];
  return nullptr;
}

Если предполагается использование STL (после первого семестра это ещё не предполагается), то явный цикл можно заменить на вызов std::find, так:

// Ah, std::find does the work in C++.
float* find_row_without_zero_4
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
    if (std::find(a[i], a[i] + cols, 0.f) == a[i] + cols)
      return a[i];
  return nullptr;
}

или даже так (ещё std::find_if, шаблон по типу данных и значение по умолчанию вместо явного нуля):

// It's not beautiful, I guess...
template <class NT>
NT* find_row_without_zero_5
  (NT **a, size_t rows, size_t cols)
{
  const auto p = std::find_if(a, a + rows,
    [cols](NT *a)
    {
      return std::find(a, a + cols, NT{}) == a + cols;
    });
  return p == a + rows? nullptr: *p;
}

Последний вариант, впрочем, наверно, избыточно «навороченный».

Ну и наконец, можно поразвлекаться в функциональном стиле (требуется C++14), отказавшись от указателей и рассматривая матрицу как диапазон диапазонов. Будем возвращать итератор в стиле STL.

using std::forward;

template <class Item>
auto equals(Item &&item)
{
  return [it = forward<Item>(item)]
  (auto &&other)
  {
    return other == it;
  };
}

template <class Pred>
auto find(Pred &&pred)
{
  return [pr = forward<Pred>(pred)]
  (auto &&rng)
  {
    using std::begin;
    using std::end;

    const auto e = end(rng);
    auto p = begin(rng);
    while (p != e && !pr(*p))
      ++p;

    return p;
  };
}

template <class Pred>
auto none(Pred &&pred)
{
  return [pr = forward<Pred>(pred)]
  (auto &&rng)
  {
    using std::end;
    return find(pr)(rng) == end(rng);
  };
}

Вышеприведённые определения позволяют в итоге записать решение в виде композиции:

template <class Rng2>
auto find_row_without_zero_6(Rng2 &&rows)
{
  return find(none(equals(0)))(rows);
}

Забавно, что такая конструкция применима не только, например, к вектору векторов, но и к статическому двумерному массиву:

float m[3][4]
{
  { 1, 1, 1, 0 },
  { 0, 1, 1, 1 },
  { 1, 0, 1, 1 },
};

assert(find_row_without_zero_6(m) == std::end(m));

Обобщённая функция удаления элементов из контейнера

Введение

Стандартная библиотека C++ не предлагает общего интерфейса для удаления элементов из произвольного контейнера. Тем не менее, средства метапрограммирования и принцип SFINAE позволяют реализовать такой интерфейс (например, в виде функции-шаблона) поверх стандартных компонент и выполнять выбор подходящего для переданного контейнера метода на этапе компиляции. Наша функция будет удалять элементы контейнера, удовлетворяющие произвольному предикату, переданному в виде функтора. Обычно это линейная по числу элементов контейнера операция, в случае же ассоциативных контейнеров, реализованных в виде деревьев, она может быть линеарифмической.

Далее предполагаем следующую «шапку».

#include <algorithm>
#include <iterator>
#include <type_traits>
#include <utility>
using namespace std;

Реализации свободной функции для разных видов контейнеров

Контейнеры-списки реализуют желаемую функцию в виде функции-члена, которой и следует переадресовать вызов.

/// Удаляет элементы, для которых pred возвращает истину,
/// из связных списков (std::list, std::forward_list).
template <class List, class Pred> inline
void list_remove_if(List &l, Pred pred)
{
  l.remove_if(pred);
}

Прочие линейные контейнеры с изменяемым размером (вектор, дек) предполагают перемещение элементов в начало (затирая удаляемые) с сохранением относительного порядка и «отрезание хвоста» контейнера.

/// Удаляет элементы, для которых pred возвращает истину,
/// из линейных контейнеров (std::vector, std::deque).
template <class Vector, class Pred> inline
void vector_remove_if(Vector &v, Pred pred)
{
  v.erase(remove_if(begin(v), end(v), pred), end(v));
}

Ассоциативные контейнеры требуют прямого обхода с удалением конкретных элементов.

/// Удаляет элементы, для которых pred возвращает истину,
/// из ассоциативных контейнеров (разные варианты set и map).
/// В случае map предикат вычисляется для пары ключ-значение.
template <class Assoc, class Pred>
void assoc_remove_if(Assoc &a, Pred pred)
{
  for (auto it = begin(a), it_end = end(a); it != it_end; )
    if (pred(*it))
      it = a.erase(it);
    else
      ++it;
}

 
Автоматический выбор правильного варианта по виду контейнера

Автоматический выбор между представленными тремя вариантами можно организовать через банальную лобовую перегрузку функции remove_if для std::list, std::vector, std::set и т.д. Но мне не нравится этот вариант, так как он не способен автоматически «принять» нестандартные контейнеры, удовлетворяющие интерфейсу того или иного STL-контейнера — для каждого вида контейнера надо будет написать явную перегрузку.

Вместо этого мы можем заставить компилятор автоматически выбирать тот или иной вариант в зависимости от того, какие операции доступны для контейнера (сделать это позволяет SFINAE).

Вариант list_remove_if годится тогда, когда контейнер предоставляет метод remove_if (сам контейнер лучше знает как эффективно удалять из него элементы). Вариант vector_remove_if годится для контейнеров с методом удаления диапазона erase(iterator, iterator) и итератором произвольного доступа (формально достаточно итератора однонаправленного последовательного доступа («forward»), но тогда будут подходить и ассоциативные контейнеры). Наконец, вариант assoc_remove_if годится для всех контейнеров с методом iterator erase(iterator).

Далее приводится код метафункций, проверяющих наличие тех или иных операций и таким образом определяющих вид контейнера.

namespace impl
{
  /// Категория итератора для условного типа "не-итератора".
  struct Not_an_iterator_tag {};
  /// Специальный тип, означающий отсутствие типа,
  /// но позволяющий создавать объекты (в отличие от void).
  struct No_type_tag {};
  /// Условный тип "не-итератор", имитирует итератор с
  /// определёнными вложенными типами std::iterator_traits.
  struct False_iterator
    : iterator<Not_an_iterator_tag, No_type_tag> {};
 
  /// Метафункция. Извлекает тип итератора коллекции X 
  /// (использует для этого [std::]begin).
  /// Если begin(X) не определена, возвращает False_iterator.
  template <class X>
  class Collection_iterator
  {
    template <class Y>
    static False_iterator check(unsigned long, Y&&);
 
    template <class Y>
    static auto check(int, Y &&y) -> decltype(begin(y));
 
  public:
    using type = decltype(check(1, declval<X>()));
  };
 
  /// Проверяет, является ли X "коллекцией".
  /// Коллекция -- объект, для которого свободные функции
  /// begin и end возвращают итераторы с категорией,
  /// приводимой к forward_iterator_tag.
  /// Проверяется только [std::]begin.
  /// В качестве коллекции может выступать любой STL контейнер,
  /// объект std::basic_string или статический массив.
  template <class X>
  struct Is_collection
    : is_convertible<typename
        iterator_traits<typename Collection_iterator<X>::type>::iterator_category,
        forward_iterator_tag> {};
 
  ///////////////////////////////////////////////////////////////////////////// 
  /// Проверка, возможен ли вызов assoc_remove_if для объекта Cont.
  template <class Cont>
  class Assoc_remove_if_possible
  {
    template <class C>
    static false_type check(unsigned long, C&&);
 
    template <class C>
    static auto check(int, C &&c) ->
      is_convertible<
        decltype(c.erase(begin(c))),
        decltype(begin(c))>;
 
    static const bool
      erase_defined = decltype(check(1, declval<Cont>()))::value,
      is_collection = Is_collection<Cont>::value;
 
  public:
    using type = integral_constant<bool, erase_defined && is_collection>;
    static const bool value = type::value;
  };
 
  /// Проверка, возможен ли вызов vector_remove_if для объекта Cont.
  template <class Cont>
  class Vector_remove_if_possible
  {
    template <class C>
    static false_type check(unsigned long, C&&);
 
    template <class C>
    static auto check(int, C &&c) -> integral_constant<bool,
         is_convertible<typename
            iterator_traits<typename Collection_iterator<C>::type>::iterator_category,
            random_access_iterator_tag>::value
      && is_convertible<
            decltype(c.erase(begin(c), end(c))),
            decltype(begin(c))>::value>;
 
  public:
    using type = decltype(check(1, declval<Cont>()));
    static const bool value = type::value;
  };
 
  /// Проверка, возможен ли вызов list_remove_if.
  template <class Cont>
  class List_remove_if_possible
  {
    template <class C>
    static false_type check(unsigned long, C&&);
 
    // Универсальный предикат.
    struct Unipred
    {
      template <class T> 
      bool operator()(const T&) { return false; }
    };
 
    template <class C>
    static auto check(int, C &&c) ->
      decltype((void)c.remove_if(Unipred{}), true_type{});
 
  public:
    using type = decltype(check(1, declval<Cont>()));
    static const bool value = type::value;
  };
 
  ///////////////////////////////////////////////////////////////////////////// 
  // Статическая диспетчеризация вызова
  // Теги для статического выбора соответствующего варианта remove_if.
  struct Assoc_remove_if_tag {};
  struct Vector_remove_if_tag {};
  struct List_remove_if_tag {};
 
  /// Выбор подходящего тега remove_if для контейнера.
  template <class Cont>
  class Remove_if_tag
  {
    static const bool
      assoc_possible  = Assoc_remove_if_possible<Cont>::value,
      vector_possible = Vector_remove_if_possible<Cont>::value,
      list_possible   = List_remove_if_possible<Cont>::value;
 
    static_assert(assoc_possible || vector_possible || list_possible,
        "remove_if can not be called for this type");
 
  public:
    using type = conditional_t<list_possible, List_remove_if_tag,
                 conditional_t<vector_possible, Vector_remove_if_tag,
                 conditional_t<assoc_possible, Assoc_remove_if_tag,
                               void>>>;
  };

  // Перенаправление вызова remove_if на варианты для различных типов контейнеров.
  template <class Cont, class Pred> inline
  void remove_if(Cont &c, Pred pred, Assoc_remove_if_tag)
  { assoc_remove_if(c, pred); }
 
  template <class Cont, class Pred> inline
  void remove_if(Cont &c, Pred pred, Vector_remove_if_tag)
  { vector_remove_if(c, pred); }
 
  template <class Cont, class Pred> inline
  void remove_if(Cont &c, Pred pred, List_remove_if_tag)
  { list_remove_if(c, pred); }
}

/// Вызов разрешается статически в подходящий вариант в зависимости от типа контейнера.
template <class Cont, class Pred> inline
void remove_if(Cont &container, Pred pred)
{
  impl::remove_if(container, pred, typename impl::Remove_if_tag<Cont>::type{});
}

Итак, remove_if(контейнер, предикат) выбирает один из трёх способов удаления элементов, удовлетворяющих предикату, в зависимости от операций, предоставляемых контейнером. Приведённый код целиком и небольшой проверочный код: ideone, Bitbucket.

UPD. Увидел, что нечто подобное предлагалось для добавления в стандартную библиотеку. Пока оно представлено в виде расширения в пространстве имён experimental.

О реализации одного алгоритма топологической сортировки

Как-то так получилось, что именно записи на в общем-то побочную для меня тему алгоритмов на графах привлекают основной объём посещений блога «со стороны». Отчасти поэтому, а также потому, что на лекции я предложил было не слишком разумный вариант реализации одного из алгоритмов топологической сортировки вершин ориентированного ациклического графа (DAG — «directed acyclic graph»), я решил написать данный текст и соответствующий код.

На изображении ниже представлен DAG, использованный для проверки реализованной программы.

Пример DAG

Пример DAG

Рассматриваемый алгоритм весьма прост. Его можно описать следующим образом.

Вход: граф G = (V, E).
Выход: последовательность R.
Пусть R = ∅.
Пока |R| < |V|:
    пусть v ∈ V \ R такая что indeg(v) = 0;
    дописать v к R;
    удалить v из G.

Хотя этот алгоритм описан как «деструктивный» по отношению к графу G, на практике совсем не обязательно «разрушать» пусть даже копию некоторого представления G. Его можно переформулировать следующим образом (предполагая представление графа в виде списков соседей).

Создать A = [ indeg(v) | v ∈ V ] — массив степений захода вершин.
Завести очередь Q и поместить туда все v ∈ V: A[v] = 0;
Пока Q ≠ ∅:
    извлечь очередную v из Q;
    поместить v в R;
    для всех w ∈ out(v):
        если A[w] = 1: поместить w в Q;
        если A[w] > 0: A[w] := A[w] – 1.

Как нетрудно видеть, данный алгоритм имеет сложность O(|V| + |E|), при предположении, что затраты на перечисление соседей вершины пропорциональны их количеству. Никакие «сложные» структуры данных не используются.

Ссылка на реализацию на C++.