Long time lurking, first time posting.
I am working with a chemical system that is detected only for a certain period of time, so I will have the reaction and the decay of the signal. The equation is given by:
Derivative(GL, t): (-k*GL) - GL/a,
Derivative(GM, t): (k*GL) - GM/b,
I have managed to fit my data by using symfit package (image below to give an idea of the system), however since I will need to do Monte Carlo simulation, I need to fit my data using scipy. 
I have tried to define the equation in this way:
def f(C, xdata):
GL = ydataScaled
GM = ydataScaled2
dGLdt = -k*GL - GL/a
dGMdt = k*GL - GM/b
return [dGLdt, dGMdt]
However, I am not able to fit neither by using optimize.minimize or odeint. What would be the right approach in this case to fit two dataset in y that share some parameters?
Full code:
import nmrglue as ng
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from scipy import integrate
from scipy.optimize import curve_fit
from scipy.integrate import odeint
from symfit import variables, parameters, Fit, ODEModel, Derivative, D, exp, sin, Model, cos, integrate
# read in the bruker formatted data
dic,data = ng.bruker.read_pdata('/opt/topspin4.1.0/NMR/2021_09_27_Glutamine/90/pdata/1')
#Bruker to NMRPipe data
C = ng.convert.converter()
C.from_bruker(dic, data)
pdic, ppdata = C.to_pipe()
#process the spectrum
ZF_Number = 16384
ppdata = ng.proc_base.di(ppdata) # discard the imaginaries
show = ppdata[2] #show the spectra number X
# determind the ppm scale
udic = ng.bruker.guess_udic(dic, data)
uc = ng.fileiobase.uc_from_udic(udic)
ppm_scale = uc.ppm_scale()
ppms = uc.ppm_scale()
#Plot the spectra
fig1 = plt.figure()
bx = fig1.add_subplot(111)
bx.plot(ppms, show)
plt.xlabel('Chemical Shift (ppm)')
plt.ylabel('Intensity')
First = 0
End = 80
#Integration for every i in the range
Area = []
Area2 = []
Area3 = [] #noise measurement, using the same chemical shift lenght as the product-peak.
#limits = [(176, 180), (180, 183)]
for i in range(First,End):
Area.append(ng.analysis.integration.integrate(ppdata[i], uc, (177.15, 177.80), unit = "ppm", noise_limits = None, norm_to_range = None, calibrate = 1.0))
NP_Area = np.asarray(Area)
for i in range(First, End):
Area2.append(ng.analysis.integration.integrate(ppdata[i], uc, (180.80, 181.10), unit = "ppm", noise_limits = None, norm_to_range = None, calibrate = 1.0))
NP_Area2 = np.asarray(Area2)
for i in range(First, End):
Area3.append(ng.analysis.integration.integrate(ppdata[i], uc, (20.0, 20.3), unit = "ppm", noise_limits = None, norm_to_range = None, calibrate = 1.0))
NP_Area3 = np.asarray(Area3)
#Plot the buildUP
fig2 = plt.figure()
cx = fig2.add_subplot(111)
cx.plot(NP_Area)
cx.plot(NP_Area2)
plt.xlabel('Time (seconds)')
plt.ylabel('Intensity')
#Fitting
d1 = dic['acqus']['D'][1]
xdata = (np.arange(First, End) - First)*d1
ydata = NP_Area[:,0]
ydata2 = NP_Area2[:,0]
ydataScaled = ydata/max(ydata) #normalized to the initial value of the Glu signal to compensate for any variations in the polarization level
ydataScaled2 = ydata2/max(ydata) # same as above
#GL, GM, t = variables('GL, GM, t')
a, b, k = parameters('a, b, k')
# Define the equation considering the enzymatic reaction Gl -> Gm with the HP decay.
def f(C, xdata):
GL = ydataScaled
GM = ydataScaled2
dGLdt = -k*GL - GL/a
dGMdt = k*GL - GM/b
return [dGLdt, dGMdt]
C0 = [1, 0]
popt, pcov = sp.optimize.minimize(f, xdata, args = (ydataScaled, ydataScaled2))
And the error:
runfile('/Users/karensantos/Desktop/Codes/Stack_question.py', wdir='/Users/karensantos/Desktop/Codes')
2
(512, 32768)
float64
/opt/anaconda3/lib/python3.8/site-packages/nmrglue/fileio/convert.py:68: UserWarning: Incompatible dtypes, conversion not recommended
warn("Incompatible dtypes, conversion not recommended")
Traceback (most recent call last):
File "/Users/karensantos/Desktop/Codes/Stack_question.py", line 112, in <module>
popt, pcov = sp.optimize.minimize(f, xdata, args = (ydataScaled, ydataScaled2))
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_minimize.py", line 612, in minimize
return _minimize_bfgs(fun, x0, args, jac, callback, **options)
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/optimize.py", line 1101, in _minimize_bfgs
sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/optimize.py", line 261, in _prepare_scalar_function
sf = ScalarFunction(fun, x0, args, grad, hess,
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 76, in __init__
self._update_fun()
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 166, in _update_fun
self._update_fun_impl()
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 73, in update_fun
self.f = fun_wrapped(self.x)
File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 70, in fun_wrapped
return fun(x, *args)
TypeError: f() takes 2 positional arguments but 3 were given
