with this code I'm trying to make a curvefitting with the scipy.optimize.curve_fit using a function that contains integrals, fitting then to some data I was given. However, when running the code it is throwing an error that seems to be related with the parameters that I'm trying to fit (the code is below). It gives me the error
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-118-13cc86454e5e> in <module>
88
89 vlc_fit = np.vectorize(lc_fit)
---> 90 popt, pcov = scipy.optimize.curve_fit(lc_fit,
91 t_data, lc_init_data)
92
~/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
782 # Remove full_output from kwargs, otherwise we're passing it in twice.
783 return_full = kwargs.pop('full_output', False)
--> 784 res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
785 popt, pcov, infodict, errmsg, ier = res
786 ysize = len(infodict['fvec'])
~/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
408 if not isinstance(args, tuple):
409 args = (args,)
--> 410 shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
411 m = shape[0]
412
~/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
22 def _check_func(checker, argname, thefunc, x0, args, numinputs,
23 output_shape=None):
---> 24 res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
25 if (output_shape is not None) and (shape(res) != output_shape):
26 if (output_shape[0] != 1):
~/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py in func_wrapped(params)
482 if transform is None:
483 def func_wrapped(params):
--> 484 return func(xdata, *params) - ydata
485 elif transform.ndim == 1:
486 def func_wrapped(params):
<ipython-input-118-13cc86454e5e> in lc_fit(t, args)
77 def lc_fit(t, args):
78
---> 79 p, a, b, T, delta, gamma1, gamma2 = args
80 # Main Loop to calculate the Flux
81
TypeError: cannot unpack non-iterable numpy.float64 object
and for what I see, I have no flaws on the code, or at least comparing it to another code given in an other stackoverflow question that fits something that uses integrals in it (here is the link to that code Fitting data with integral function)
Here is what I think is the part of my code that is having trouble:
import scipy
from scipy.misc import derivative
import scipy.integrate as integrate
import numpy as np
import matplotlib.pyplot as plt
import timeit
t_data = np.delete(tce1_time[0:200], np.argwhere(np.isnan(tce1_lc_init[0:200]))[0:], axis=None) #take away the time values where we have a nan value in lc_init
lc_init_data = np.delete(tce1_lc_init[0:200], np.argwhere(np.isnan(tce1_lc_init[0:200]))[0:], axis=None) #take away the nan values in the data array of lc_init
def lc_fit(t, args):
p, a, b, T, delta, gamma1, gamma2 = args
# Main Loop to calculate the Flux
Integral_1 = integrate.quad(integrand_1, 0.0001, 1., args=(t, p, a, b, T, delta, gamma1, gamma2))[0]
Integral_2 = integrate.quad(integrand_2, 0.0001, 1., args=(gamma1, gamma2))[0]
F= 1. - Integral_1/Integral_2
return F
popt, pcov = scipy.optimize.curve_fit(lc_fit,
t_data, lc_init_data)
and here is the entire code in case there is something else you want to check
import scipy
from scipy.misc import derivative
import scipy.integrate as integrate
import numpy as np
import matplotlib.pyplot as plt
import timeit
t_data = np.delete(tce1_time[0:200], np.argwhere(np.isnan(tce1_lc_init[0:200]))[0:], axis=None) #take away the time values where we have a nan value in lc_init
lc_init_data = np.delete(tce1_lc_init[0:200], np.argwhere(np.isnan(tce1_lc_init[0:200]))[0:], axis=None) #take away the nan values in the data array of lc_init
def L(r, P, Z): #it is already multiplied by r² as it is in the derivate into the integral
'''
Obstruction function
'''
p = P/r
z0 = Z/r
z = np.zeros(len(z0))
for i in range(len(z0)):
# Reflects the information w.r.t. the y axis
if z0[i]>0:
z[i] = z0[i]
if z0[i]<0:
z[i] = -z0[i]
if z[i]>1+p:
return 0.
elif (abs(1-p)<z[i] and z[i]<=1+p):
k0 = np.arccos((p**2 + z[i]**2 -1)/(2*p*z[i]))
k1 = np.arccos((1-p**2 + z[i]**2)/(2*z[i]))
L0 = k0*p**2
L2 = np.sqrt((4*z[i]**2- (1+z[i]**2-p**2)**2)/4)
return (L0 + k1 - L2)*r**2/np.pi
elif (z[i]<=1-p):
return p**2*r**2
elif z[i]<= p-1:
return 1.*r**2
def I_function(r, gamma1, gamma2):
'''
Quadratic limb-darkening function
'''
mu = np.sqrt(1-r**2)
return 1. - gamma1*(1-mu) - gamma2*(1-mu)**2
def integrand_1(r, args1):
'''
Integrand in the numerator
'''
t, p, a, b, T, delta, gamma1, gamma2 = args1
# p = 0.1 #radius ratio
# b = 0.7 #impact. parameter
# gamma1 = 0.296 #linear limb darkening
# gamma2 = 0.34 #quadratic limb darkening
# a = 2. #normalized semi-major axis (normalized with the star radius)
# T = 10. #orbital period
# delta = 0. #orbital phase
omega = 2*np.pi/T #angular velocity
x = a*np.cos(omega*t+delta)
z = np.sqrt(x**2 + b**2)
T1 = derivative(L, r, dx=1e-6, args=(p,z))
T2 = I_function(r, gamma1, gamma2)
return T1*T2
def integrand_2(r, args2):
'''
Integrand in the denominator
'''
gamma1, gamma2 = args2
return I_function(r, gamma1, gamma2)*2*r
def lc_fit(t, args):
p, a, b, T, delta, gamma1, gamma2 = args
# Main Loop to calculate the Flux
Integral_1 = integrate.quad(integrand_1, 0.0001, 1., args=(t, p, a, b, T, delta, gamma1, gamma2))[0]
Integral_2 = integrate.quad(integrand_2, 0.0001, 1., args=(gamma1, gamma2))[0]
F= 1. - Integral_1/Integral_2
return F
popt, pcov = scipy.optimize.curve_fit(lc_fit,
t_data, lc_init_data)
