Understanding and translating this Julia code to Python

Question

I'm attempting to translate this Julia code, an implementation of the multivariate Newton-Raphson method, to Python, and expand it to accept a system of four equations in four variables.

The code:

using SymPy
x, y = Sym("x, y")
function SymPyDerivatives(f::String, g::String)
    #
    # Returns the string representations of all
    # partial derivatives of the expression in x and y,
    # given in the strings f and g.
    #
    formf = parse(f)
    evalformf = eval(formf)
    fx = diff(evalformf, x)
    fy = diff(evalformf, y)
    formg = parse(g)
    evalformg = eval(formg)
    gx = diff(evalformg, x)
    gy = diff(evalformg, y)
    return [string(fx) string(fy); string(gx) string(gy)]
end

function NewtonStep(fun::Array{SymPy.Sym,1},
    jac::Array{SymPy.Sym,2},
    x0::Float64, y0::Float64)
    #
    # Runs one step with Newton’s method
    #
    valfun = -SymPyFun(fun, x0, y0)
    nfx = norm(valfun)
    valmat = SymPyMatrixEvaluate(jac, x0, y0)
    update = valmat\valfun
    ndx = norm(update)
    x1 = x0 + update[1]
    y1 = y0 + update[2]
    sfx = @sprintf("%.2e", nfx)
    sdx = @sprintf("%.2e", ndx)
    sx1 = @sprintf("%.16e", x1)
    sy1 = @sprintf("%.16e", y1)
    println(" $sfx $sdx $sx1 $sy1 ")
    return [x1, y1]
end

function main()
#
# Prompts the user for expressions in x and y,
# calls SymPy to take the derivatives with respect
# to x and y of the two given expressions.
#
    println("Reading two expressions, hit enter for examples")
    print("Give a first expression in x and y : ")
    one = readline(STDIN)
    if one == ""
        one = "exp(x) - y"
        two = "x*y - exp(x)"
        x0, y0 = 0.9, 2.5
    else
        print("Give a second expression in x and y : ")
        two = readline(STDIN)
        print("Give a value for x0 : ")
        x0 = parse(Float64(readline(STDIN)))
        print("Give a value for y0 : ")
        y0 = parse(Float64(readline(STDIN)))
    end
    derivatives = SymPyDerivatives(one, two)
    println("The Jacobian matrix :")
    println(derivatives)
    fx = derivatives[1,1]
    fy = derivatives[1,2]
    gx = derivatives[2,1]
    gy = derivatives[2,2]
    jacmat = SymPyMatrix(fx, fy, gx, gy)
    valmat = SymPyMatrixEvaluate(jacmat, x0, y0)
    vecfun = SymPyExpressions(one, two)
    for i=1:5
        xsol, ysol = NewtonStep(vecfun, jacmat, x0, y0)
        x0, y0 = xsol, ysol
    end
end

main()

My problem being, I've effectively never used Julia and, while I understand Newton's method fairly well, I don't quite follow this code.

Here's what I've got so far:

Inside function NewtonStep()

Not sure what's going on here: valfun = -SymPyFun(fun, x0, y0)

We take the norm of that matrix nfx = norm(valfun)

Not sure what's going on here: valmat = SympPyMatrixEvaluate(jax,c, x0, y0)

Update the matrix:

update = valmat\valfun
ndx = norm(update)

Not sure what's going on afterwards. And why is function SymPyDerivatives() never used? Surely it must be somehow.

Not really sure what's going on inside function SymPyDerivatives(), outside of the lines taking the partial derivatives.

Update:

Just trying to pull this apart piece-by-piece, this is what I have so far:

from sympy import symbols, diff
import numpy as np

class Localize:

    receivers: list

    def jacobian(f, g, h, k):
        x, y, z = symbols('x y z', real=True)

        fx = diff(f, x);    gx = diff(g, x)
        fy = diff(f, y);    gy = diff(g, y)

        hx = diff(h, x);    kx = diff(k, x)
        hy = diff(h, y);    ky = diff(k, y)

    def newtonStep(x0, x1):

        # Magic happens

        x1 = x0 + update[1]
        y1 = y0 + update[2]

Can't seem to find any documentation on SymPyFun(), which is bizarre.

Edit:

According to my own research and that of @steamsy, it would seem that SymPyFun() does not exist, which is odd. Perhaps someone has some sense of what this was intended to accomplish, and what that might look like in Python? Not really sure what it's supposed to do here.

Update:

I was able to find function SymPyFun()

function SymPyFun(fun::Array{SymPy.Sym,1},
                  xval::Float64,yval::Float64)
#
# Given an array of SymPy expressions,
# evaluates the expressions at (xval, yval).
#
    result1 = Float64(subs(subs(fun[1], x, xval), y, yval))
    result2 = Float64(subs(subs(fun[2], x, xval), y, yval))
    return [result1; result2]
end

Answer 1

Can't seem to find any documentation on SymPyFun(), which is bizarre.

Because it doesn't exist. This documentation on SymPy.jl and my scraping of the GitHub source code never mentions SymPyFun.

I'm assuming you're following along with this PDF from UIC, this is the only thing I could find that references the method.