Author Archive: evetro

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

На днях пришлось столкнуться с «методом главных компонент» (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;
}

Храните деньги в банке

Храните деньги в банке. А банку — под кроватью. В предыдущей записи я неосмотрительно сравнил слово банка (с вариантом баночка) со словом доска, заявив в нём корень банк-. Ну действительно, вроде в русском нет подходящего слова без суффикса (**бана?), да и само слово похоже на недавнее заимствование… Естественно, возникает вопрос, родственны ли слова банк (очевидное позднее заимствование, вероятно, из немецкого Bank) и банка?

Значение «банк» возникло из итальянского banco (видимо, через значения «скамья», «прилавок»), заимствованного из германских языков. Для германского корня Wiktionary.org приводит два слова: *bankiz (откуда англ. bench) и *bankô (откуда англ. bank «берег реки», «отмель»). К модели «изгиб» → «возвышение; вздутие», «возвышение» → «скамья; берег, холм» вполне подходило бы дополнение вида «вздутие» → «медицинская банка» → «банка», в предположении, что медицинские банки пришли в Россию из немецкоговорящих краёв. Но увы в германских языках производные bank- для этой цели, похоже, не используются, а медицинские банки известны весьма продолжительное время.

Ну ладно, посмотрим, что, например, написано в Wiktionary о русском слове: «From Proto-Slavic *banъka, from *banja +‎ *-ъka.» Вот это поворот! Банка — уменьшительное от баня. Но ведь в русском и так есть такое слово: собственно банька. Откуда мог взяться вариант с твёрдым -н-? Сомнительно, чтобы вариант *банъка где-то был зафиксирован непосредственно, но если посмотреть ЭССЯ (т.1, с.152), то, оказывается, в ряде славянских языков есть вариант с мягким -н- и значением, идентичным или близким к «банка». Более того, если посмотреть список слов из предыдущей статьи *ban’a (баня), то там можно найти аналогичное значение («сосуд для воды» или даже «все пузатое») без уменьшительного суффикса. Твёрдость -н- не объясняется, но я могу предложить следующие возможные объяснения: а) в данном значение слово могло быть заимствовано в русский из другого славянского языка, где уже произносилось «твёрдо»; б) когда-то был (но куда-то сплыл) также вариант слова баня с твёрдым исходом; в) повлияло другое слово, например, жбан. Последний вариант можно «вывернуть наизнанку»: чьбанъжбанкабанка под влиянием слова баня, тем более, в ЭССЯ в числе производных праславянского *čьbanъ дано такое слово как слвц. dbanka.

Итак, с этимологической позиции слово банка в значении «сосуд для воды» имеет корень бан-, и я был не прав, поставив его в один ряд со словом доска. (Однако, в русском wiktionary.org указан корень банк-.) Впрочем, следует отметить существование позднего заимствования германского слова «банка» (в значении «скамья» и др.), в котором должен быть корень банк-, но склоняется это слово тоже так, как если бы в нём был суффикс -к-.

Доска, desk, discus

Disclaimer: автор дилетант в лингвистике и просит не судить строго.

Как-то довелось мне встретить рекомендацию писать досочка вместо досточка в качестве уменьшительной формы слова доска. Забавно, мне бы и в голову не пришло писать досточка, но правда ли досочка лучше? Вариант досочка предполагает, что в слове доска корень дос и суффикс к, с чем сложно увязать наличие слов дощатый и дощечка. Ага, вопрос «что правильно досточка или досочка?» некорректен! Правильно — дощечка! Здесь мы имеем чередование в корне: доск/дощ (щ- — результат палатализации ск-). Но точно такое же чередование возможно для корней, оканчивающихся на -ст, что объясняет возникновение варианта досточка то ли из дощатый, то ли из дощечка (депалатализация при замене «смягчающей» -е- на «несмягчающую» -о-). Есть закрепившийся в литературном языке результат аналогичного перехода: пуститьпущатьпускать. Впрочем, есть и примеры закреплённого переосмысления части корня в качестве суффикса (в пользу досочки): банка (корень банк-, upd: см. также) — баночка, а также фляжка (корень фляжк-, ср. англ. flask, нем. Flasche) — фляга (ещё и депалатализация -ж-, почему там вообще -ж-, а не -ш-?). В общем-то, уже форма мн.ч. р.п. досок предполагает переосмысление -к- как суффикса.

Ну ладно, корень доск-. А какова его этимология? Те, кто учил английский, могли заметить подозрительно похожее слово: desk. Этимологи давно заметили это сходство (как и подобные слова в других германских языках), и общим мнением является заимствование из германского *disk(az) славянского *dъska (и в словаре Фасмера, и в ЭССЯ, т.5, с.184), продолжаемого др.-русск. и ст.-слав. дъска. Германское слово заимствовано из латинского discus, а то — в свою очередь, из др.-греческого. Правда, есть технические сложности: германское слово мужского рода и должно было в славянском дать **dьskъ (также другая гласная в корне), чего мы не наблюдаем.

Впрочем, история этого корня в германских языках тоже не так проста: др.-англ. disc [diʃ] не является прообразом desk. В английском есть его прямое продолжение: dish. Это первый «продукт» discus в английском. Третий, поздний, продукт: disc или disk (собственно «диск», с орфографией всё никак не могут определиться). А desk — второй продукт, en.wiktionary.org предоставляет его этимологию: «From Middle English deske, desque, from Medieval Latin desca, modified from Old Italian desco, from Latin discus.» По-моему, это прекрасно. Вариант desca («дэска») вполне годится в качестве прообраза слова дъска. Тем более, что мы встречаем его, например, в ряде славянских языков (по сведениям wiktionary и ЭССЯ). Тут бы хорошо знать момент фиксации слова desca, конечно. Но, думаю, вполне вероятно то, что оно могло появиться много раньше VIII в.