Документация разработчика

Веса зависят от интерполяционной схемы и позиции x, y..., но не от коэффициентов интерполируемого массива. Поэтому возможно множество ситуаций, когда может потребоваться повторно использовать ранее вычисленные веса, и пакет Interpolations.jl спроектирован с учетом этого факта.

Ключевой концепцией на этом шаге является тип Interpolations.WeightedIndex, но приводить здесь всю информацию из его подробного docstring не имеет смысла. Достаточно будет добавить, что WeightedIndex — это по сути абстрактный тип с двумя конкретными подтипами:

Способ вычисления WeightedIndex зависит от интерполяционной схемы (например, Linear или Quadratic), а также от того, вычисляются ли значения, градиенты (gradient) или гессианы (hessian). Порядок обработки производных более подробно описывается ниже.

Общие массивы AbstractArray можно индексировать с помощью индексов WeightedIndex, и в результате получается интерполированное значение. Иными словами, конечным результатом является просто itp.coefs[wis...], где wis — кортеж индексов WeightedIndex. Чтобы такая перегрузка действовала, мы заключаем coefs в InterpGetindex, то есть InterpGetindex(itp.coefs)[wis...]

Производные по определенной оси можно вычислять, просто подставляя вместо компонента wis тот, который предназначен для вычисления производных, а не значений.

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

Если оставить в стороне такие детали, как проверка границ, главный вызов — это вызов метода Interpolations.weightedindexes:

Можно заметить, что каждый компонент wis соответствует определенной позиции: 1,2, 1,4 и 1,7 соответственно. Мы можем обратиться к A по индексу wis, и в результате будет возвращено значение itp(x...), которое в данном случае будет таким:

Было вычислено значение itp при x..., потому что метод weightedindexes был вызван с единственной функцией, Interpolations.value_weights (а это означает, что «для весов требовалось вычислить значение»).

Для вычисления производных мы также передаем дополнительные функции, например Interpolations.gradient_weights:

В этом случае можно заметить, что wis — 3-элементный кортеж из 3-элементных кортежей. A[wis[i]...] можно использовать для вычисления i-го компонента градиента.

Если внимательно посмотреть на каждый элемент в wis, можно заметить, что в каждом «внутреннем» 3-элементном кортеже повторяются два из трех элементов компонента wis, полученного выше при вызове weightedindexes только с функцией value_weights. wis[1] заменяет первый элемент взвешенным индексом с весами (-1.0, 1.0), что соответствует вычислению углового коэффициента по измерению. Аналогичным образом, wis[2] и wis[3] заменяют соответственно второй и третий индексы значений таким же вычислением углового коэффициента.

Гессиан вычисляется очень похожим образом. Отличие в том, что иногда необходимо заменить два разных индекса или один индекс на набор весов, соответствующих второй производной.

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

Код для выполнения такой замены несколько сложен из-за необходимости поддерживать произвольную размерность, чтобы в Julia можно было успешно выполнять вывод типов. В нем применяются особые операции с кортежами, иногда называемые «программированием кортежей в стиле Лиспа». Дополнительные советы, касающиеся такого программирования, можно найти на форуме обсуждения Julia. Реализовать такую замену можно было бы и с помощью генерируемых функций, но это существенно увеличило бы время компиляции и могло бы привести к проблемам с «возрастом мира».

В настоящее время Interpolations.jl поддерживает использование интерполянтов на GPU посредством трансляции.

Базовый рабочий процесс выглядит так:

Для этой цели в ITP <: AbstractInterpolation должен быть определен собственный метод Adapt.adapt_structure(to, itp::ITP), который создает новый объект ITP с адаптированными полями (adapt(to, itp.fieldname)) itp. Адаптацию полей можно пропустить, если известно, что объект уже совместим с GPU, как в случае, например, с диапазоном isbit.

Кроме того, для всех совместимых с GPU типов AbstractInterpolation должен быть определен собственный метод Interpolations.root_storage_type. Эта функция позволяет изменять механизм трансляции, перегружая BroadcastStyle по умолчанию. Дополнительные сведения см. в разделе Настройка трансляции.