Scipy using PyCall in Julia

Question

I have a moderately large system of non-linear equations that I want to solve using scipy.optimize in Julia. The problem is that I store the equations in a vector before I pass them to the solver and PyCall doesn't accept that. For example these methods both work:

using PyCall
@pyimport scipy.optimize as so

function F(x)
 f1=1- x[1] - x[2]
 f2=8 - x[1] - 3*x[2]
 return f1, f2
end

x0 = [1,1]
x = so.fsolve(F, x0)

function G(x)
 f=[1 - x[1] - x[2],
    8 - x[1] - 3*x[2]]
 return f
end

x0 = [1,1]
x = so.fsolve(G, x0)

However this does not:

function H(x)
 f[1]=1 - x[1] - x[2]
 f[2]=8 - x[1] - 3*x[2]
 return f
end

x0 = [1,1]
x = so.fsolve(H, x0)

Neither does this:

function P(x)
 f[1]= 1 - x[1] - x[2]
 f[2]= 8 - x[1] - 3*x[2]
 return f[1], f[2]
end

x0 = [1,1]
x = so.fsolve(P, x0)

I don't think it's feasible to not use a loop because of the nature of the problem. Is there any way to return the vector in a way that fsolve can accept it?

Answer 1

The second two methods never create f which is the issue. You have to create the array first.

function H(x)
 f = similar(x)
 f[1]=1 - x[1] - x[2]
 f[2]=8 - x[1] - 3*x[2]
 return f
end

which will automatically match the size and type of x, or you can use the constructor:

function H(x)
 f = Vector{Float64}(2)
 f[1]=1 - x[1] - x[2]
 f[2]=8 - x[1] - 3*x[2]
 return f
end

However you want to do it, you need to make the array. P has the same problem.

Also, you should checkout NLSolve.jl. It allows a pre-allocated form:

function H(x,f)
 f[1]=1 - x[1] - x[2]
 f[2]=8 - x[1] - 3*x[2]
 return nothing
end

which should allocate less and do better. Roots.jl is another good Julia option.