RBF +++ Python-пакет, содержащий инструменты для применения функций с радиальной базисной функцией (RBF). Применения включают интерполяцию/сглаживание рассеянных данных и решение уравнений в частных производных (PDE) на нерегулярных областях. Этот пакет был значительным образом вдохновлен книгами "Meshfree Approximation Methods with Matlab" Грегори Фассхауэра и "A Primer on Radial Basis Functions with Applications to the Geosciences" Бенкта Форнберга и Наташи Флайер. Полная документация для этого пакета доступна `здесь `_. Особенности ============ * Функции для оценки RBF и их точных производных. * Класс для интерполяции RBF, который используется для интерполяции и сглаживания рассеянных, шумных, N-мерных данных. * Алгоритм для генерации весов Radial Basis Function Finite Difference (RBF-FD). Этот алгоритм используется для решения больших систем уравнений PDE на нерегулярных областях. * Алгоритм генерации узлов, который может быть использован для решения PDE с помощью спектрального метода RBF или метода RBF-FD. * Абстракция для гауссовских процессов. Гауссовские процессы здесь主要用于高斯过程回归（GPR），这是一种非参数贝叶斯插值/平滑方法。 * Генератор последовательностей Халтона. * Функции вычислительной геометрии (например, тестирование внутренних точек многоугольника) для двумерного и трёхмерного пространства. Установка ========= RBF требует следующие Python-пакеты: numpy, scipy, sympy, cython и rtree. Эти зависимости могут быть удовлетворены базовым Anaconda Python-пакетом (https://www.continuum.io/downloads). Загрузка пакета RBF .. code-block:: bash $ git clone http://github.com/treverhines/RBF.git Компиляция и установка .. code-block:: bash $ cd RBF $ python setup.py install Тестирование всех функций на работоспособность .. code-block:: bash $ cd test $ python -m unittest discover Демонстрация ============ Сглаживание рассеянных данных ------------------------------ .. code-block:: python ''' В этом примере мы генерируем синтетические рассеянные данные с добавленным шумом, а затем аппроксимируем их сглаженным RBF-интерполянтом. Интерполянт в этом примере эквивалентен тонкому пластинчатому сингуляру. ''' import numpy as np from rbf.interpolate import RBFInterpolant import matplotlib.pyplot as plt np.random.seed(1) # точки наблюдения x_obs = np.random.random((100, 2)) # значения в точках наблюдения u_obs = np.sin(2 * np.pi * x_obs[:, 0]) * np.cos(2 * np.pi * x_obs[:, 1]) u_obs += np.random.normal(0.0, 0.1, 100) # создаем интерполянт тонкого пластинчатого сингуляра, где данные считаются # шумовыми I = RBFInterpolant(x_obs, u_obs, sigma=0.1, phi='phs2', order=1) # создаем точки интерполяции и оцениваем интерполянт x1, x2 = np.linspace(0, 1, 200), np.linspace(0, 1, 200) x_itp = np.reshape(np.meshgrid(x1, x2), (2, 200 * 200)).T u_itp = I(x_itp) # строим результаты plt.tripcolor(x_itp[:, 0], x_itp[:, 1], u_itp, vmin=-1.1, vmax=1.1, cmap='viridis') plt.scatter(x_obs[:, 0], x_obs[:, 1], s=100, c=u_obs, vmin=-1.1, vmax=1.1, cmap='viridis', edgecolor='k') plt.xlim((0.05, 0.95)) plt.ylim((0.05, 0.95)) plt.colorbar() plt.tight_layout() plt.show() .. figure:: docs/figures/interpolate.a.png График, сгенерированный вышеуказанным кодом. Наблюдения показаны как точки рассеяния, а гладкий интерполянт — как цветовое поле.Решение уравнений в частных производных --------------------------------------- Существуют два метода решения уравнений в частных производных с помощью функций радиального базиса (RBF): спектральный метод и метод RBF-FD. Спектральный метод известен своей выдающейся точностью; однако он применим только для малых задач и требует хорошего выбора параметра формы. Метод RBF-FD привлекателен тем, что он может использоваться для больших задач, не требует настройки параметра формы (при условии использования полигармонических сплайнов для генерации весов), и более высокий порядок точности может быть достигнут путем увеличения размера шаблона или увеличения порядка полинома, используемого для генерации весов. Кратко, метод RBF-FD всегда предпочтителен перед спектральным методом RBF. Примеры двух методов приведены ниже... code-block:: python ''' В этом примере мы решаем уравнение Пуассона на L-образной области с фиксированными граничными условиями. Мы используем мультиквадратический RBF (`mq`) ''' import numpy as np from rbf.basis import mq from rbf.pde.geometry import contains from rbf.pde.nodes import poisson_disc_nodes import matplotlib.pyplot as plt # Определяем область задачи с помощью линий. vert = np.array([[0.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0], [1.0, 2.0], [0.0, 2.0]]) smp = np.array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]]) spacing = 0.07 # приближенное расстояние между узлами eps = 0.3/spacing # параметр формы # генерируем узлы. `nodes` это (N, 2) массив типа float, `groups` это # словарь, идентифицирующий, к какому группе принадлежит каждый узел nodes, groups, _ = poisson_disc_nodes(spacing, (vert, smp)) N = nodes.shape[0] # создаем "левую часть" матрицу A = np.empty((N, N)) A[groups['interior']] = mq(nodes[groups['interior']], nodes, eps=eps, diff=[2, 0]) A[groups['interior']] += mq(nodes[groups['interior']], nodes, eps=eps, diff=[0, 2]) A[groups['boundary:all']] = mq(nodes[groups['boundary:all']], nodes, eps=eps) # создаем "правую часть" вектор d = np.empty(N) d[groups['interior']] = -1.0 # вынуждающая сила d[groups['boundary:all']] = 0.0 # граничное условие ''' # Решаем для коэффициентов RBF coeff = np.linalg.solve(A, d) # интерполируем решение на сетке xg, yg = np.meshgrid(np.linspace(0.0, 2.02, 100), np.linspace(0.0, 2.02, 100)) points = np.array([xg.flatten(), yg.flatten()]).T u = mq(points, nodes, eps=eps).dot(coeff) # маскировать точки, находящиеся вне области u[~contains(points, vert, smp)] = np.nan # сворачиваем решение в сетку ug = u.reshape((100, 100)) # создаем контурную карту решения fig, ax = plt.subplots() p = ax.contourf(xg, yg, ug, np.linspace(0.0, 0.16, 9), cmap='viridis') ax.plot(nodes[:, 0], nodes[:, 1], 'ko', markersize=4) for s in smp: ax.plot(vert[s, 0], vert[s, 1], 'k-', lw=2) ax.set_aspect('equal') ax.set_xlim(-0.05, 2.05) ax.set_ylim(-0.05, 2.05) fig.colorbar(p, ax=ax) fig.tight_layout() plt.show() .. figure:: docs/figures/basis.a.png :align: center :width: 50% :alt: Результат решения Результат решения.. code-block:: python ''' В этом примере мы решаем уравнение Пуассона на L-образной области с фиксированными граничными условиями. Мы используем метод RBF-FD. Метод RBF-FD предпочтителен перед спектральным методом RBF, так как он масштабируем и не требует от пользователя указания параметра формы (при условии использования нечетного порядка полигармонических сплайнов для вычисления весов). ''' import numpy as np from scipy.sparse import coo_matrix from scipy.sparse.linalg import spsolve import matplotlib.pyplot as plt from rbf.sputils import add_rows from rbf.pde.fd import weight_matrix from rbf.pde.geometry import contains from rbf.pde.nodes import poisson_disc_nodes # Определяем область задачи с помощью линий. vert = np.array([[0.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0], [1.0, 2.0], [0.0, 2.0]]) smp = np.array([[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]]) # Размер узла составляет 0.03 в точке [1, 1] и увеличивается по мере # удаления от этой точки def spacing(x): return 0.04 + 0.08 * np.linalg.norm(x - 1.0, axis=1) n = 25 # размер шаблона. Увеличение этого значения обычно улучшает точность phi = 'phs3' # радиальная базовая функция, используемая для вычисления весов. # Нечетные порядки полигармонических сплайнов (например, phs3) # всегда хорошо работают для меня и не требуют от пользователя # настройки параметра формы. Используйте более высокие порядки # полигармонических сплайнов для более высоких порядков уравнений. order = 2 # Порядок добавленных полиномов. Это должно быть не меньше # порядка решаемого уравнения ПДУ (2 в данном случае). # Большие значения могут улучшить точность. # сгенерировать узлы nodes, groups, _ = poisson_disc_nodes(spacing, (vert, smp)) N = nodes.shape[0] # создать матрицу "левой части" # создать компонент, который оценивает PDE A_interior = weight_matrix(nodes[groups['interior']], nodes, n, diffs=[[2, 0], [0, 2]], phi=phi, order=order) # создать компонент для фиксированных граничных условий A_boundary = weight_matrix(nodes[groups['boundary:all']], nodes, 1, diffs=[0, 0]) # Добавить компоненты в соответствующие строки `A` A = coo_matrix((N, N)) A = add_rows(A, A_interior, groups['interior']) A = add_rows(A, A_boundary, groups['boundary:all']) # создать вектор "правой части" d = np.zeros((N,)) d[groups['interior']] = -1.0 d[groups['boundary:all']] = 0.0 # найти решение в узлах u_soln = spsolve(A, d) # Создать сетку для интерполяции решения xg, yg = np.meshgrid(np.linspace(0.0, 2.02, 100), np.linspace(0.0, 2.02, 100)) points = np.array([xg.flatten(), yg.flatten()]).T # Можно использовать любой метод рассеянной интерполяции (например, # scipy.interpolate.LinearNDInterpolator). Здесь мы используем метод RBF-FD # для интерполяции с высоким порядком точности u_itp = weight_matrix(points, nodes, n, diffs=[0, 0]).dot(u_soln) # маскировать точки вне области u_itp[~contains(points, vert, smp)] = np.nan ug = u_itp.reshape((100, 100)) # вернуть в сетку # создать график контурного решения fig, ax = plt.subplots() p = ax.contourf(xg, yg, ug, np.linspace(-1e-6, 0.16, 9), cmap='viridis') ax.plot(nodes[:, 0], nodes[:, 1], 'ko', markersize=4) for s in smp: ax.plot(vert[s, 0], vert[s, 1], 'k-', lw=2) ax.set_aspect('equal') ax.set_xlim(-0.05, 2.05) ax.set_ylim(-0.05, 2.05) fig.colorbar(p, ax=ax) fig.tight_layout() plt.show() # figure:: docs/figures/fd.i.png # Заметка: Текст для описания фигуры отсутствовал в исходном документе.