Блог

Я пытаюсь преобразовать фрагмент кода python в Julia 1.1.0. Код python:

 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 
  

и моя попытка кода Julia:

 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
  

когда я запускаю коды Julia, для 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, теперь я могу преобразовать символьные значения в 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 
  

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

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

Что касается первого, вы сталкиваетесь с проблемой Julia с целочисленными основаниями и нелитеральными, отрицательными целочисленными степенями. Вы можете изменить эту строку, чтобы сделать степени символьными, например, 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 выше.