Инверсия Scipy interp1d при использовании массива: как найти корень массива выходных функций? - PullRequest
Новая
       4

Инверсия Scipy interp1d при использовании массива: как найти корень массива выходных функций?

0 голосов
/ 05 октября 2018

Допустим, у меня есть следующий 3d-Numpy-массив A: 

array([[[ 1,  2,  3],
        [ 1,  4,  9]],

       [[ 1,  8, 27],
        [ 1, 16, 81]]])

Я хочу интерполировать данные, используя interp1d с axis=2, т.е. интерполировать значения функции 1,2,3, 1,4,9, 1,8,27 и 1,16,81 одновременно с x=np.array([1,2,3]) (значения представляют f (x) = x, x ^ 2, x ^ 3 и x ^ 4 соответственно).

К счастью, этовозможно, так как interp1d может принимать ND-массивы в качестве второго аргумента, если ось интерполяции равна длине x: 

f = interp1d(x,A,kind='linear',axis=2,fill_value='extrapolate')

Пока это работает довольно хорошо, f (1)производит 

array([[1., 1.],
       [1., 1.]]),

f (2) производит 

array([[ 2.,  4.],
       [ 8., 16.]])

и интерполяция f (1.5) производит 

array([[1.5, 2.5],
       [4.5, 8.5]])

т.е. я получаю массив интерполяционных функций 2x2f (x), оцененное в x.

Теперь возникает проблема: я хочу инвертировать эти функции интерполяции, чтобы получить x для конкретного значения функции - поэлементно для каждой записи массива 2x2.

При работе с 1D-функцией g это обычно делается путем нахождения корня функции интерполяции минус запрошенное значение функции, скажем a: 

g_subtracted = lambda x: f(x) - a

и нахождения нуля, дляПример использования scipy.optimize.newton: 

optimize.newton(g_subtracted,1.0)

, где 1.0 - начальное предположение.

Вот мой актуальный вопрос: как мне найти ноль моего массива функций 2x2 interp2d для элемента f? Когда я просто делаю 

f_subtracted = lambda x: f(x) - a

(это работаетхорошо для моего массива функций), а затем 

optimize.newton(f_subtracted,1.0)

я получаю следующую ошибку: 

File "<ipython-input-187-4cf34581a978>", line 1, in <module>
newton(f,1.0)

File "/anaconda3/lib/python3.6/site-packages/scipy/optimize/zeros.py", line 201, in newton
if q1 == q0:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Мне кажется, что optimize.newton не может работать с «массивами функций» каквход для оценки их поэлементно.Кто-нибудь знает способ, как это сделать?

Заранее большое спасибо!

Ответы [ 2 ]

0 голосов
/ 06 октября 2018

В общем, interp1d не очень рекомендуется.Лучше использовать более конкретный интерполятор.Например, вот кубический сплайн: 

In [13]: x = [1, 2, 3]

In [14]: A = np.array([[[ 1,  2,  3],
    ...:         [ 1,  4,  9]],
    ...: 
    ...:        [[ 1,  8, 27],
    ...:         [ 1, 16, 81]]])

In [15]: from scipy.interpolate import CubicSpline

In [17]: spl = CubicSpline(x, A.T)   # <<<< Not the transpose: interpolation is along axis=0

In [19]: spl.roots()
Out[19]: 
array([[array([ 0.]), array([], dtype=float64)],
       [array([], dtype=float64), array([ 1.19999999,  1.20000001])]], dtype=object)

In [20]: r = spl.roots()

In [21]: r.shape
Out[21]: (2, 2)

Однако, обычным способом выполнения таких действий является обратная интерполяция: интерполируйте y против x и вычисляйте значение интерполятора в нуле,Для этого вам нужно будет зациклить входной массив вручную, так как scipy interolators не принимает массив n-dim x.

0 голосов
/ 05 октября 2018

В самых общих чертах это сложная задача NP, и обратное обычно не определяется.Инверсия вашей произвольной функции не гарантируется, даже если бы она была 1-й вместо 4-й.Например, скажем, моя произвольная функция была 1D и задана таблицей A = [1,2,1,0].Я могу легко интерполировать из A при любом x, например, f = interp1d([0,1,2,3],A,kind='linear') дает f(0.5) == 1.5, но также обратите внимание, что f(1.5) == 1.5.Есть два значения x, для которых f(x) равно 1,5!Поэтому инвертирование f (или действительно A) невозможно для всех x.Поэтому вы должны разделить свой домен x на области, в которых существует обратное, и каким-то образом знать, из какого из них вы хотели бы получить обратное.

Теперь представьте более высокие измерения для A и f, как в вашем примере.«Подразделение x» становится проблемой выбора n-мерных объемов, в которых A может быть инвертировано.Вы можете понять, насколько это сложная проблема, для которой нет общего решения.

При этом, если вы знаете (или предполагаете), что ваш f обратим, можно построить инверсию для каждого измерения вашего результата.Таким образом, вместо f_subtracted = lambda x: f(x) - a вы могли бы сделать 

f_subtracted00 = lambda x: f(x)[0,0] - a
f_subtracted01 = lambda x: f(x)[0,1] - a

и т. Д.

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

...