Интеграция 2D данных по нерегулярной сетке в Python - PullRequest
Новая
       15

Интеграция 2D данных по нерегулярной сетке в Python

2 голосов
/ 06 марта 2019

Итак, у меня есть 2D-функция, которая нерегулярно дискретизируется по домену, и я хочу вычислить объем под поверхностью.Данные организованы в виде [x,y,z] на простом примере: 

def f(x,y):
    return np.cos(10*x*y) * np.exp(-x**2 - y**2)

datrange1 = np.linspace(-5,5,1000)
datrange2 = np.linspace(-0.5,0.5,1000)

ar = []
for x in datrange1:
    for y in datrange2:
        ar += [[x,y, f(x,y)]]


for x in xrange2:
    for y in yrange2:
        ar += [[x,y, f(x,y)]] 

val_arr1 = np.array(ar)

data = np.unique(val_arr1)


xlist, ylist, zlist = data.T

, где np.unique сортирует данные в первом столбце, а затем во втором.Данные упорядочены таким образом, что мне нужно более интенсивно производить выборку вокруг начала координат, поскольку есть острая особенность, которую необходимо разрешить.

Теперь я задумался о построении двухмерной интерполяционной функции с использованием scipy.interpolate.interp2d, а затем об интегрированиичерез это dblquad.Оказывается, это не только не элегантно и медленно, но также выдает ошибку: 

RuntimeWarning: No more knots can be added because the number of B-spline
coefficients already exceeds the number of data points m.

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

1 Ответ

0 голосов
/ 06 марта 2019

Если вы можете произвести выборку данных с достаточно высоким разрешением вокруг интересующей вас функции, то, более разреженно, везде, где возникает проблема, становится определение области под каждой выборкой. Это легко сделать с обычными прямоугольными выборками, и, вероятно, это можно сделать пошагово с шагом разрешения вокруг начала координат. Подход, который я использовал, состоит в том, чтобы генерировать двумерные клетки Вороного для каждого образца, чтобы определить их площадь. Я извлек большую часть кода из этого ответа, поскольку в нем уже были почти все необходимые компоненты. 

import numpy as np
from scipy.spatial import Voronoi

#taken from: # https://stackoverflow.com/questions/28665491/getting-a-bounded-polygon-coordinates-from-voronoi-cells
#computes voronoi regions bounded by a bounding box
def square_voronoi(xy, bbox): #bbox: (min_x, max_x, min_y, max_y)
    # Select points inside the bounding box
    points_center = xy[np.where((bbox[0] <= xy[:,0]) * (xy[:,0] <= bbox[1]) * (bbox[2] <= xy[:,1]) * (bbox[2] <= bbox[3]))]
    # Mirror points
    points_left = np.copy(points_center)
    points_left[:, 0] = bbox[0] - (points_left[:, 0] - bbox[0])
    points_right = np.copy(points_center)
    points_right[:, 0] = bbox[1] + (bbox[1] - points_right[:, 0])
    points_down = np.copy(points_center)
    points_down[:, 1] = bbox[2] - (points_down[:, 1] - bbox[2])
    points_up = np.copy(points_center)
    points_up[:, 1] = bbox[3] + (bbox[3] - points_up[:, 1])
    points = np.concatenate((points_center, points_left, points_right, points_down, points_up,), axis=0)
    # Compute Voronoi
    vor = Voronoi(points)
    # Filter regions (center points should* be guaranteed to have a valid region)
    # center points should come first and not change in size
    regions = [vor.regions[vor.point_region[i]] for i in range(len(points_center))]
    vor.filtered_points = points_center
    vor.filtered_regions = regions
    return vor

#also stolen from: https://stackoverflow.com/questions/28665491/getting-a-bounded-polygon-coordinates-from-voronoi-cells
def area_region(vertices):
    # Polygon's signed area
    A = 0
    for i in range(0, len(vertices) - 1):
        s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
        A = A + s
    return np.abs(0.5 * A)

def f(x,y):
    return np.cos(10*x*y) * np.exp(-x**2 - y**2)

#sampling could easily be shaped to sample origin more heavily
sample_x = np.random.rand(1000) * 10 - 5 #same range as example linspace
sample_y = np.random.rand(1000) - .5
sample_xy = np.array([sample_x, sample_y]).T

vor = square_voronoi(sample_xy, (-5,5,-.5,.5)) #using bbox from samples
points = vor.filtered_points
sample_areas = np.array([area_region(vor.vertices[verts+[verts[0]],:]) for verts in vor.filtered_regions])
sample_z = np.array([f(p[0], p[1]) for p in points])

volume = np.sum(sample_z * sample_areas)

Я точно не проверял это, но принцип должен работать, и математика проверяется.

...