Вычисление символического выражения в SymPy и преобразование их в числовое значение в Julia - PullRequest
Новая
       17

Вычисление символического выражения в SymPy и преобразование их в числовое значение в Julia

0 голосов
/ 03 апреля 2019

Я пытаюсь преобразовать часть кода Python в Julia 1.1.0. Код питона: 

import numpy as np
import sympy as sp
x, y = sp.symbols('x y')
data = np.random.randn(1000, 2)
a,b = data[:,0], data[:,1]
M = len(data[:,0]) 
m = 5
n = round(m*(m+1)/2) 
L = np.zeros((M,n)) 
l = sp.zeros(1,n)
k = 0
for i in range(1,m+1):
    for j in range(1,i+1):
        l[0,k]=((i-j)*(i-j-1)*x**(i-j-2)*y**(j-1))

        f = sp.lambdify([x,y], l[0,k], "numpy")
        L[:,k] = np.squeeze(f(a,b))
        k=k+1

и моя попытка Юлии код: 

using SymPy
data = rand(1000,2)
a = data[:,1];
b = data[:,2];
M = length(data[:,1]) 
x = symbols("x");
y = symbols("y");
m = 5;
n = Int(m*(m+1)/2) 
L = zeros(M,n)
l = zeros(1,n)
k = 1;
for i in 1:m
      for j in 1:i
            l[1,k] = ((i-j)*(i-j-1)*x^(i-j-2)*y^(j-1))
            f = l[1,k]
            L[:,k] = f.subs([(x, a), (y, b)])
            k=k+1  
    end
end

когда я запускаю коды Джулии, для l[1,k] я получил следующую ошибку 

DomainError with -2:
Cannot raise an integer x to a negative power -2.
Convert input to float.
Stacktrace:
 [1] throw_domerr_powbysq(::Sym, ::Int64) at ./intfuncs.jl:173
 [2] power_by_squaring(::Sym, ::Int64) at ./intfuncs.jl:196
 [3] ^(::Sym, ::Int64) at ./intfuncs.jl:221
 [4] top-level scope at ./In[80]:14

Кроме того, я не уверен насчет следующих кодов 

    f = l[1,k]
    L[:,k] = f.subs([(x, a), (y, b)])

Буду признателен, если кто-нибудь поможет мне перевести коды python в коды julia.

Обновление:

Основываясь на сообщении @jverzani, теперь я могу конвертировать значения sym в float с помощью следующих кодов 

using SymPy
data = rand(1000,2)
a = data[:,1];
b = data[:,2];
M = length(data[:,1])
x = symbols("x");
y = symbols("y");
m = 5;
n = Int(m*(m+1)/2)
LL = zeros(M,n)
L = zeros(Sym, M,n)
l = zeros(Sym, 1,n)
k = 1;
for i in 1:m
      for j in 1:i
            l[1,k] = ((i-j)*(i-j-1)*x^Sym(i-j-2)*y^Sym(j-1))
            f = l[1,k]
            L[:,k] .= f.subs([(x, a), (y, b)])
            global k=k+1
    end
end
for s in 1:M
    for r in 1:n
        LL[s,r] = float(subs(L[s,r],(x,a[s]),(y,b[s])))
    end
end

Но на этот раз коды очень медленные. Как я могу оптимизировать коды.

1 Ответ

0 голосов
/ 03 апреля 2019

SymPy только что получил серьезное изменение, унаследованное от PyCall, что заставляет последний бит с subs работать как напечатано.

Что касается первого, вы сталкиваетесь с проблемой Джулии с целочисленными основаниями и не буквальными, отрицательными целочисленными степенями.Вы можете изменить эту строку, чтобы сделать полномочия символическими, например x^Sym(i-j-2) вместо x^(i-j-2).(Возможно, лучше обойти это.)

Следующая редакция пытается отразить дух кода Python.Он превращается в функцию для проверки скорости на размере N и позволяет избежать необходимости global в назначении k.

using SymPy
import PyCall
np = PyCall.pyimport("numpy") #import numpy as np
const sp = sympy 



function fill_L(N=1000)

    x, y = sp.symbols("x y")
    offset = 1
    data = np.random.randn(N, 2)
    a,b = data[:,0+offset], data[:,1+offset]
    M = length(data[:,0+offset])
    m = 5
    n = round(Int, m*(m+1)/2)
    L = np.zeros((M,n))
    l = sp.zeros(1,n)
    k = 0
    for i in 1:m
        for j in 1:i
            l[1,k+offset]=((i-j)*Sym(i-j-1)*x^Sym(i-j-2)*y^Sym(j-1))
            f = lambdify(l[1,k+offset], (x, y))
            #f = sp.lambdify([x,y], l[0,k], "numpy")
            L[:,k+offset] = f.(a,b)
            k=k+1
        end
    end
    L
end

Производительность в настоящее время является разумной из-за трансляции f выше.

...