Rewrite Mathematica output in Python

Question

I want to convert the piecewise-function output of a calculation in Mathematica into Python. Taking inspiration from this page for the Mathematica->Python conversion and this page for writing piecewise functions, I have

from numpy import linspace, vectorize, array
from numpy import arctan, log
import matplotlib.pyplot as plt
from sympy.parsing.mathematica import parse

def fun(x,a,b,c):
    # string inside parse('') is my mathematica output
    if a == 0:
        out = parse('a (-I b + x) ArcTan[(b - a)/c]')
    else:
        out = parse('4 c a^2 Log[c^2 + (x - a)^2]')
    return out

a = 0.17
b = 0.44
c = 0.29
x = linspace(0,50,int(1e3))

vfun = vectorize(fun)
y = vfun(x,a,b,c)
yp = 4*c*a**2*log(c**2 + (x - a)**2)

plt.figure()
plt.plot(x,y,'.-')
plt.title('mathematica -> python conversion')

plt.figure()
plt.plot(x,yp,'.-')
plt.title('expected')

plt.show()

The plot looks like:

whereas it should be:

Have I done something wrong when converting Mathematica to Python? Or is there a problem when assigning numerical values to a, b, and c? (Note that this is a MWE, and the Mathematica code that I want to convert is much longer and complicated than what is shown above.)

Answer 1

A much easier solution is using eval. Now I know eval is very dangerous but most of the time when it takes input from user, but here the input is already defined for us.

Now onto the answer.
You are not getting the expected output because your vectorized array does not contain any floats but only contains output from mathematica's parser which return unevaluated string, so we can use .evalf() as answer given by @David, lambdify() which uses eval() internally or we can directly use eval().
Here is the documentation of the methods used : https://docs.sympy.org/latest/tutorial/basic_operations.html

from numpy import linspace, vectorize, array
from numpy import arctan, log
import matplotlib.pyplot as plt
from sympy.parsing.mathematica import MathematicaParser

def fun(x,a,b,c):
    obj = MathematicaParser()
    if a == 0:
        out = obj.parse('a (-I b + x) ArcTan[(b - a)/c]')
    else:
        out = obj.parse('4 c a^2 Log[c^2 + (x - a)^2]')

    return out

a = 0.17
b = 0.44
c = 0.29
x = linspace(0,50,int(1e3))
yp = 4*c*a**2*log(c**2 + (x - a)**2)
vfun = vectorize(fun)
y = vfun(x,a,b,c)
#Here y is a type <class 'numpy.ndarray'> containing 1000 <class 'numpy.str_'>
#['4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#'4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#'4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#'4*c*a**2*log(c**2+(x-a)**2)' '4*c*a**2*log(c**2+(x-a)**2)'
#....
y = eval(y[0]) #y[0] is '4*c*a**2*log(c**2+(x-a)**2)'
#When we evaluate y[0] we again get y as <class 'numpy.ndarray'> conatining 1000 <class 'numpy.float64'>
#Because x is <class 'numpy.ndarray'> it evaluates the first string over
#all the points in x.
#[-0.07309464 -0.07770262 -0.08110382 -0.0828403  -0.08263539 -0.08052339
#-0.07683235 -0.07203573 -0.06659307 -0.06086366 -0.05509179 -0.04942739
#-0.04395413 -0.03871304 -0.03371924 -0.0289728  -0.02446552 -0.02018495
#.....
plt.figure()
plt.plot(x, y,'.-')
plt.title('mathematica -> python conversion')

plt.figure()
plt.plot(x,yp,'.-')
plt.title('expected')

plt.show()

Output :

Answer 2

This is a dirty solution:

import numpy as np
import matplotlib.pyplot as plt
from sympy.parsing.mathematica import mathematica
from sympy import symbols

def fun(x, a, b, c):
    # string inside parse('') is my mathematica output
    if a == 0:
        out = mathematica('a (-I b + x) ArcTan[(b - a)/c]')
    else:
        out = mathematica('4 c a^2 Log[c^2 + (x - a)^2]')
    return out

aa = 0.17
bb = 0.44
cc = 0.29
II = 1
xx = np.linspace(0, 50, int(1e3))

x, a, b, c, I = symbols('x, a, b, c, I')

fun1 = fun(x, a, b, c)
fun2 = fun1.subs({a: aa, c: cc})
print(fun2.evalf())

y_list = []

for item in xx:
    y_list.append(fun2.subs({x:item}).evalf())

print(y_list[:10])

plt.figure()
plt.plot(xx, y_list,'.-')
plt.title('mathematica -> python conversion')
plt.show()

I will explain it later how all of this works.

Update

As you may see, when a == 0, the function is 0, you need to check that.

When you use mathematica, the type of fun1 is a sympy function (sympy.core.mul.Mul), that's why you have to use symbols and fun1.subs() and fun2.evalf(), in other words, you need to know who to use the basic functions of sympy.

The way you evaluate a function in order to plot it, well, you use a list:

y_list = []

for item in xx:
    y_list.append(fun2.subs({x:item}).evalf())

By the way, I am using sympy version 1.4.