Getting back the impedance expression or at least the rational polynomial equivalent expression for a complex data.

[richinex]

I am trying to back the impedance expression or at least the rational polynomial equivalent from which a complex data was generated. However, the final expression (best model) doesn't seem close to the actual. Although the fit is not that bad, it's not perfect either. I also noticed that the "s" variable was having a complex exponent which is not supposed to be. Therefore, I would like to constrain the exponent to always be real. The loss function used is the classic loss function used in complex nonlinear least squares in impedance spectroscopy. I would also like to know how I could further choose hyperparameters which might help. Below is my working example:

import numpy as np
import pysr
import matplotlib.pyplot as plt
from sympy import symbols

# Actual impedance expression
# Define the variables
R_s, C_dl, R_ct, C_1, R_1, s = symbols('R_s C_dl R_ct C_1 R_1 s')
# Define the expression
expr = R_s + 1/(1 / (R_ct + 1/(1/R_1 + C_1 * s)) + s * C_dl)


# Pade approximant expression
s, a_2, a_1, a_0, b_2, b_1 = symbols('s a_2 a_1 a_0 b_2 b_1')
# Define the expression
expr = (s**2 * a_2 + s * a_1 + a_0) / (s**2 * b_2 + s * b_1 + 1)



# Actual parameters and their values
# parameters = [R_s, C_dl, R_ct, C_1, R_1]
# parameter_values = np.array([5.26589219e+00, 7.46288724e-06, 8.27089860e-01,
#              1.99066599e+01, 3.40764484e-03, 2.19277541e+01])

# s  # basically the laplace variable (s) which is equivalent to sqrt(1j*2*pi*f) where f is the frequency array
# ydata # the impedance data

s = np.array([5.60499122e-01+5.60499122e-01j, 6.39067324e-01+6.39067324e-01j,
       7.28648858e-01+7.28648858e-01j, 8.31354548e-01+8.31354548e-01j,
       9.47890025e-01+9.47890025e-01j, 1.07959755e+00+1.07959755e+00j,
       1.23182355e+00+1.23182355e+00j, 1.40348801e+00+1.40348801e+00j,
       1.60110574e+00+1.60110574e+00j, 1.82485293e+00+1.82485293e+00j,
       2.08216182e+00+2.08216182e+00j, 2.37138163e+00+2.37138163e+00j,
       2.70553338e+00+2.70553338e+00j, 3.08529184e+00+3.08529184e+00j,
       3.51822044e+00+3.51822044e+00j, 4.01060524e+00+4.01060524e+00j,
       4.57073201e+00+4.57073201e+00j, 5.21294317e+00+5.21294317e+00j,
       5.93176514e+00+5.93176514e+00j, 6.77253665e+00+6.77253665e+00j,
       7.72594722e+00+7.72594722e+00j, 8.80893515e+00+8.80893515e+00j,
       1.00421673e+01+1.00421673e+01j, 1.14594316e+01+1.14594316e+01j,
       1.30609525e+01+1.30609525e+01j, 1.48928319e+01+1.48928319e+01j,
       1.69730388e+01+1.69730388e+01j, 1.93351888e+01+1.93351888e+01j,
       2.20668725e+01+2.20668725e+01j, 2.51913024e+01+2.51913024e+01j,
       2.86896719e+01+2.86896719e+01j, 3.27304613e+01+3.27304613e+01j,
       3.73058380e+01+3.73058380e+01j, 4.25388924e+01+4.25388924e+01j,
       4.84758837e+01+4.84758837e+01j, 5.52880607e+01+5.52880607e+01j,
       6.29158704e+01+6.29158704e+01j, 7.17789102e+01+7.17789102e+01j,
       8.19939527e+01+8.19939527e+01j, 9.34538794e+01+9.34538794e+01j,
       1.06494833e+02+1.06494833e+02j, 1.21513314e+02+1.21513314e+02j,
       1.38433071e+02+1.38433071e+02j, 1.57937474e+02+1.57937474e+02j,
       1.79884419e+02+1.79884419e+02j, 2.05176367e+02+2.05176367e+02j,
       2.33802721e+02+2.33802721e+02j, 2.67047099e+02+2.67047099e+02j,
       3.04428946e+02+3.04428946e+02j, 3.46876057e+02+3.46876057e+02j,
       3.95539270e+02+3.95539270e+02j, 4.50844812e+02+4.50844812e+02j,
       5.14317121e+02+5.14317121e+02j, 5.85178263e+02+5.85178263e+02j])

ydata = np.array([45.1  -1.09j, 47.5  -1.43j, 46.8  -1.77j, 46.2  -2.29j,
       46.2  -2.97j, 47.2  -3.8j , 47.   -4.85j, 45.1  -5.99j,
       45.8  -7.33j, 42.3  -9.05j, 42.6 -10.2j , 36.5 -10.8j ,
       34.5 -11.2j , 32.1 -10.2j , 30.   -9.18j, 29.4  -8.j  ,
       27.3  -6.64j, 26.7  -5.18j, 25.3  -4.12j, 25.4  -3.26j,
       25.2  -2.51j, 24.9  -1.94j, 24.9  -1.64j, 25.4  -1.35j,
       25.5  -1.24j, 24.8  -1.1j , 24.7  -1.03j, 23.9  -1.04j,
       25.2  -1.1j , 24.9  -1.27j, 25.   -1.46j, 25.4  -1.65j,
       24.4  -1.98j, 24.5  -2.34j, 24.5  -2.91j, 23.8  -3.47j,
       22.9  -4.13j, 22.3  -4.91j, 20.9  -5.66j, 20.3  -6.03j,
       18.4  -6.96j, 17.6  -7.24j, 16.5  -7.74j, 14.3  -7.42j,
       12.7  -7.17j, 11.2  -6.76j,  9.85 -5.89j,  8.68 -5.38j,
        7.92 -4.53j,  7.2  -3.83j,  6.81 -3.2j ,  6.65 -2.67j,
        6.11 -2.16j,  5.86 -1.77j], dtype=complex)

model = pysr.PySRRegressor(
    binary_operators=["+", "-", "*", "/", "^"], niterations=100, populations=100, 
    loss="f(x, y) = abs(((((x.re) - (y.re))^2 + ((x.im) - (y.im))^2)/abs(x)^2))", 
)
model.fit(s[:, None], ydata)

The expression found is shown below (where j represents the imaginary unit)
expr = (((((1.618935e0 + 0.96540135e0j) + s) ^ (-0.112029515e0 - 0.0105559165e0j)) / (0.020140855e0 - 0.0013685775e0j)) + ((s ^ (0.033040527e0 - 0.73647106e0j)) - ((-1.6239787e0 - 0.12344298e0j) / ((x0 ^ (0.033040527e0 - 0.73647106e0j)) ^ (0.75328916e0 - 0.108941875e0j)))))


# Zpred = model.predict(s[:, None])
# plt.plot(ydata.real, -ydata.imag, 'bo', label='Data')
# plt.plot(Zpred.real, -Zpred.imag, 'r-', label='Predicted') # nit bad but not perfect

[MilesCranmer]

If your expression is polynomials only up to degree 2, then you are likely better off avoiding ^ and just relying on * or perhaps square.

If you still need the power law, maybe you could define a custom operator? For example

binary_operators=[..., "mypow(x,y)=x^(y.re)"],
extra_sympy_mappings ={"mypow": lambda x,y: x**re(y)} # give sympy equivalent

If you want only integers you could further use

mypow(x, y) = x^ceil(Int, y.re)

for further tuning advice I would check out https://astroautomata.com/PySR/tuning/

[MilesCranmer]

Also - for your loss function, this is a simpler equivalent:

f(x, y) = abs2(x - y)/abs2(x)

Note that the form of this is f(prediction, target). Perhaps you meant to divide by abs2(y) instead?

[richinex]

Ah yes , I meant to divide by the absolute value of the target. Thanks. I'll also try the other suggestions for the real exponent.