symfit special functions, exponential integral

[1000NightsAKnight]

I asked this question first on stackoverflow, though this might be the more appropriate place. My objective function contains exponential integrals, Ei, but this results in a name error. i.e. 'Name Error: name 'Ei' is not defined,' which I have not been able to solve.

import numpy as np
import scipy.optimize
from scipy.special import expi
import numpy as np
import symfit as sf

y1 = np.array(data1)

y3 = np.array(ydata3)

x1 = np.array(xdata1)

x3 = np.array(xdata3)

t1 = 50     # AOM 1st OFF Time
t2 = 100    # AOM 2nd ON Time
t3 = 150    # AOM 2nd OFF Time 
t4 = 450

def mod1(data, a_decay, b_decay, constant, on1_a_amplitude, on1_b_amplitude, ARE_A, ARE_B): 
    first_term_a =  sf.Ei(-data/a_decay) 
    second_term_a = sf.Ei(-data*(1+ARE_A*a_decay)/a_decay)) 
    first_term_b =  sf.Ei(-data/b_decay)) 
    second_term_b = (sf.Ei(-data*(1+ARE_B*b_decay)/b_decay))
 return constant+on1_a_amplitude*(first_term_a-second_term_a+sf.log(1+ARE_A*a_decay)) + on1_b_amplitude*(first_term_b-second_term_b+sf.log(1+ARE_B*b_decay))

def mod2(data, a_decay, b_decay, constant, on2_a_amplitude, on2_b_amplitude, ARE_A, ARE_B): # not all parameters are used here
    first_term_a =  (sf.Ei((t2-t1-data)/a_decay)) - (sf.Ei((t2-data)/a_decay))
    second_term_a = (sf.Ei((t2-data)*(1+ARE_A*a_decay)/a_decay)) - (sf.Ei((t2-t1-data)*(1+ARE_A*a_decay)/a_decay))
    third_term_a = (sf.Ei((t2-data)/a_decay)) -(sf.Ei((t2-data)*(1+ARE_A*a_decay)/a_decay)) + sf.log(ARE_A*a_decay)

    first_term_b =  (sf.Ei((t2-t1-data)/b_decay)) - (sf.Ei((t2-data)/b_decay))
    second_term_b = (sf.Ei((t2-data)*(1+ARE_B*b_decay)/b_decay)) - (sf.Ei((t2-t1-data)*(1+ARE_B*b_decay)/a_decay))
    third_term_b = (sf.Ei((t2-data)/b_decay)) - (sf.Ei((t2-data)*(1+ARE_B*b_decay)/b_decay)) + sf.log(ARE_B*b_decay)

return constant+on2_a_amplitude*(sf.exp((t1-t2)/a_decay)*(first_term_a+second_term_a)+third_term_a)+on2_b_amplitude*(sf.exp((t1-t2)/b_decay)*(first_term_b+second_term_b)+third_term_b)  

x_1, x_3, y_1, y_3 = sf.variables('x_1, x_3, y_1, y_3')

a_dec = sf.Parameter(value=64.32, min=0.0, max=100)
b_dec = sf.Parameter(value=53.88, min=0.0, max=100)
con = sf.Parameter(value=2.911E-4, min=0.0, max=1)
off1_a_amp = sf.Parameter(value=0.015, min=0.0, max=1)
off1_b_amp = sf.Parameter(value=0.015, min=0.0, max=1)
off2_a_amp = sf.Parameter(value=0.015, min=0.0, max=1)
off2_b_amp = sf.Parameter(value=0.015, min=0.0, max=1)
AR_A = sf.Parameter(value=0.032, min=0.0, max=1)
AR_B = sf.Parameter(value=0.1, min=0.0, max=1)

model = sf.Model({
y_1: mod1(x_1, a_dec, b_dec, con, off1_a_amp, off1_b_amp, AR_A, AR_B),
y_3: mod2(x_3, a_dec, b_dec, con, off2_a_amp, off2_b_amp, AR_A, AR_B ),
})

fit = sf.Fit(model, x_1=x1, x_3=x3, y_1=y1, y_3=y3)
fit_result = fit.execute()
print(fit_result)

Regards, Ricardo

[pckroon]

The main issue is that we need a symbolic expression of your function, with (ideally) defined derivatives. We use sympy for all symbolic operations, so look there to find a symbolic equivalent function, or open an issue/feature request there. Otherwise you could try going for a https://symfit.readthedocs.io/en/stable/module_docs.html#symfit.core.models.CallableNumericalModel, but I don't know the required syntax from the top of my head.

[1000NightsAKnight]

Thanks for the reply @pckroon. As you advised, I did look at sympy to find a symbolic equivalent function. It tells me that Ei(z), as used in my my code is that function. However, I still get the following error: 'NameError: name 'Ei' is not defined'.

[pckroon]

You do need to import the Ei from sympy (not symfit) before using it of course.

[1000NightsAKnight]

I tried to follow your advice, @pckroon, by adding to my header the line: import sympy as sp

I then modified the original functions mod1 and mod2 to:

def mod1(data, a_decay, b_decay, constant, on1_a_amplitude, on1_b_amplitude, ARE_A, ARE_B): # not all parameters are used here
    first_term_a =  sp.Ei(-data/a_decay) 
    second_term_a = (sp.Ei(-data*(1+ARE_A*a_decay)/a_decay)) 
    first_term_b =  (sp.Ei(-data/b_decay)) 
    second_term_b = (sp.Ei(-data*(1+ARE_B*b_decay)/b_decay))
    return constant+on1_a_amplitude*(first_term_a-second_term_a+sf.log(1+ARE_A*a_decay)) + on1_b_amplitude*(first_term_b-second_term_b+sf.log(1+ARE_B*b_decay))

def mod2(data, a_decay, b_decay, constant, on2_a_amplitude, on2_b_amplitude, ARE_A, ARE_B): # not all parameters are used here
    first_term_a =  (sp.Ei((t2-t1-data)/a_decay)) - (sp.Ei((t2-data)/a_decay))
    second_term_a = (sp.Ei((t2-data)*(1+ARE_A*a_decay)/a_decay)) - (sp.Ei((t2-t1-data)*(1+ARE_A*a_decay)/a_decay))
    third_term_a = (sp.Ei((t2-data)/a_decay)) -(sp.Ei((t2-data)*(1+ARE_A*a_decay)/a_decay)) + sf.log(ARE_A*a_decay)

    first_term_b =  (sp.Ei((t2-t1-data)/b_decay)) - (sp.Ei((t2-data)/b_decay))
    second_term_b = (sp.Ei((t2-data)*(1+ARE_B*b_decay)/b_decay)) - (sp.Ei((t2-t1-data)*(1+ARE_B*b_decay)/a_decay))
    third_term_b = (sp.Ei((t2-data)/b_decay)) - (sp.Ei((t2-data)*(1+ARE_B*b_decay)/b_decay)) + sf.log(ARE_B*b_decay)

    return constant+on2_a_amplitude*(sf.exp((t1-t2)/a_decay)*(first_term_a+second_term_a)+third_term_a)+on2_b_amplitude*(sf.exp((t1-t2)/b_decay)*(first_term_b+second_term_b)+third_term_b)  

I still get the same 'NameError: name 'Ei' is not defined' error.

[pckroon]

I did a bit more digging. The issue is that Ei cannot be lambdified by sympy lambdify (which we use internally). Consider the following snippet:

from sympy import Symbol, Ei, lambdify
x = Symbol('x')
expr = Ei(x)
f = lambdify([x], expr, dummify=False)  
f(3)

This gives the NameError.

We could monkeypatch this in symfit.core.printing, but I'm really not enthusiastic about that. I'd much rather see this fixed in sympy.

[1000NightsAKnight]

Hello @pckroon. Thanks for continuing to dig. I believe I have encountered the same complaint against sympy here. Is there a work around or do I have to open a sympy issue/feature request?

[pckroon]

Yes, that's the same (or at least very similar) issue. Please do open an issue at sympy regarding Ei, since there seems to be a clear scipy implementation.

[1000NightsAKnight]

@pckroon I did open the issue at sympy and got the response that the Ei can be lambdified. Can you suggest what else the problem might be?

[1000NightsAKnight]

@pckroon My attempt to open an issue at sympy regarding Ei doesn't seem to have gone anywhere. Do we have further recourse?

[pckroon]

You got a reply on your sympy issue that you never responded to. Summary is that this is fixed in sympy==1.12

[1000NightsAKnight]

@pckroon I did indeed get a reply. But I am not sure what kind of response I could have given since my problem is with symfit. In other words: the problem of lambdifying Ei is fixed in sympy==1.12. Is it therefore fixed in symfit?

[pckroon]

Yes.

[1000NightsAKnight]

@pckroon Thanks! I will try it shortly.