Fitting multidimensional data with Python symfit ODEModel

I am trying to fit the parameters of an ODE to data with two dimensions, which should generally be possible, according to the example Fitting multidimensional datasets.

This is my failed attempt so far

import symfit as sf 
import numpy as np

# data 
x = np.arange(0,19) 
data = 10e-4 * np.array([8,10,12,11,10,15,25,37,46,40,43,35,27,14,8,10,13,9,10]) 
data2 = 10e-3 * np.array([0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0]) 

# model 
t, mp, inp = sf.variables('t, mp, inp') 
tau, w = sf.parameters('tau, w') 
model = {sf.D(mp, t): (-mp + inp * w) / tau } 

# fitting 
ode_model = sf.ODEModel(model, initial={inp: data2[0], t: 0.0, mp: data[0]})
fit = sf.Fit(ode_model, inp=data2, t=x, mp=data) 
fit_result = fit.execute() 

It seems I'm not correctly defining inp in the model definition. I'm getting the following error TypeError: got an unexpected keyword argument 'inp' . I suspect that I'm making a mistake in providing inp as named_data to sf.Fit() as ist does not appear to be an independent variable in the model documentation sf.Fit()

This is the full error message:

/usr/local/lib/python3.10/dist-packages/symfit/core/fit.py in __init__(self, model, objective, minimizer, constraints, absolute_sigma, *ordered_data, **named_data)
    372         # Bind as much as possible the provided arguments.
    373         signature = self._make_signature()
--> 374         bound_arguments = signature.bind_partial(*ordered_data, **named_data)
    375 
    376         # Select objective function to use. Has to be done before calling

/usr/lib/python3.10/inspect.py in bind_partial(self, *args, **kwargs)
   3191         Raises `TypeError` if the passed arguments can not be bound.
   3192         """
-> 3193         return self._bind(args, kwargs, partial=True)
   3194 
   3195     def __reduce__(self):

/usr/lib/python3.10/inspect.py in _bind(self, args, kwargs, partial)
   3173                 arguments[kwargs_param.name] = kwargs
   3174             else:
-> 3175                 raise TypeError(
   3176                     'got an unexpected keyword argument {arg!r}'.format(
   3177                         arg=next(iter(kwargs))))

TypeError: got an unexpected keyword argument 'inp'

Could someone help? Thank you so so much :-)

inp has to be defined as an expression and integrated into the expression d mp / dt. To do so, the data for inp has to be fit so as to reproduce data2. Since data2 looks like a square wave, a Fourier series is used to fit the data. The following is the code implementing these ideas:

import symfit as sf 
import numpy as np
from functools import reduce
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

from scipy.integrate import solve_ivp

def square_wave_symbolic(x, L, shift, n=8):
    return reduce(lambda a, b: a + b, [sf.sin(2*(1+2 *k) *sf.pi * (x+shift) / L)/(1+2* k) for k in range(n)])

def square_wave_numeric(x, L, shift, n=8):
    return reduce(lambda a, b: a + b, [np.sin(2*(1+2 *k) *np.pi * (x+shift) / L)/(1+2* k) for k in range(n)])

# data 
x = np.arange(0,19) 
data = 10e-4 * np.array([8,10,12,11,10,15,25,37,46,40,43,35,27,14,8,10,13,9,10]) 
data2 = 10e-3 * np.array([0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0]) 

# Fit square wave
alpha0 = 5E-3
beta0 = 5E-3
L0 = 18
shift0 = -4
n = 9

plt.figure()
plt.title('Approximation of inp using Fourier Series')
result = curve_fit(lambda x, alpha, beta, shift, L: alpha + beta * square_wave_numeric(x, L, shift, n=n), x, data2, p0=(alpha0, beta0, shift0, L0), full_output=True)
alpha, beta, shift, Lf = result[0]
plt.scatter(x, data2, label='Original data')
plt.plot(x, alpha + beta * square_wave_numeric(x, Lf, shift, n=n), label='Fit data')
plt.legend()


# model 
t, mp = sf.variables('t, mp') 
tau, w = sf.parameters('tau, w') 
inp = alpha + beta * square_wave_symbolic(t, Lf, shift, n=n)
model = {sf.D(mp, t): (-mp + inp * w) / tau } 

# fitting 
ode_model = sf.ODEModel(model, {t: 0.0, mp: data[0]})
fit = sf.Fit(ode_model, t=x, mp=data) 
fit_result = fit.execute() 
print(fit_result.params)
w = fit_result.params['w']
tau = fit_result.params['tau']

# Verify parameters by solving numerically
def f(t, x, w, tau):
    mp = x[0]
    inp = alpha + beta * square_wave_numeric(t, Lf, shift, n=n)
    return np.array([(-mp + inp * w) / tau])
    
t0 = x[0]
tf = x[-1]
t_eval = np.linspace(t0, tf, 200)
ode_result = solve_ivp(f, (t0, tf), (data[0],), t_eval=t_eval, args=(w, tau), method='Radau')

plt.figure()
plt.scatter(x, data, label='Original data')
plt.plot(ode_result.t, ode_result.y[0], label='ODE solver data')
plt.legend()

I end up with the following value for fit_result.params:

OrderedDict([('tau', 1.191491195628205), ('w', 3.4675371620653133)])

The following are the plot for inp variable fit:

ODE fit plot:

Notes:

1. The square wave function describing inp is not smooth. I am not sure about where the data comes from, but if this is just some test data, the real data may have to be fit using a different function.

2. It would make sense to do the ODE and inp function fit together. Doing it using a minimize function is feasible. I am not sure how this could be done with symfit as I am not familiar with this library.