Author Archive: evetro

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

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}

О методе главных компонент

На днях пришлось столкнуться с «методом главных компонент» (principal component analysis, PCA). Ниже даю обретённое мною понимание оного.

Дано: набор точек { pi ∈ ℝd | i∈1:n }, dn.

Предположение: точки получены аффинным преобразованием реализации случайного процесса y: pi = μ + yiW, где μ ∈ ℝd, W — некая матрица, а yi удовлетворяют многомерному нормальному распределению с нулевым матожиданием и единичной матрицей ковариации.

Задача: найти μ, W и yi.

Замечание 1. Размерность yi не обязана совпадать с размерностью pi.

Замечание 2. В качестве оценки матожидания μ возьмём арифметическое среднее pi: μ = Σi pi / n. Обозначим: xi = pi − μ.

Замечание 3. Так как покоординатно независимое нормальное распределение порождает «шарообразное» скопление точек yi, то следует ожидать, что xi формируют скопление, визуально напоминающее повёрнутый d-осный эллипсоид. Иными словами, можно искать матрицу W как комбинацию анизотропного масштабирования и поворота: W = TR, где T может иметь ненулевые элементы только на главной диагонали, а R ортогональна.

Через X обозначим n×d-матрицу, составленную из строк xi.

«Поточечная» матрица ковариации XX имеет размеры n×n. «Покоординатная» матрица ковариации XX имеет размеры d×d. Обе эти матрицы симметричны, но в общем случае не диагональны.

Через Y обозначим n×r-матрицу, составленную из строк yi. По предположению имеем YY = 1n или YY = 1r. Выберем пока первый вариант. Из него немедленно следует rn. Логично положить r = n, тогда Y ортогональна.

Итак, X = YTR, где Y и R ортогональны, а T имеет ненулевые элементы только на главной диагонали. Данное разложение матрицы X есть ни что иное, как разложение сингулярных значений (singular value decomposition, SVD).

Возьмём любую реализацию SVD и построим разложение: X = USV*. В общем случае матрицы U и V унитарны, звёздочкой обозначено эрмитово сопряжение. В нашем случае они ортогональны и сопряжение есть транспонирование, поэтому можно положить:

Y = U, T = S, R = V, W = SV.

Замечание 4. На диагонали S стоят сингулярные значения, в нашей задаче имеющие геометрический смысл длин полуосей эллипсоида. Они неотрицательны и по построению отсортированы по величине от больших к меньшим. В конце могут стоять нули, означающие, что по соответствующим направлениям (которым отвечают столбцы V) эллипсоид вырожден, и эти измерения можно отбросить. С практической точки зрения может быть разумным отбросить все измерения, для которых полуось (сингулярное значение) меньше некоторой заданной величины. WW = SS — диагональная матрица, составленная из квадратов сингулярных значений.

Шаг 1. A = XV. Так как V есть матрица, обратная к R, то A = YT. Это эллипсоид, повёрнутый так, чтобы его оси совпали с осями координат стандартного базиса, причём разброс данных на каждой следующей оси не больше, чем на предыдущей. Кому-то этого результата уже может быть достаточно.

Шаг 2. Пусть задано ρ ≥ 0 — предельная «чувствительность» по измерению. Исходная полная матрица S может быть неквадратной и содержать нулевой блок под квадратом d×d, здесь я его проигнорирую и положу S = diag(sj), 1 ≤ jd.

Пусть r = max{ j | sj > ρ }. Положим B = A[1:r], т.е. первые r столбцов матрицы A. Строки bi матрицы B есть точки в пространстве уменьшенной размерности r ≤ d.

Шаг 3. Составим усечённую псевдообратную матрицу S+ = diag(1/sj), 1 ≤ jr.
Положим Z = BS+. Так как S+ есть матрица, обратная к усечённой T, то мы получим, что Z = Y[1:r], т.е. первые r столбцов Y. Если мы уже вычислили U, то можно просто взять из неё первые r столбцов: Z = U[1:r]. Имеем ZZ = 1r.

Точки zi имеют уменьшенную размерность r ≤ d и отмасштабированы по осям координат.

Гладкая интерполяция между синусоидой и прямоугольной волной

Представим, что есть задача получить нечто среднее между синусоидой

\mathrm{Sin}\,(t) \triangleq \sin{2\pi t}

и прямоугольной волной

\mathrm{Rect}\,(t) \triangleq \mathrm{sgn}\,\mathrm{Sin}\,(t).

Очевидно, что банальная поточечная линейная интерполяция с параметром \alpha \in [0, 1]:

(1-\alpha)\,\mathrm{Sin}\,t + \alpha\,\mathrm{Rect}\,t

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

R_\alpha\,(t) \triangleq \dfrac{\mathrm{arctg}\,(\alpha\,\mathrm{Sin}\,(t))}{\mathrm{arctg}\,\alpha}.

Здесь уже \alpha \in (0, +\infty).
Легко видеть, что

\lim\limits_{\alpha\to 0} R_\alpha\,(t) = \mathrm{Sin}\,(t),\quad \lim\limits_{\alpha\to +\infty} R_\alpha\,(t) = \mathrm{Rect}\,(t).

Графики R_\alpha\,(t), t \in [0, 2], \alpha \in \{\,1, 10, 1000\,\}:

Один алгоритм раскраски графа 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]

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

«Мы» или «вы»

В данном посте я продолжаю заниматься лингвистическим дилетантизмом. Впрочем, здесь я, в общем-то, не выдвигаю новых гипотез, а просто привожу подборку материала для давно выдвинутой гипотезы (по крайней мере, Т.В.Гамкрелидзе и В.В.Ивановым, см. ЭССЯ *my, т. 21, с. 22, см. также Википедию), касающейся происхождения индоевропейских личных местоимений. Эта гипотеза мне кажется интересной.

Легко заметить, что слово «мы» может иметь различный смысл:

  • «я/мы и вы»;
  • «мы, но не вы».

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

Вначале приведу две таблицы. Формы 2л. ед.ч. приведены для полноты.

Формы глагола «есть» (дефис отделяет корень, реконструкция взята с en.wiktionary.org):

форма ст.-слав. ПИЕ
1л, ед.ч. (я) ес-мь *h₁és-mi
1л. мн.ч. (мы) ес-мъ *h₁s-mós
2л. ед.ч. (ты) ес-и *h₁és-(s)i
2л. мн.ч. (вы) ес-те *h₁s-té

Местоимения «я», «мы», «вы» в разных ИЕ-языках + реконструкции для праиндоевропейского (ПИЕ) и прауральского (ПУ) (по данным en.wiktionary.org, по поводу сходства ПИЕ и ПУ см. также). Столбец «форма» (ф.) расшифровывается как лицо-число-падеж:

ф. лит. лтш. лат. санскр. др.-англ. др.-арм. ПИЕ (Sihler) ПУ
1еи я àš es ego aham iċ, ih es *eǵoH *mun, *minä
1ер меня manę̃s manis meī mama mīn im *mé-me, *m-os ?
1ев меня manè mani mām meċ, mē is *mé ?
1ми мы mẽs mēs nōs vayam mekʿ *we-i *me
1мр нас mū́sų mūsu nostrī asmākam, naḥ ūser, ūre mer *n̥s-óm, *no-s ?
1мв нас mùs mūs nōs asmān, naḥ ūsiċ, ūs mez *n̥s-mé, *no-s ?
2еи ты tu tvám þū du *ti/*tu *tun, *tinä
2ер тебя tavę̃s tevis tuī táva, te þīn kʿo *té-we, *t-os ?
2ев тебя tavè tevi tvā́m, tvā þeċ, þē kʿez *t-wé, *te ?
2ми вы jū̃s jūs vōs yūyám ġē dukʿ *yūs *te
2мр вас jū́sų jūsu vestrī yuṣmā́kam, vaḥ ēower jer *us-óm, *wo-s ?
2мв вас jùs jūs vōs yuṣmā́n, vaḥ ēowic, ēow jez *us-mé, *wo-s ?

ПИЕ реконструкция 1л.мн.ч. им. п. основана на данных санскрита, германских и анатолийских (нет в таблице) языков. Вообще, влияние санскрита на индоевропеистику было всеобъемлющим. Наверное, не будет преувеличением сказать, что именно ознакомление европейцев с санскритом положило начало современной индоевропеистике, и первые реконструкции праиндоевропейского языка были очень близки санскриту. Вероятно, несколько избыточная сила этого влияния сохраняется до сих пор. В данном аспекте традиционная реконструкция подкрепляется ещё и древнейшим известным письменным ИЕ языком — хеттским. Впрочем, общая «инаковость» анатолийских языков приводит к тому, что их порой отделяют от основного набора ИЕ языков в близкородственную семью.

Очевидный параллелизм форм в балтийских, армянском, греческом (нет в таблице) и латыни позволяет толковать его как вторичную потерю особых форм им.п. по аналогии с формами косвенных падежей (что, очевидно, произошло в латыни), однако есть маленькая «неприятность» в виде славянских языков: «мы» vs «нас». Авторы ЭССЯ объясняют форму слав. *my как результат влияния форм *vy и *ny (2л.мн.ч. вин. п.) на исходное ПИЕ *mes, сохранившееся в балтийских. Эта форма толкуется как мн. ч. исходного *me («я»), что хорошо сочетается с глагольными окончаниями.

Итак, гипотеза состоит в следующем: *mes представляет собой эксклюзивное «мы», в то время как *wei или *wes представляет собой инклюзивное «мы», происходя от мн.ч. местоимения 2л. Можно развить идею родства и форм косвенных падежей: редукция *yūs > *us-, *mes > *m̥s- > *n̥s-, аналогичный последнему переход может наблюдаться в окончании мн.ч. вин.п. существительных.

Ввиду отсутствия известных ИЕ языков с противопоставлением инклюзивного и эксклюзивного «мы», а также вообще с одновременным наличием продолжений *mes и *wes/*wei, предположение о наличии такого противопоставления в древности находится на зыбком основании. Что, впрочем, никак не вредит гипотезе о первичности исходной праформы *mes, дополняемой различными объяснениями возникновения *wes/*wei (например, от дв.ч. «мы (оба)»).

Кубический сплайн перестройки частоты

Предположим, требуется сгенерировать «синусоиду» с плавной перестройкой частоты. Частота зависит от времени? Пожалуйста:

S(t) \triangleq \sin(2\pi f(t) t),

где f(t) — функция, управляющая частотой (в том, что это ещё не есть частота per se убедимся на примере).

Итак, пусть заданы два числа f_1 и f_2 — начальная и конечная частоты. Ну что же, ничтоже сумняшеся строим f(t) как кубический сплайн с условиями f(0) = f_1, f(1) = f_2, f'(0) = f'(1) = 0. Из этих условий легко получить, что

f(t) = f_1 + (f_2 - f_1)(3t^2 - 2t^3).

Результаты аналитических манипуляций надобно проверять практикой. Подставим в синус и построим перестройку с 110Гц на 330Гц и обратно (1с 110Гц, 1с перестройка, 1с 330Гц, 1с перестройка обратно, 1с 110Гц): звук 1. (Звук внедрён в pdf, поскольку бесплатный план WordPress.com не позволяет размещать звуковые файлы.)

Мда, что-то тут не так… Если построить график f(t), то он вроде бы выглядит идеально. Очевидно, считать эту функцию частотой (высотой звука в данном контексте) ошибочно. Я не хочу лезть глубже элементарного анализа (никаких преобразований Фурье и иже с ними), поэтому рассмотрю задачу с другой стороны. Если частота k постоянна, то имеем:

(\sin(2\pi k t))' = 2\pi k \cos(2\pi k t)).

Что если частоту определять через производную функции S(t)?

S'(t) = 2\pi (f'(t) t + f(t)) \cos(2\pi f(t) t)).

Определим новую функцию

F(t) \triangleq (f'(t) t + f(t)) = 4at^3 + 3bt^2 + 2ct + d.

Опять будем искать f(t) в форме кубического четырёхчлена f(t) \triangleq at^3 + bt^2 + ct + d, но условия изменим. Пусть теперь F(0) = f_1, F(1) = f_2, F'(0) = F'(1) = 0. Получим:

f(t) = f_1 + (f_2 - f_1)(t^2 - \frac{1}{2} t^3).

Подставим в синус и построим такую же перестройку с 110Гц на 330Гц и обратно: звук 2. Ну что же, теперь всё в порядке!

Конкатенация и транспонирование tuple

В недавно вышедшей Boost 1.66 beta 1 присутствует новая библиотека метапрограммирования Mp11. Я начал читать отсюда и обратил внимание на два примера.

Первый — конкатенация произвольного набора кортежей типов. Мне показалось очень сложным представленное решение Эрика Ниблера, поскольку «на глаз» было видно, что стандартными средствами это сделать довольно легко, ЕСЛИ ограничиться поддержкой std::tuple. Естественно, я попробовал написать свою версию, которую и привожу ниже.

Вначале конкатенация двух tuple:

template <class A, class B>
struct Tcat;

template <class... T, class... S>
struct Tcat<std::tuple<T...>, std::tuple<S...>>
{
  using type = std::tuple<T..., S...>;
};

template <class A, class B>
using tcat = typename Tcat<A, B>::type;

Теперь не составит труда определить конкатенацию произвольного количества tuple рекурсивно:

template <class... A> struct Tcats;

template <class... T>
struct Tcats<std::tuple<T...>>
{
  using type = std::tuple<T...>;
};

template <class A, class B, class... Other>
struct Tcats<A, B, Other...>
{
  using type = typename Tcats<tcat<A, B>, Other...>::type;
};

template <class... T>
using tcats = typename Tcats<T...>::type;

Второй — transform бинарной операцией двух кортежей (zip, если в терминах функциональных языков программирования). Конкретный случай для двух кортежей не столь интересен. Легко видеть, что при обобщении на произвольное число кортежей мы получаем ни что иное как транспозицию матрицы (F — применяемая операция):


<T11, T12, T13, ..., T1m>,
<T21, T22, T23, ..., T2m>,
...
<Tn1, Tn2, Tn3, ..., Tnm>
--->
F<T11, T21, ..., Tn1>,
F<T12, T22, ..., Tn2>,
F<T13, T23, ..., Tn3>,
...
F<T1m, T2m, ..., Tnm>

Насколько просто сделать такую штуку на шаблонах C++? Результат, который у меня получился, мне не особенно понравился. Думаю, можно это сделать лучше, однако, не факт, что это «лучше» будет означать «короче».

Итак, метафункция, порождающая из столбца исходной матрицы строку новой матрицы с префиксом F (n-арной операцией):

template <class I,
  template <class...> class F, class... T>
struct Tcol
{
  using type = F<std::tuple_element_t<I::value, T>...>;
};

template <class I,
  template <class...> class F, class... T>
using tcol = typename Tcol<I, F, T...>::type;

Здесь I — тип, определяющий зависимое имя value, раскрываемое в целочисленную константу — индекс столбца.

Теперь надо применить tcol к диапазону индексов столбцов. Так как готовых iota, map и bind у меня нет, я написал рекурсивную реализацию аналога «for» «в лоб»:

template <class S, class T,
  template <class...> class F, class... Other>
struct Tfor;

template <class Int, Int s,
  template <class...> class F, class... Other>
struct Tfor<
  std::integral_constant<Int, s>,
  std::integral_constant<Int, s>,
  F, Other...>
{
  using type = std::tuple<>;
};

template <class Int, Int s, Int t,
  template <class...> class F, class... Other>
struct Tfor<
  std::integral_constant<Int, s>,
  std::integral_constant<Int, t>,
  F, Other...>
{
  static_assert(s < t, "");
  using S  = std::integral_constant<Int, s>;
  using S1 = std::integral_constant<Int, s+1>;
  using T  = std::integral_constant<Int, t>;
  
  using type = tcat<std::tuple<F<S, Other...>>,
    typename Tfor<S1, T, F, Other...>::type>;
};

template <class S, class T,
  template <class...> class F, class... Other>
using tfor = typename Tfor<S, T, F, Other...>::type;

Здесь S и T — типы, задающие начало и конец диапазона (::value — соответствующие целочисленные константы), пак Other — исходная матрица типов.

Теперь записать обобщённый zip не составит труда:

template <
  template <class...> class F,
  class... T>
struct Tzip;

template <
  template <class...> class F>
struct Tzip<F> { using type = std::tuple<>; };

template <
  template <class...> class F,
  class T0, class... T>
struct Tzip<F, T0, T...>
{
  template <class I, class... Y>
  using tcolF = tcol<I, F, Y...>;
  using type = tfor<
    std::integral_constant<std::size_t, 
      0>,
    std::integral_constant<std::size_t,
      std::tuple_size<T0>::value>,
    tcolF, T0, T...>;
};

template <
  template <class...> class F,
  class... T>
using tzip = typename Tzip<F, T...>::type;

Пришлось вынуть первый кортеж (строку матрицы) — T0, чтобы узнать количество столбцов. Представляю, какое веселье будет в выводе компилятора, если скормить последовательность кортежей разных размеров :)

Тестировочный код:

template <class... T>
struct Tmax;

template <class T>
struct Tmax<T> { using type = T; };

template <class T0, class... T>
struct Tmax<T0, T...>
{
  using tail_max = typename Tmax<T...>::type;
  using type = std::conditional_t<
    (T0::value < tail_max::value), tail_max, T0>;
};

template <class... T>
using tmax = typename Tmax<T...>::type;

template <class T>
using Sizeof =
  std::integral_constant<std::size_t, sizeof(T)>;

template <class... T>
using tmax_sizeof = tmax<Sizeof<T>...>;

int main()
{
  using namespace std;
  // Concat.
  {
    using t0 = tuple<int, int*>;
    using t1 = tuple<float>;
    using t2 = tuple<char*, char**, void*, void**>;
    using t3 = tuple<>;
    using ex = tuple<
      int, int*, float, char*, char**, void*, void**>;
    static_assert(is_same<tcats<t0, t1, t2, t3>, ex>::value);
  }
  
  // Generic zip.
  {
    using r0 = tuple<char, short, int, long>;
    using r1 = tuple<unsigned, void*, void**, void***>;
    using exr = tuple<
      tuple<char, unsigned>,
      tuple<short, void*>,
      tuple<int, void**>,
      tuple<long, void***>>;
    //print_t(tzip<tuple, r0, r1>{});
    static_assert(is_same_v<exr, tzip<tuple, r0, r1>>);
    
    using t0 = tuple<int[11], int[12], int[13]>;
    using t1 = tuple<int[21], int[22], int[23]>;
    using t2 = tuple<int[31], int[32], int[33]>;
    using t3 = tuple<int[41], int[42], int[43]>;
    using ext = tuple<
      tuple<int[11], int[21], int[31], int[41]>,
      tuple<int[12], int[22], int[32], int[42]>,
      tuple<int[13], int[23], int[33], int[43]>>;
    static_assert(
      is_same_v<ext, tzip<tuple, t0, t1, t2, t3>>);
    
    using exsz = tuple<
      Sizeof<int[41]>, Sizeof<int[42]>, Sizeof<int[43]>>;
    static_assert(
      is_same_v<exsz, tzip<tmax_sizeof, t0, t1, t2, t3>>);
  }
}

Вычисление массы молекулы по формуле вещества с помощью 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: избежать создания лишних строк.

Проверка возможности записи через итератор вывода

Периодически возникает желание иметь статическую проверку согласованности типа значения итератора вывода и типа записываемого по итератору значения. К сожалению, в том случае, когда в качестве итератора передаётся настоящий итератор вывода, может оказаться, что iterator_traits<OutIt>::value_type есть void. Таковы, в частности, стандартные итераторы вывода, например, back_insert_iterator и ostream_iterator. Причина такого поведения заключается в том, что данные итераторы допускают только присваивание через разыменование и инкремент (который может ничего не делать):

*out = val;

Тип значения, возвращаемого оператором * обычно не является ссылкой на тип записываемого значения (это может быть прокси-объект, в качестве которого может выступать и сам итератор). Подстановка void в качестве value_type приводит к скорым ошибкам компиляции при попытке некорректно использовать этот тип.

Если требуется только проверить доступность такой операции записи для заданного типа значения, то это следует сделать напрямую — проверить возможность указанного выше присваивания:

#include <type_traits>

/// Check if a value type can be written via an iterator.
template <class Val, class It>
using Is_writable_to = typename std::is_assignable<
    decltype(*std::declval<It>()), Val>::type;
     
/// Check if a value type can be written via an iterator.
template <class Val, class It>
constexpr bool is_writable_to 
    = Is_writable_to<Val, It>::value;

Теперь следующий код выводит «110» (double можно привести к int, а void* — нельзя):

#include <iostream>
#include <iterator>
#include <vector>

int main()
{
  using namespace std;
  vector<int> a;
  auto bi = back_inserter(a);
   
  cout << is_writable_to<int, decltype(bi)>;
  cout << is_writable_to<double, decltype(bi)>;
  cout << is_writable_to<void*, decltype(bi)>;
  return 0;
}