NaN and division by 0 in solve_for

[MarkNahabedian]

I originally posted this to the Julia #helpdesk slack on October 9th, but didn't get an answer then. I'm feeling like it's a bug.

I'm having trouble using Symbolics.solve_for.

I'm working with vectors in three dimensional space. Because @variables is for real numbers, I've defined

eltname(name, i) = Symbol("($name)_$i")

function symbolic_vector(name::Symbol, n::Integer)
    map(variable, map(i -> eltname(name, i), 1:n))
end

This gives me a vector of n Symbolics variables.

I think this is the approach suggested to me by (I think) Sashi at a recent CAJUN.

As a convenience, when I want to substitute numeric values for one of these vectors, I define

function vsubs(v::Vector{Symbolics.Num}, values...)
    Pair.(v, [1,2,3])
end

I should be able to wrap the results of multiple calls to vsubs in a Dict to pass to `Symbolics.substitute'.

I can define a function that will give me the parametric formula for points on a line if I know two points on the line:

function line(p1, p2, param)
    p1 + param * (p2 - p1)
end

Suppose I want to find the point where two lines intersect:

let
    P1 = symbolic_vector(:P1, 3)
    P2 = symbolic_vector(:P2, 3)
    @variables r
    l1 = line(P1, P2, r)

    P3 = symbolic_vector(:P3, 3)
    P4 = symbolic_vector(:P4, 3)
    @variables s
    l2 = line(P3, P4, s)

    eq = Equation(l1, l2)
    solve_for(eq, s)
end

Yes, there is no guarantee that these lines intersect, but I don't see how that is a justification to divide by zero:

Symbolics.Num[
    (1//0)*(var"(P1)_1"
            - var"(P3)_1"
            + (-var"(P1)_1"
               + var"(P2)_1")*r
            - (-var"(P3)_1"
               + var"(P4)_1")*s),
    (1//0)*(var"(P1)_2"
            - var"(P3)_2"
            + (-var"(P1)_2"
               + var"(P2)_2")*r
            - (-var"(P3)_2"
               + var"(P4)_2")*s),
    (1//0)*(var"(P1)_3"
            - var"(P3)_3"
            + (-var"(P1)_3"
               + var"(P2)_3")*r
            - (-var"(P3)_3"
               + var"(P4)_3")*s)]
[1:3][1:3]

If I substitute concrete values for my points

let
    P1 = symbolic_vector(:P1, 3)
    P2 = symbolic_vector(:P2, 3)
    @variables r
    l1 = line(P1, P2, r)

    P3 = symbolic_vector(:P3, 3)
    P4 = symbolic_vector(:P4, 3)
    @variables s
    l2 = line(P3, P4, s)

    subs = Dict([
        vsubs(P1, 0, 0, 0)...,
        vsubs(P2, 10, 10, 0)...,
        vsubs(P3, 4, 0, 0)...,
        vsubs(P4, 0, 4, 0)...
    ])

    eq = Equation(l1, l2)
    # Expected solution: r = 0.2, s = 0.5.
    eq2 = substitute(eq, subs)    

    solve_for(eq2, s)
end

I can just look at eq2 and verify that r = 0.2 and s = 0.5, but the result of solve_for still has zero denominators:

Symbolics.Num[(1//0)*(-4 + 10r + 4s), (1//0)*(10r - 4s), NaN][1:3][1:3]

Am I doing something wrong?