Fitting custom model resulted using sympy

[karimzadehz]

Hi Symfit, I'm trying to implement a custom model using Fit . I read the doc to see how to do that. Unfortunately, there are some errors in my code and it seems to me that might be related to my model that is driven from sympy. Here, I try to show very simple one of my problem:

import numpy as np
import sympy as sp
from symfit import Fit
from symfit import parameters, variables, Model

L_vals = np.array([0.299, 0.295, 0.290, 0.284, 0.279, 0.273, 0.268, 0.262, 0.256, 0.250])
K_vals = np.array([2.954, 3.056, 3.119, 3.163, 3.215, 3.274, 3.351, 3.410, 3.446, 3.416])
VA_vals = np.array([0.919, 0.727, 0.928, 0.629, 0.656, 0.854, 0.955, 0.981, 0.908, 0.794])

gamma, alpha, beta, eta, L, K, VA = sp.symbols('gamma, alpha, beta, eta, L, K, VA')
y = gamma - (1 / eta) * sp.log(alpha * L ** -eta + beta * K ** -eta)

gamma, alpha, beta, eta = parameters('gamma, alpha, beta, eta')
L, K, VA = variables('L, K, VA')
model = Model({VA: y})
fit = Fit(model, L=L_vals, K=K_vals, VA=VA_vals)
fit_result = fit.execute()

Traceback (most recent call last):
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 286, in __get__
    return getattr(obj, self.cache_attr)
AttributeError: 'Model' object has no attribute '_cached_connectivity_mapping'. Did you mean: 'connectivity_mapping'?

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "F:\Zohreh\MainZohreh\postdoc-field\CSU\pythonProject\symfit_ask.py", line 18, in <module>
    model = Model({VA: f_Vi_est })
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 916, in __init__
    super(HessianModel, self).__init__(*args, **kwargs)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 865, in __init__
    super(GradientModel, self).__init__(*args, **kwargs)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 143, in __init__
    self._init_from_dict(model)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 692, in _init_from_dict
    super(BaseCallableModel, self)._init_from_dict(model_dict)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 327, in _init_from_dict
    ordered = list(toposort(self.connectivity_mapping))
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 289, in __get__
    setattr(obj, self.cache_attr, self.fget(obj))
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 401, in connectivity_mapping
    vars, params = seperate_symbols(expr)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 86, in seperate_symbols
    for symbol in func.free_symbols:
AttributeError: 'function' object has no attribute 'free_symbols'

Please kindly let me know how deal with it. Regards, ZK

[pckroon]

I expect the issue is in here:

gamma, alpha, beta, eta, L, K, VA = sp.symbols('gamma, alpha, beta, eta, L, K, VA')
y = gamma - (1 / eta) * sp.log(alpha * L ** -eta + beta * K ** -eta)

gamma, alpha, beta, eta = parameters('gamma, alpha, beta, eta')
L, K, VA = variables('L, K, VA')
model = Model({VA: y})

The model you made (y) is constructed with sympy symbols rather than symfit parameters/variables. Try the following instead:

gamma, alpha, beta, eta = parameters('gamma, alpha, beta, eta')
L, K, VA = variables('L, K, VA')
y = gamma - (1 / eta) * sp.log(alpha * L ** -eta + beta * K ** -eta)
model = Model({VA: y})

[karimzadehz]

    import numpy as np
    import sympy as sp
    from symfit import Fit
    from symfit import parameters, variables, Model

    L_vals = np.array([0.299, 0.295, 0.290, 0.284, 0.279, 0.273, 0.268, 0.262, 0.256, 0.250])
    K_vals = np.array([2.954, 3.056, 3.119, 3.163, 3.215, 3.274, 3.351, 3.410, 3.446, 3.416])
    VA_vals = np.array([0.919, 0.727, 0.928, 0.629, 0.656, 0.854, 0.955, 0.981, 0.908, 0.794])

    gamma, alpha, beta, eta = parameters('gamma, alpha, beta, eta')
    L, K, VA = variables('L, K, VA')
    y = gamma - (1 / eta) * sp.log(alpha * L ** -eta + beta * K ** -eta)

    model = Model({VA: y})

    fit = Fit(model, L=L_vals, K=K_vals, VA=VA_vals)
    fit_result = fit.execute()
    print(fit_result)

still gives:

Traceback (most recent call last):
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 286, in __get__
    return getattr(obj, self.cache_attr)
AttributeError: 'Model' object has no attribute '_cached_numerical_components'. Did you mean: 'numerical_components'?

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "F:\Zohreh\MainZohreh\postdoc-field\CSU\pythonProject\symfit_ask.py", line 17, in <module>
    fit_result = fit.execute()
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\fit.py", line 584, in execute
    minimizer_ans = self.minimizer.execute(**minimize_options)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 427, in wrapped_func
    return func(*bound_args.args, **bound_args.kwargs)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\minimizers.py", line 413, in execute
    return super(ScipyGradientMinimize, self).execute(jacobian=jacobian, **minimize_options)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 427, in wrapped_func
    return func(*bound_args.args, **bound_args.kwargs)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\minimizers.py", line 350, in execute
    ans = minimize(
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_minimize.py", line 687, in minimize
    res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_optimize.py", line 1296, in _minimize_bfgs
    sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_optimize.py", line 263, in _prepare_scalar_function
    sf = ScalarFunction(fun, x0, args, grad, hess,
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_differentiable_functions.py", line 158, in __init__
    self._update_fun()
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_differentiable_functions.py", line 251, in _update_fun
    self._update_fun_impl()
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_differentiable_functions.py", line 155, in update_fun
    self.f = fun_wrapped(self.x)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\scipy\optimize\_differentiable_functions.py", line 137, in fun_wrapped
    fx = fun(np.copy(x), *args)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 427, in wrapped_func
    return func(*bound_args.args, **bound_args.kwargs)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\objectives.py", line 312, in __call__
    evaluated_func = super(LeastSquares, self).__call__(
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\objectives.py", line 93, in __call__
    result = self.model(**key2str(parameters))._asdict()
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 708, in __call__
    return ModelOutput(self.keys(), self.eval_components(*args, **kwargs))
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 649, in eval_components
    components = dict(zip(self, self.numerical_components))
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 289, in __get__
    setattr(obj, self.cache_attr, self.fget(obj))
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 856, in numerical_components
    components.append(sympy_to_py(expr, ordered))
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 127, in sympy_to_py
    lambdafunc = lambdify(args, func, dummify=False)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\utilities\lambdify.py", line 875, in lambdify
    funcstr = funcprinter.doprint(funcname, iterable_args, _expr, cses=cses)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\utilities\lambdify.py", line 1126, in doprint
    argstrs, expr = self._preprocess(args, expr)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\utilities\lambdify.py", line 1195, in _preprocess
    s = self._argrepr(arg)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\printing\codeprinter.py", line 150, in doprint
    lines = self._print(expr).splitlines()
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\printing\printer.py", line 309, in _print
    return getattr(expr, self.printmethod)(self, **kwargs)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\argument.py", line 76, in _sympystr
    return printer.doprint(self.name)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\printing\codeprinter.py", line 143, in doprint
    expr = self._handle_UnevaluatedExpr(expr)
  File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\sympy\printing\codeprinter.py", line 103, in _handle_UnevaluatedExpr
    return expr.replace(re, lambda arg: arg if isinstance(
TypeError: Basic.replace() missing 1 required positional argument: 'value'

[pckroon]

Interesting. print(model(L=0.225, K=3, alpha=1, beta=1, gamma=1, eta=1)) gets the same traceback on Python 3.9 on the master branch. On the other hand model2 = Model({VA: K + L * alpha}); print(model2(K=np.arange(5), L=np.arange(1), alpha=1)) DOES work.

I'll dig a bit further, but I'm also rather swamped with work at the moment.

[pckroon]

model2 = Model({VA: L**alpha + K}); print(model2(L=1, alpha=1, K=1)) Works model2 = Model({VA: L**alpha + K**beta}); print(model2(L=1, alpha=1, K=1, beta=1)) Breaks

... and print(sp.lambdify((K, L, alpha, beta), L**alpha + K**beta)) Breaks.

So there's a (compatibility?) issue with sympy printing. Could you try with an older sympy version? I really don't have the time at the moment.

[karimzadehz]

I gave it a try:

`import numpy as np
import sympy as sp
from symfit import Fit
from symfit import parameters, variables, Model

L_vals = np.array([0.299, 0.295, 0.290, 0.284, 0.279, 0.273, 0.268, 0.262, 0.256, 0.250])
K_vals = np.array([2.954, 3.056, 3.119, 3.163, 3.215, 3.274, 3.351, 3.410, 3.446, 3.416])
VA_vals = np.array([0.919, 0.727, 0.928, 0.629, 0.656, 0.854, 0.955, 0.981, 0.908, 0.794])

gamma, alpha, beta = parameters('gamma, alpha, beta')
L, K, VA = variables('L, K, VA')
y = gamma * L - sp.log(alpha) + beta + K ** 0.5
from sympy.printing.numpy import NumPyPrinter
code = NumPyPrinter().doprint(y)
model = Model({VA: code})

fit = Fit(model, L=L_vals, K=K_vals, VA=VA_vals)
fit_result = fit.execute()
print(fit_result)`

no effect!

[pckroon]

With which sympy version?

[karimzadehz]

Name: sympy Version: 1.10.1

[karimzadehz]

similar issue have been solved for scipy based on the following:

`import numpy as np
import sympy as sp
from scipy.optimize import minimize, curve_fit

L_vals = np.array([0.299, 0.295, 0.290, 0.284, 0.279, 0.273, 0.268, 0.262, 0.256, 0.250])
K_vals = np.array([2.954, 3.056, 3.119, 3.163, 3.215, 3.274, 3.351, 3.410, 3.446, 3.416])
VA_vals = np.array([0.919, 0.727, 0.928, 0.629, 0.656, 0.854, 0.955, 0.981, 0.908, 0.794])

gamma, alpha, beta, eta, L, K, VA = sp.symbols('gamma, alpha, beta, eta, L, K, VA')

Vi_est = gamma - (1 / eta) * sp.log(alpha * L ** -eta + beta * K ** -eta)

# Note that the first entry in the args tuple is a 4 tuple
# That is needed to unpack the arguments from an array
f_Vi_est = sp.lambdify(((gamma, alpha, beta, eta), L, K, VA), Vi_est)

def f(param):
    Vi_est_vals = f_Vi_est(param, L_vals, K_vals, VA_vals)
    return np.sum((np.log(VA_vals) - Vi_est_vals) ** 2)

bnds = [(1, np.inf), (0, 1), (0, 1), (-1, np.inf)]
x0 = (1, 0.01, 0.98, 1)

result = minimize(f, x0, bounds=bnds)

print(result.message)
print(result.x[0], result.x[1], result.x[2], result.x[3])`

[pckroon]

Could you try the original issue with sympy 1.9 or even 1.8?

[pckroon]

pip install sympy==1.8, or something along those lines.

[karimzadehz]

what I get by 1.8 version Traceback (most recent call last): File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 286, in get return getattr(obj, self.cache_attr) AttributeError: 'Model' object has no attribute '_cached_connectivity_mapping'. Did you mean: 'connectivity_mapping'?

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "F:\Zohreh\MainZohreh\postdoc-field\CSU\pythonProject\symfit_ask.py", line 17, in model = Model({VA: code}) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 916, in init super(HessianModel, self).init(*args, *kwargs) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 865, in init super(GradientModel, self).init(args, **kwargs) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 143, in init self._init_from_dict(model) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 692, in _init_from_dict super(BaseCallableModel, self)._init_from_dict(model_dict) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 327, in _init_from_dict ordered = list(toposort(self.connectivity_mapping)) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 289, in get setattr(obj, self.cache_attr, self.fget(obj)) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\models.py", line 401, in connectivity_mapping vars, params = seperate_symbols(expr) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 86, in seperate_symbols for symbol in func.free_symbols: AttributeError: 'str' object has no attribute 'free_symbols'

[karimzadehz]

good news by sympy 1.11 I got:

:2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) :2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) :2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) :2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) :2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) :2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) :2: RuntimeWarning: invalid value encountered in log return (beta + gamma*L + sqrt(K) - log(alpha)) Parameter Value Standard Deviation alpha -2.389480e+03 nan beta -3.726690e+03 nan gamma -1.025736e+03 nan Status message Desired error not necessarily achieved due to precision loss. Number of iterations 2 Objective Minimizer Goodness of fit qualifiers: chi_squared nan objective_value nan r_squared nan

[karimzadehz]

and with.....exp instaed of log in y: y = gamma * L - sp.exp(alpha) + beta + K ** 0.5 here is what I got: RuntimeWarning: invalid value encountered in sqrt return np.sqrt(self.variance(param))

Parameter Value Standard Deviation alpha 7.937665e-01 nan beta 9.858380e-01 nan gamma 9.485058e-01 2.767540e+00 Status message Desired error not necessarily achieved due to precision loss. Number of iterations 19 Objective <symfit.core.objectives.LeastSquares object at 0x000002ACD7F11420> Minimizer <symfit.core.minimizers.BFGS object at 0x000002ACD7F113F0>

Goodness of fit qualifiers: chi_squared 0.13416606704535364 objective_value 0.06708303352267682 r_squared 0.07428839797345088

[karimzadehz]

while in my real code still it gives error, data.xlsx

`import numpy
from numpy import pi
import sympy as sp
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import itertools

from sympy import lambdify

anion = {'Cl': 1}
cation = {'Na': 1}
neutral = {'water': 0}
alfa = 13
alfa1 = 0.0
t = 298.15
m1 = 18.02
ro = 14.0405
nu_c = 1.0
nu_a = 1.0

xa = {'XCl': []}
xc = {'XNa': []}
xn = {'Xwater': []}
ionic_dic = dict()
for d in (cation, anion):
    ionic_dic.update(d)

L=['BNaCl', 'B1NaCl', 'WwaterNaCl', 'UwaterNaCl', 'VwaterNaCl', 'XCl', 'XNa', 'Xwater']

for s in range(len(L)):
    locals()[L[s]] = sp.symbols(L[s])
# z = np.concatenate((zc, za))
# x = np.concatenate((xc, xa))
# ec = np.zeros((len(cation)), dtype=np.float64)
# ea = np.zeros((len(anion)), dtype=np.float64)

def ea(x, z):
    xzc = 0
    for k1, k2 in zip(xa, anion):
        xzc += sp.symbols(k1) * anion[k2]
    return x * z / xzc

def ec(x, z):
    xzc = 0
    for k1, k2 in zip(xc, cation):
        xzc += sp.symbols(k1) * cation[k2]
    return x * z / xzc

def f():
    return 1 / (0.5 * (sum(sp.symbols(f'X{str(j)}') * cation[j] for j in cation) + \
                       sum(sp.symbols(f'X{str(k)}') * anion[k] for k in anion)))

'''Here sqrt numpy is used'''

def ionicstr():
    Ionic = 0.0
    sqrt_I =0.0
    for k1 in ionic_dic:
        Ionic += 0.5 * (sp.symbols(f'X{str(k1)}') * (ionic_dic[k1] ** 2))
    sqrt_I += sp.sqrt(Ionic)
    return Ionic, sqrt_I

# def dieletricconstant_water(t):
#     # for TK: 273-372
#     return 0.24921e3 - 0.79069 * t + 0.72997e-3 * t ** 2
#
#
# # @numba.njit
# def density_water(t):
#     # for t: 273-372
#     return (
#             0.183652
#             + 0.00724987 * t
#             - 0.203449e-4 * t ** 2
#             + 1.73702e-8 * t ** 3
#     )

# print(density_water(298.15))

def afix(t, m1):
    ds = 997.00  # kg/m³, water density
    e = 1.60218 * (10 ** (-19))  # coulombs, electron charge
    na = 6.023 * (10 ** 23)  # Nav = 6.023 x 10**23 atom/gm.mole
    ep0 = 8.85419 * (10 ** (-12))

    kb = 1.38064852 * (10 ** (-23))  # 1.380649×10−23 J K-1, 1.380649 × 10-23 m2 kg s-2 K-1
    epr = 78.41
    # dieletric = 305.7*sp.exp(-sp.exp(-12.741+0.01875*t)-(t/219))
    # dieletric = 0.24921e3 - 0.79069 * t + 0.72997e-3 * t ** 2
    # ds = -58.05 - 0.01098*t + 3053/t - 1.75e5/t**2 + 21.8*sp.log10(t)
    a = (1.0 / 3.0) * ((2.0 * pi * na * ds) ** 0.5) * (((e * e) / (4.0 * pi * ep0 * epr * kb * t)) ** (3.0 / 2.0))
    # dieletric = 78.2
    # a = 1400684 * (ds / dieletric * t) ** (1.5)
    '''ax should be 2.917 but now is 2.915170'''
    Ax = ((1000 / m1) ** 0.5) * a
    #print(a, Ax)
    return Ax

# print(afix(298.15, 18.02))

def cationshortrange():
    term_wnMa = 0
    for k in anion:
        for j in cation:
            if j == next(iter(cation)):
                for i in neutral:
                    term_wnMa += sp.symbols(f'X{str(i)}') * ea(sp.symbols(f'X{str(k)}'), anion[k]) * \
                                 ((cation[j] + anion[k]) / (anion[k])) * sp.symbols(f'W{str(i)}{str(j)}{str(k)}')
    zmterm = 0
    for j in cation:
        if j == next(iter(cation)):
            zmterm += (cation[j] / 2) + (1 / f())

    term_w_solvent1 = 0
    for k in anion:
        for j in cation:
            if j == next(iter(cation)):
                for i in neutral:
                    if i == next(iter(neutral)):
                        term_w_solvent1 += - ea(sp.symbols(f'X{str(k)}'), anion[k]) * ((1 - ec(sp.symbols(f'X{str(j)}'),\
                                            cation[j]) / 2) * (cation[j] + anion[k]) / anion[k])\
                                           * sp.symbols(f'W{str(i)}{str(j)}{str(k)}')
    zmm = 0
    for j in cation:
        if j == next(iter(cation)):
            zmm += cation[j] / 2

    term_w_solvent2 = 0
    for k in anion:
        for j in cation:
            if j != next(iter(cation)):
                for i in neutral:
                    if i == next(iter(neutral)):
                        term_w_solvent2 += - ea(sp.symbols(f'X{str(k)}'), anion[k]) * (- zmm) * \
                                           ec(sp.symbols(f'X{str(j)}'), cation[j])\
                                           * ((cation[j] + anion[k]) / (cation[j] * anion[k])) * \
                                           sp.symbols(f'W{str(i)}{str(j)}{str(k)}')

    term_wnca = 0
    for k in anion:
        for j in cation:
            for i in neutral:
                term_wnca += term_wnMa - sp.symbols(f'X{str(i)}') * ea(sp.symbols(f'X{str(k)}'), anion[k]) * zmterm * \
                             ec(sp.symbols(f'X{str(j)}'), cation[j]) * \
                            sp.symbols(f'W{str(i)}{str(j)}{str(k)}') * ((cation[j] + anion[k]) / (cation[j] * anion[k]))
    term_unMa = 0
    for k in anion:
        for j in cation:
            if j == next(iter(cation)):
                for i in neutral:
                    term_unMa += sp.symbols(f'X{str(i)}') * sp.symbols(f'X{str(k)}') * \
                                 ((cation[j] + anion[k]) ** 2 / (cation[j] * anion[k])) * \
                                 sp.symbols(f'U{str(i)}{str(j)}{str(k)}')

    term_unca = 0
    for k in anion:
        for j in cation:
            for i in neutral:
                term_unca += term_unMa - 2 * sp.symbols(f'X{str(j)}') * ((cation[j] + \
                            anion[k]) ** 2 / (cation[j] * anion[k])) * sp.symbols(f'U{str(i)}{str(j)}{str(k)}')
    term_vnMa = 0
    for k in anion:
        for j in cation:
            if j == next(iter(cation)):
                for i in neutral:
                    term_vnMa += 4 * sp.symbols(f'X{str(i)}') ** 2 * sp.symbols(f'X{str(k)}') * \
                                 sp.symbols(f'V{str(i)}{str(j)}{str(k)}')
    term_vnca = 0
    for k in anion:
        for j in cation:
            for i in neutral:
                term_vnca += term_vnMa - 12 * sp.symbols(f'X{str(j)}') * sp.symbols(f'X{str(i)}') ** 2 \
                             * sp.symbols(f'X{str(k)}') * sp.symbols(f'V{str(i)}{str(j)}{str(k)}')
    lnfMs = term_wnca + term_unca + term_vnca + term_w_solvent1 + term_w_solvent2
    return lnfMs

def anionshortrange():
    termx_w_solvent1 = 0
    for k in anion:
        if k == next(iter(anion)):
            for j in cation:
                for i in neutral:
                    if i == next(iter(neutral)):
                        termx_w_solvent1 += - ec(sp.symbols(f'X{str(j)}'), cation[j]) * ((1 - ea(sp.symbols(f'X{str(k)}'), anion[k]) / 2) *\
                                          (cation[j] + anion[k]) / cation[j]) * sp.symbols(f'W{str(i)}{str(j)}{str(k)}')
    zxx = 0
    for k in anion:
        if k == next(iter(anion)):
            zxx += anion[k] / 2

    termx_w_solvent2 = 0
    for k in range(len(anion)):
        if k != next(iter(anion)):
            for j in range(len(cation)):
                for i in range(len(neutral)):
                    if i == next(iter(neutral)):
                        termx_w_solvent2 += -ec(sp.symbols(f'X{str(j)}'), cation[j]) * (- zxx) * \
                                            ea(sp.symbols(f'X{str(k)}'), anion[k]) * ((cation[j] + anion[k]) / \
                                               (cation[j] * anion[k])) * sp.symbols(f'W{str(i)}{str(j)}{str(k)}')
    term_wncX = 0
    for k in anion:
        if k == next(iter(anion)):
            for j in cation:
                for i in neutral:
                    term_wncX += sp.symbols(f'X{str(i)}') * ec(sp.symbols(f'X{str(j)}'), cation[j]) * \
                                 ((cation[j] + anion[k]) / (cation[j])) * \
                                 sp.symbols(f'W{str(i)}{str(j)}{str(k)}')
    zxterm = 0
    for k in anion:
        if k == next(iter(anion)):
            zxterm += (anion[k] / 2) + (1 / f())
    term_wnca = 0
    for k in anion:
        for j in cation:
            for i in neutral:
                term_wnca += term_wncX - sp.symbols(f'X{str(i)}') * ec(sp.symbols(f'X{str(j)}'), cation[j]) * zxterm * \
                             ea(sp.symbols(f'X{str(k)}'), anion[k]) * \
                             ((cation[j] + anion[k]) / cation[j]) * sp.symbols(f'W{str(i)}{str(j)}{str(k)}')

    term_uncX = 0
    for k in anion:
        if k == next(iter(anion)):
            for j in cation:
                for i in neutral:
                    term_uncX += sp.symbols(f'X{str(i)}') * sp.symbols(f'X{str(j)}') * \
                                 ((cation[j] + anion[k]) ** 2 / (cation[j] * anion[k])) * \
                                 sp.symbols(f'U{str(i)}{str(j)}{str(k)}')

    term_unca = 0
    for k in anion:
        for j in cation:
            for i in neutral:
                term_unca += term_uncX - 2 * sp.symbols(f'X{str(i)}') * sp.symbols(f'X{str(j)}') * \
                             sp.symbols(f'X{str(k)}') * ((cation[j] + anion[k]) ** 2 / (cation[j] * anion[k])) \
                             * sp.symbols(f'U{str(i)}{str(j)}{str(k)}')

    term_vncX = 0
    for k in anion:
        if k == next(iter(anion)):
            for j in cation:
                for i in neutral:
                    term_vncX += 4 * sp.symbols(f'X{str(i)}') ** 2 * sp.symbols(f'X{str(j)}') * \
                                 sp.symbols(f'V{str(i)}{str(j)}{str(k)}')

    term_vnca = 0
    for k in anion:
        for j in cation:
            for i in neutral:
                term_vnca += term_vncX - 4 * 3 * sp.symbols(f'X{str(k)}') * sp.symbols(f'X{str(i)}') ** 2 \
                             * sp.symbols(f'V{str(i)}{str(j)}{str(k)}') * sp.symbols(f'X{str(j)}')
    lnfXs = term_wnca + term_unca + term_vnca + termx_w_solvent1 + termx_w_solvent2
    return lnfXs

'''Here exp numpy is used'''

def g(x):
    return 2 * (1 - (1 + x) * sp.exp(-x)) / x ** 2

def cationlongrange():
    #global term1_c_l,  term2_c_l,  term3_c_l,  term4_c_l,  term5_c_l
    for j in cation:
        if j == next(iter(cation)):
            '''Here log numpy is used'''
            term1_c_l = -cation[j] ** 2 * afix(t, m1) * ((2 / ro) * sp.log(1 + ro * ionicstr()[1]) + ionicstr()[1] * (
                    1 - 2 * (ionicstr()[0]) / cation[j] ** 2) / (1 + ro * ionicstr()[1]))
    term2_c_l = 0.0
    for k in anion:
        for j in cation:
            if j == next(iter(cation)):
                term2_c_l += sp.symbols(f'X{str(k)}') * sp.symbols(f'B{str(j)}{str(k)}') * \
                             g(alfa * ionicstr()[1])

    term3_c_l = 0.0
    for j in cation:
        if j == next(iter(cation)):
            '''Here log numpy is used'''
            term3_1 = ((cation[j] ** 2 * (g(alfa * ionicstr()[1])) / (2 * ionicstr()[0])) + (
                    (1 - cation[j] ** 2) / (2 * ionicstr()[0])) * sp.exp(-alfa * ionicstr()[1]))
    for k in anion:
        for j in cation:
            term3_c_l += -sp.symbols(f'X{str(j)}') * sp.symbols(f'X{str(k)}') * \
                         sp.symbols(f'B{str(j)}{str(k)}') * term3_1

    term4_c_l = 0.0
    for j in cation:
        if j == next(iter(cation)) and cation[j] > 1:
            for k in anion:
                term4_c_l += sp.symbols(f'X{str(k)}') * sp.symbols(f'B1{str(j)}{str(k)}') * g(alfa1 * ionicstr()[1])

    term5_c_l = 0.0
    for j in cation:
        if j == next(iter(cation)) and cation[j] > 1:
            '''Here log numpy is used'''
            term5_1 = (cation[j] ** 2 * (g(alfa1 * ionicstr()[1]) / (2 * ionicstr()[0])) + ((1 - cation[j] ** 2) /\
                      (2 * ionicstr()[0])) * sp.exp(-alfa1 * ionicstr()[1]))
    for k in anion:
        for j in cation:
            if cation[j] > 1:
                term5_c_l += -sp.symbols(f'X{str(j)}') * sp.symbols(f'X{str(k)}') * \
                             sp.symbols(f'B1{str(j)}{str(k)}') * term5_1
    lnfMDH = term1_c_l + term2_c_l + term3_c_l + term4_c_l + term5_c_l
    return lnfMDH

def anionlongrange():
    #global term1_a_l, term2_a_l, term3_a_l, term4_a_l, term5_a_l
    for k in anion:
        if k == next(iter(anion)):
            '''Here log numpy is used'''
            term1_a_l = -anion[k] ** 2 * afix(t, m1) * ((2 / ro) * sp.log(1 + ro * ionicstr()[1]) + ionicstr()[1] * ( \
                       1 - 2 * (ionicstr()[0]) / anion[k] ** 2) / (1 + ro * ionicstr()[1]))

    term2_a_l = 0.0
    for k in anion:
        if k == next(iter(anion)):
            for j in cation:
                term2_a_l += sp.symbols(f'X{str(j)}') * sp.symbols(f'B{str(j)}{str(k)}') * \
                             g(alfa * ionicstr()[1])

    term3_a_l = 0.0
    for k in anion:
        if k == next(iter(anion)):
            '''Here log numpy is used'''
            term3_1 = ((anion[k] ** 2 * (g(alfa * ionicstr()[1])) / (2 * ionicstr()[0])) + (
                    (1 - anion[k] ** 2) / (2 * ionicstr()[0])) * sp.exp(-alfa * ionicstr()[1]))
    for k in anion:
        for j in cation:
            term3_a_l += -sp.symbols(f'X{str(j)}') * sp.symbols(f'X{str(k)}') * \
                         sp.symbols(f'B{str(j)}{str(k)}') * term3_1

    term4_a_l = 0.0
    for k in anion:
        if k == next(iter(anion)):
            for j in cation:
                if cation[j] > 1:
                    term4_a_l += sp.symbols(f'X{str(j)}') * sp.symbols(f'B1{str(j)}{str(k)}') * \
                                g(alfa1 * ionicstr()[1])

    term5_a_l = 0.0
    for k in anion:
        if k == next(iter(anion)):
            '''Here log numpy is used'''
            term5_1 = (anion[k] ** 2 * (g(alfa1 * ionicstr()[1]) / (2 * ionicstr()[0])) + (
                    (1 - anion[k] ** 2) / (2 * ionicstr()[0])) * sp.exp(-alfa1 * ionicstr()[1]))
    for k in anion:
        for j in cation:
            if cation[j] > 1:
                term5_a_l += -sp.symbols(f'X{str(j)}') * sp.symbols(f'X{str(k)}') * \
                            sp.symbols(f'B1{str(j)}{str(k)}') * term5_1
    lnfXDH = term1_a_l + term2_a_l + term3_a_l + term4_a_l + term5_a_l
    return lnfXDH

def lnfca():
    lnFca = (1 / (nu_c + nu_a)) * (anionlongrange() + cationshortrange() + cationlongrange() + anionshortrange())
    return lnFca

from sympy.printing.numpy import NumPyPrinter

expr = lnfca()
code = NumPyPrinter().doprint(expr)
from scipy.optimize import minimize, curve_fit
import pandas as pd
#exp_NaCl path: F:\Zohreh\MainZohreh\postdoc-field\CSU\Duplicat_Pure
df = pd.read_excel(r'F:\Zohreh\MainZohreh\postdoc-field\CSU\Duplicat_Pure\data.xlsx', sheet_name='NaCl_exp')
XNa_vals = df['XNa'].to_numpy()
XCl_vals = df['XCl'].to_numpy()
Xwater_vals = df['Xwater'].to_numpy()
ln_gama_x_vals = df['ln_gama_x'].to_numpy()

#WwaterNaCl, UwaterNaCl, VwaterNaCl, BNaCl, XCl, XNa, Xwater, gama_x = sp.symbols('WwaterNaCl, UwaterNaCl, VwaterNaCl, BNaCl, XCl, XNa, Xwater,  gama_x')
y = lnfca()
# ln_gama_x_est = code
from symfit import Fit
from symfit import parameters, variables, Model
# BNaCl = Parameter('BNaCl')
# UwaterNaCl = Parameter('UwaterNaCl')
# VwaterNaCl = Parameter('VwaterNaCl')
# WwaterNaCl = Parameter('WwaterNaCl')
# XCl = Variable('XCl')
# XNa = Variable('XNa')
# Xwater = Variable('Xwater')
# y = Variable('y')
# model_dict = Model({y: lnfca()})
# fit = Fit(model_dict, XCl=XCl_vals, XNa=XNa_vals, Xwater=Xwater_vals, y=ln_gama_x_vals)
# fit_result = fit.execute()
# print(fit_result)
#****************************************
BNaCl, WwaterNaCl, UwaterNaCl, VwaterNaCl = parameters('BNaCl,WwaterNaCl, UwaterNaCl, VwaterNaCl')
XCl, XNa, Xwater, y = variables('XCl, XNa, Xwater, y')
model = Model({y: lnfca()})
fit = Fit(model, XCl=XCl_vals, XNa=XNa_vals, Xwater=Xwater_vals, y=ln_gama_x_vals)
fit_result = fit.execute()
#*****************************************************************`

Traceback (most recent call last): File "F:\Zohreh\MainZohreh\postdoc-field\CSU\pythonProject\psc_ac_pure_symfit.py", line 404, in fit = Fit(model, XCl=XCl_vals, XNa=XNa_vals, Xwater=Xwater_vals, y=ln_gama_x_vals) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 427, in wrapped_func return func(*bound_args.args, bound_args.kwargs) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\fit.py", line 394, in init super(Fit, self).init(self.model, absolute_sigma=absolute_sigma, File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\support.py", line 427, in wrapped_func return func(*bound_args.args, *bound_args.kwargs) File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\fit.py", line 75, in init raise err File "C:\Users\Zohreh\AppData\Roaming\Python\Python310\site-packages\symfit\core\fit.py", line 63, in init bound_arguments = signature.bind(ordered_data, named_data) File "C:\Program Files\Python310\lib\inspect.py", line 3179, in bind return self._bind(args, kwargs) File "C:\Program Files\Python310\lib\inspect.py", line 3094, in _bind raise TypeError(msg) from None TypeError: missing a required argument: 'BNaCl'

[pckroon]

Ok, great! As far as I can tell symfit is behaving as intended. I won't go through your code though, that's on you. If you have a minimal non-working example I could take a look, but I'll close this issue for now.