Python scipy.optimize.leastsq() Examples

def gaussian_fit_curve(x, y, mu0=0, sigma0=1, a0=None, return_all=False,
        **kwargs):
    """Gaussian fit of curve (x,y).
    If a0 is None then only (mu,sigma) are fitted (to a gaussian density).
    `kwargs` are passed to the leastsq() function.

    If return_all=False then return only the fitted (mu,sigma) values
    If return_all=True (or full_output=True is passed to leastsq)
    then the full output of leastsq is returned.
    """
    if a0 is None:
        gauss_pdf = lambda x, m, s: np.exp(-((x-m)**2)/(2*s**2))/\
                                    (np.sqrt(2*np.pi)*s)
        err_fun = lambda p, x, y: gauss_pdf(x, *p) - y
        res = leastsq(err_fun, x0=[mu0, sigma0], args=(x, y), **kwargs)
    else:
        gauss_fun = lambda x, m, s, a: a*np.sign(s)*np.exp(-((x-m)**2)/(2*s**2))
        err_fun = lambda p, x, y: gauss_fun(x, *p) - y
        res = leastsq(err_fun, x0=[mu0, sigma0, a0], args=(x, y), **kwargs)

    if 'full_output' in kwargs:
        return_all = kwargs['full_output']
    mu, sigma = res[0][0], res[0][1]
    if return_all: return res
    return mu, sigma 

def _fitGivenSmoothness(y, t, ode, theta, s):
    # p = ode.getNumParam()
    # d = ode.getNumState()

    splineList = interpolate(y, t, s=s)
    interval = np.linspace(t[1], t[-1], 1000)
    # xApprox, fxApprox, t = _getApprox(splineList, interval)

    # g2 = partial(residual_sample, ode, fxApprox, xApprox, interval)
    # g2J = partial(jac_sample, ode, fxApprox, xApprox, interval)
    g2 = partial(residual_interpolant, ode, splineList, interval)
    g2J = partial(jac_interpolant, ode, splineList, interval)

    res = leastsq(func=g2, x0=theta, Dfun=g2J, full_output=True)

    loss = 0
    for spline in splineList:
        loss += spline.get_residual()
    # approximate the integral using fixed points
    r = np.reshape(res[2]['fvec']**2, (len(interval), len(splineList)), 'F')

    return (r.sum())*(interval[1] - interval[0]) + loss 

def error(params, xs, ys):
        """A callback function for use with :func:`scipy.optimize.leastsq` to compute
        the residuals given a set of parameters, x, and y

        Parameters
        ----------
        params : list
            The model's parameters, a list of floats
        xs : :class:`np.ndarray`
            The time array to predict with respect to.
        ys : :class:`np.ndarray`
            The intensity array to predict against

        Returns
        -------
        :class:`np.ndarray`:
            The residuals of ``ys - self.shape(params, x)`` with optional penalty terms
        """
        raise NotImplementedError() 

def guess(xs, ys):
        """Get crude estimates of roughly where to start fitting parameters

        The results are by no means accurate, but will serve as a reasonable starting point
        for :func:`scipy.optimize.leastsq`.

        Parameters
        ----------
        xs : :class:`np.ndarray`
            The time array to predict with respect to.
        ys : :class:`np.ndarray`
            The intensity array to predict against

        Returns
        -------
        list
        """
        raise NotImplementedError() 

def test_regression_2639(self):
        # This test fails if epsfcn in leastsq is too large.
        x = [574.14200000000005, 574.154, 574.16499999999996,
             574.17700000000002, 574.18799999999999, 574.19899999999996,
             574.21100000000001, 574.22199999999998, 574.23400000000004,
             574.245]
        y = [859.0, 997.0, 1699.0, 2604.0, 2013.0, 1964.0, 2435.0,
             1550.0, 949.0, 841.0]
        guess = [574.1861428571428, 574.2155714285715, 1302.0, 1302.0,
                 0.0035019999999983615, 859.0]
        good = [ 5.74177150e+02, 5.74209188e+02, 1.74187044e+03, 1.58646166e+03,
                 1.0068462e-02, 8.57450661e+02]

        def f_double_gauss(x, x0, x1, A0, A1, sigma, c):
            return (A0*np.exp(-(x-x0)**2/(2.*sigma**2))
                    + A1*np.exp(-(x-x1)**2/(2.*sigma**2)) + c)
        popt, pcov = curve_fit(f_double_gauss, x, y, guess, maxfev=10000)
        assert_allclose(popt, good, rtol=1e-5) 

def _do_fit_fano(self, amplitudes_sq):
        # initial guess
        bw = 1e6
        q  = 1 #np.sqrt(1-amplitudes_sq).min()  # 1-Amp_sq = 1-1+q^2  => A_min = q
        fr = self._fit_frequency[np.argmin(amplitudes_sq)]
        a  = amplitudes_sq.max()

        p0 = [q, bw, fr, a]

        def fano_residuals(p,frequency,amplitude_sq):
            q, bw, fr, a = p
            err = amplitude_sq-self._fano_reflection(frequency,q,bw,fr=fr,a=a)
            return err

        p_fit = leastsq(fano_residuals,p0,args=(self._fit_frequency,np.array(amplitudes_sq)))
        #print(("q:%g bw:%g fr:%g a:%g")% (p_fit[0][0],p_fit[0][1],p_fit[0][2],p_fit[0][3]))
        return p_fit[0] 

def _fit_circle_iter(self,z_data, xc, yc, rc):
        '''
        this is the radial weighting procedure
        it improves your fitting value for the radius = Ql/Qc
        use this to improve your fit in presence of heavy noise
        after having used the standard algebraic fir_circle() function
        the weight here is: W=1/sqrt((xc-xi)^2+(yc-yi)^2)
        this works, because the center of the circle is usually much less
        corrupted by noise than the radius
        '''
        xdat = z_data.real
        ydat = z_data.imag
        def fitfunc(x,y,xc,yc):
            return np.sqrt((x-xc)**2+(y-yc)**2)
        def residuals(p,x,y):
            xc,yc,r = p
            temp = (r-fitfunc(x,y,xc,yc))
            return temp
        p0 = [xc,yc,rc]
        p_final = spopt.leastsq(residuals,p0,args=(xdat,ydat))
        xc,yc,rc = p_final[0]
        return xc,yc,rc 

def fit_spot(spot):
    size = spot.shape[0]
    size_half = int(size / 2)
    grid = _np.arange(-size_half, size_half + 1, dtype=_np.float32)
    model_x = _np.empty(size, dtype=_np.float32)
    model_y = _np.empty(size, dtype=_np.float32)
    model = _np.empty((size, size), dtype=_np.float32)
    residuals = _np.empty((size, size), dtype=_np.float32)
    # theta is [x, y, photons, bg, sx, sy]
    theta0 = _initial_parameters(spot, size, size_half)
    args = (spot, grid, size, model_x, model_y, model, residuals)
    result = _optimize.leastsq(
        _compute_residuals, theta0, args=args, ftol=1e-2, xtol=1e-2
    )  # leastsq is much faster than least_squares
    """
    model = compute_model(result[0], grid, size, model_x, model_y, model)
    plt.figure()
    plt.subplot(121)
    plt.imshow(spot, interpolation='none')
    plt.subplot(122)
    plt.imshow(model, interpolation='none')
    plt.colorbar()
    plt.show()
    """
    return result[0] 

def fit_random(self, ntries=10, rvs_generator=None, nparams=None):
        '''fit with random starting values

        this could be replaced with a global fitter

        '''

        if nparams is None:
                nparams = self.nparams
        if rvs_generator is None:
            rvs = np.random.uniform(low=-10, high=10, size=(ntries, nparams))
        else:
            rvs = rvs_generator(size=(ntries, nparams))

        results = np.array([np.r_[self.fit_minimal(rv),  rv] for rv in rvs])
        #selct best results and check how many solutions are within 1e-6 of best
        #not sure what leastsq returns
        return results 

def fit_random(self, ntries=10, rvs_generator=None, nparams=None):
        '''fit with random starting values

        this could be replaced with a global fitter

        '''

        if nparams is None:
                nparams = self.nparams
        if rvs_generator is None:
            rvs = np.random.uniform(low=-10, high=10, size=(ntries, nparams))
        else:
            rvs = rvs_generator(size=(ntries, nparams))

        results = np.array([np.r_[self.fit_minimal(rv),  rv] for rv in rvs])
        #selct best results and check how many solutions are within 1e-6 of best
        #not sure what leastsq returns
        return results 

def _residual(param, radial, profile, previous):
    """ `scipy.optimize.leastsq` residuals function.

        Evaluate the difference between a radial-scaled intensity profile
        and its adjacent "previous" angular slice.

    """

    radial_scaling, amplitude = param[0], param[1]

    newradial = radial * radial_scaling
    spline_prof = UnivariateSpline(newradial, profile, s=0, ext=3)
    newprof = spline_prof(radial) * amplitude

    # residual cf adjacent slice profile
    return newprof - previous 

def _fit_entire_model(self,f_data,z_data,fr,absQc,Ql,phi0,delay,a=1.,alpha=0.,maxiter=0):
        '''
        fits the whole model: a*exp(i*alpha)*exp(-2*pi*i*f*delay) * [ 1 - {Ql/Qc*exp(i*phi0)} / {1+2*i*Ql*(f-fr)/fr} ]
        '''
        def funcsqr(p,x):
            fr,absQc,Ql,phi0,delay,a,alpha = p
            return np.array([np.absolute( ( a*np.exp(np.complex(0,alpha))*np.exp(np.complex(0,-2.*np.pi*delay*x[i])) * ( 1 - (Ql/absQc*np.exp(np.complex(0,phi0)))/(np.complex(1,2*Ql*(x[i]-fr)/fr)) )  ) )**2 for i in range(len(x))])
        def residuals(p,x,y):
            fr,absQc,Ql,phi0,delay,a,alpha = p
            err = [np.absolute( y[i] - ( a*np.exp(np.complex(0,alpha))*np.exp(np.complex(0,-2.*np.pi*delay*x[i])) * ( 1 - (Ql/absQc*np.exp(np.complex(0,phi0)))/(np.complex(1,2*Ql*(x[i]-fr)/fr)) )  ) ) for i in range(len(x))]
            return err
        p0 = [fr,absQc,Ql,phi0,delay,a,alpha]
        (popt, params_cov, infodict, errmsg, ier) = spopt.leastsq(residuals,p0,args=(np.array(f_data),np.array(z_data)),full_output=True,maxfev=maxiter)
        len_ydata = len(np.array(f_data))
        if (len_ydata > len(p0)) and params_cov is not None:  #p_final[1] is cov_x data  #this caculation is from scipy curve_fit routine - no idea if this works correctly...
            s_sq = (funcsqr(popt, np.array(f_data))).sum()/(len_ydata-len(p0))
            params_cov = params_cov * s_sq
        else:
            params_cov = np.inf
        return popt, params_cov, infodict, errmsg, ier
    
    # 

def fit(self):
        """Fit E(V) model, fill ``self.params``."""
        # Quadratic fit to get an initial guess for the parameters.
        # Thanks: https://github.com/materialsproject/pymatgen
        # -> pymatgen/io/abinitio/eos.py
        a, b, c = np.polyfit(self.volume, self.energy, 2)
        v0 = -b/(2*a)
        e0 = a*v0**2 + b*v0 + c
        b0 = 2*a*v0
        b1 = 4  # b1 is usually a small number like 4
        if not self.volume.min() < v0 and v0 < self.volume.max():
            raise Exception('The minimum volume of a fitted parabola is not in the input volumes')

        # need to use lst2dct and dct2lst here to keep the order of parameters
        pp0_dct = dict(e0=e0, b0=b0, b1=b1, v0=v0)
        target = lambda pp, v: self.energy - self.func(v, self.func.lst2dct(pp))
        pp_opt, ierr = leastsq(target,
                               self.func.dct2lst(pp0_dct),
                               args=(self.volume,))
        self.params = self.func.lst2dct(pp_opt) 

def gaussian_fit_pdf(s, mu0=0, sigma0=1, a0=1, return_all=False,
        leastsq_kwargs={}, **kwargs):
    """Gaussian fit of samples s using a fit to the empirical PDF.
    If a0 is None then only (mu,sigma) are fitted (to a gaussian density).
    `kwargs` are passed to get_epdf().
    If return_all=False then return only the fitted (mu,sigma) values
    If return_all=True (or full_output=True is passed to leastsq)
    then the full output of leastsq and the PDF curve is returned.
    """
    ## Empirical PDF
    epdf = get_epdf(s, **kwargs)

    res = gaussian_fit_curve(epdf[0], epdf[1], mu0, sigma0, a0, return_all,
            **leastsq_kwargs)
    if return_all: return res, epdf
    return res 

def gaussian_fit_cdf(s, mu0=0, sigma0=1, return_all=False, **leastsq_kwargs):
    """Gaussian fit of samples s fitting the empirical CDF.
    Additional kwargs are passed to the leastsq() function.
    If return_all=False then return only the fitted (mu,sigma) values
    If return_all=True (or full_output=True is passed to leastsq)
    then the full output of leastsq and the histogram is returned.
    """
    ## Empirical CDF
    ecdf = [np.sort(s), np.arange(0.5, s.size+0.5)*1./s.size]

    ## Analytical Gaussian CDF
    gauss_cdf = lambda x, mu, sigma: 0.5*(1+erf((x-mu)/(np.sqrt(2)*sigma)))

    ## Fitting the empirical CDF
    err_func = lambda p, x, y: y - gauss_cdf(x, p[0], p[1])
    res = leastsq(err_func, x0=[mu0, sigma0], args=(ecdf[0], ecdf[1]),
            **leastsq_kwargs)
    if return_all: return res, ecdf
    return res[0] 

def two_gaussian_fit_curve(x, y, p0, return_all=False, verbose=False, **kwargs):
    """Fit a 2-gaussian mixture to the (x,y) curve.
    `kwargs` are passed to the leastsq() function.

    If return_all=False then return only the fitted paramaters
    If return_all=True then the full output of leastsq is returned.
    """
    if kwargs['method'] == 'leastsq':
        kwargs.pop('method')
        kwargs.pop('bounds')
        def err_func(p, x, y):
            return (y - two_gauss_mix_pdf(x, p))
        res = leastsq(err_func, x0=p0, args=(x, y), **kwargs)
        p = res[0]
    else:
        def err_func(p, x, y):
            return ((y - two_gauss_mix_pdf(x, p))**2).sum()
        res = minimize(err_func, x0=p0, args=(x, y), **kwargs)
        p = res.x

    if verbose:
        print(res, '\n')
    if return_all: return res
    return reorder_parameters(p) 

def gaussian2d_fit(sx, sy, guess=[0.5,1]):
    """2D-Gaussian fit of samples S using a fit to the empirical CDF."""
    assert sx.size == sy.size

    ## Empirical CDF
    ecdfx = [np.sort(sx), np.arange(0.5,sx.size+0.5)*1./sx.size]
    ecdfy = [np.sort(sy), np.arange(0.5,sy.size+0.5)*1./sy.size]

    ## Analytical gaussian CDF
    gauss_cdf = lambda x, mu, sigma: 0.5*(1+erf((x-mu)/(np.sqrt(2)*sigma)))

    ## Fitting the empirical CDF
    fitfunc = lambda p, x: gauss_cdf(x, p[0], p[1])
    errfunc = lambda p, x, y: fitfunc(p, x) - y
    px,v = leastsq(errfunc, x0=guess, args=(ecdfx[0],ecdfx[1]))
    py,v = leastsq(errfunc, x0=guess, args=(ecdfy[0],ecdfy[1]))
    print("2D Gaussian CDF fit", px, py)

    mux, sigmax = px[0], px[1]
    muy, sigmay = py[0], py[1]
    return mux, sigmax, muy, sigmay 

def fit(self, task=0, s=None, t=None, full_output=1, warn=False):
        """Fit the function to the data"""
        if self.function == 'spline':
            self.set_weight(self.yerr)
            self.results = interpolate.splrep(self.x, self.y, w=self.weight,
                                              task=0, s=None, t=None,
                                              k=self.order,
                                              full_output=full_output)
            # w=None, k=self.order, s=s, t=t, task=task,
            # full_output=full_output)
            self.set_coef(self.results[0])

        else:
            self.results = optimize.leastsq(self.erf, self.coef,
                                            args=(self.x, self.y, self.yerr),
                                            full_output=full_output)
            self.set_coef(self.results[0]) 

def LinearB(Xi, Yi):
    X = np.asfarray(Xi)
    Y = np.asfarray(Yi)

    # we want a function y = m * x + b
    def fp(v, x):
        return x * v[0] + v[1]

    # the error of the function e = x - y
    def e(v, x, y):
        return (fp(v, x) - y)

    # the initial value of m, we choose 1, because we thought YODA would
    # have chosen 1
    v0 = np.array([1.0, 1.0])

    vr, _success = leastsq(e, v0, args=(X, Y))

    # compute the R**2 (sqrt of the mean of the squares of the errors)
    err = np.sqrt(sum(np.square(e(vr, X, Y))) / (len(X) * len(X)))

#    print vr, success, err
    return vr, err 

def Regresion(self):
        t=array(self.KEq_Tab.getColumn(0)[:-1])
        k=array(self.KEq_Tab.getColumn(1)[:-1])
        if len(t)>=4:
            if 4<=len(t)<8:
                inicio=r_[0, 0, 0, 0]
                f=lambda par, T: exp(par[0]+par[1]/T+par[2]*log(T)+par[3]*T)
                resto=lambda par, T, k: k-f(par, T)
            else:
                inicio=r_[0, 0, 0, 0, 0, 0, 0, 0]
                f=lambda par, T: exp(par[0]+par[1]/T+par[2]*log(T)+par[3]*T+par[4]*T**2+par[5]*T**3+par[6]*T**4+par[7]*T**5)
                resto=lambda par, T, k: k-f(par, T)

            ajuste=leastsq(resto,inicio,args=(t, k))
            kcalc=f(ajuste[0], t)
            error=(k-kcalc)/k*100
            self.KEq_Dat.setColumn(0, ajuste[0])
            self.KEq_Tab.setColumn(2, kcalc)
            self.KEq_Tab.setColumn(3, error)

            if ajuste[1] in [1, 2, 3, 4]:
                self.ajuste=ajuste[0] 

def _fit_skewed_lorentzian(self,f_data,z_data):
        amplitude = np.absolute(z_data)
        amplitude_sqr = amplitude**2
        A1a = np.minimum(amplitude_sqr[0],amplitude_sqr[-1])
        A3a = -np.max(amplitude_sqr)
        fra = f_data[np.argmin(amplitude_sqr)]
        def residuals(p,x,y):
            A2, A4, Ql = p
            err = y -(A1a+A2*(x-fra)+(A3a+A4*(x-fra))/(1.+4.*Ql**2*((x-fra)/fra)**2))
            return err
        p0 = [0., 0., 1e3]
        p_final = spopt.leastsq(residuals,p0,args=(np.array(f_data),np.array(amplitude_sqr)))
        A2a, A4a, Qla = p_final[0]
    
        def residuals2(p,x,y):
            A1, A2, A3, A4, fr, Ql = p
            err = y -(A1+A2*(x-fr)+(A3+A4*(x-fr))/(1.+4.*Ql**2*((x-fr)/fr)**2))
            return err
        p0 = [A1a, A2a , A3a, A4a, fra, Qla]
        p_final = spopt.leastsq(residuals2,p0,args=(np.array(f_data),np.array(amplitude_sqr)))
        #A1, A2, A3, A4, fr, Ql = p_final[0]
        #print(p_final[0][5])
        return p_final[0] 

def test_regression_2639(self):
        # This test fails if epsfcn in leastsq is too large.
        x = [574.14200000000005, 574.154, 574.16499999999996,
             574.17700000000002, 574.18799999999999, 574.19899999999996,
             574.21100000000001, 574.22199999999998, 574.23400000000004,
             574.245]
        y = [859.0, 997.0, 1699.0, 2604.0, 2013.0, 1964.0, 2435.0,
             1550.0, 949.0, 841.0]
        guess = [574.1861428571428, 574.2155714285715, 1302.0, 1302.0,
                 0.0035019999999983615, 859.0]
        good = [5.74177150e+02, 5.74209188e+02, 1.74187044e+03, 1.58646166e+03,
                1.0068462e-02, 8.57450661e+02]

        def f_double_gauss(x, x0, x1, A0, A1, sigma, c):
            return (A0*np.exp(-(x-x0)**2/(2.*sigma**2))
                    + A1*np.exp(-(x-x1)**2/(2.*sigma**2)) + c)
        popt, pcov = curve_fit(f_double_gauss, x, y, guess, maxfev=10000)
        assert_allclose(popt, good, rtol=1e-5) 

def _profileObtainAndVerify(f, df, x0, full_output=False):
    '''
    Find the solution of the profile likelihood and check
    that the algorithm has converged.
    '''
    x, cov, infodict, mesg, ier = leastsq(func=f, x0=x0, Dfun=df,
                                          maxfev=10000, full_output=True)

    if ier not in (1, 2, 3, 4):
        raise EstimateError("Failure in estimation of the profile likelihood: "
                            + mesg)

    if full_output:
        output = dict()
        output['cov'] = cov
        output['infodict'] = infodict
        output['mesg'] = mesg
        output['ier'] = ier
        return x, output
    else:
        return x 

def gaussian_fit_hist(s, mu0=0, sigma0=1, a0=None, bins=np.r_[-0.5:1.5:0.001],
        return_all=False, leastsq_kwargs={}, weights=None, **kwargs):
    """Gaussian fit of samples s fitting the hist to a Gaussian function.
    If a0 is None then only (mu,sigma) are fitted (to a gaussian density).
    kwargs are passed to the histogram function.
    If return_all=False then return only the fitted (mu,sigma) values
    If return_all=True (or full_output=True is passed to leastsq)
    then the full output of leastsq and the histogram is returned.
    `weights` optional weights for the histogram.
    """
    histogram_kwargs = dict(bins=bins, density=True, weights=weights)
    histogram_kwargs.update(**kwargs)
    H = np.histogram(s, **histogram_kwargs)
    x, y = 0.5*(H[1][:-1] + H[1][1:]), H[0]
    #bar(H[1][:-1], H[0], H[1][1]-H[1][0], alpha=0.3)

    res = gaussian_fit_curve(x, y, mu0, sigma0, a0, return_all,
            **leastsq_kwargs)
    if return_all: return res, H, x, y
    return res 

def Ajustar_Curvas_Caracteristicas(self):
        """Define the characteristic curve of pump, all input arrays must be
        of same dimension
            Q: volumetric flow, m3/s
            h: head, m
            Pot: power, hp
            NPSHr: net power suption head requered to avoid pump cavitation
        """
        Q = r_[self.curvaActual[2]]
        h = r_[self.curvaActual[3]]
        Pot = r_[self.curvaActual[4]]
        NPSH = r_[self.curvaActual[5]]

        # Function to fix
        def funcion(p, x):
            return p[0]*x**2+p[1]*x+p[2]

        # Residue
        def residuo(p, x, y):
            return funcion(p, x) - y

        inicio = r_[1, 1, 1]

        ajuste_h, exito_h = optimize.leastsq(residuo, inicio, args=(Q, h))
        self.CurvaHQ = ajuste_h

        ajuste_P, exito_P = optimize.leastsq(residuo, inicio, args=(Q, Pot))
        self.CurvaPotQ = ajuste_P

        def funcion_NPSH(p, x):
            return p[0]+p[1]*exp(p[2]*x)

        def residuo_NPSH(p, x, y):
            return funcion_NPSH(p, x) - y
        ajuste_N, ex = optimize.leastsq(residuo_NPSH, inicio, args=(Q, NPSH))
        self.CurvaNPSHQ = ajuste_N 

def fit_minimal(self, start_value):
        '''minimal fitting with no extra calculations'''
        func = self.geterrors
        res = optimize.leastsq(func, start_value, full_output=0, **kw)
        return res 

def __call__(self, eq):
        """
        Apply constraints to Equilibrium eq
        """
        
        tokamak = eq.getMachine()
        start_currents = tokamak.controlCurrents()
        
        end_currents, _ = optimize.leastsq(self.psi_difference, start_currents, args=(eq,))
        
        tokamak.setControlCurrents(end_currents)

        # Ensure that the last constraint used is set in the Equilibrium
        eq._constraints = self 

def __call__(self, eq):
        """
        Apply constraints to Equilibrium eq
        """
        
        tokamak = eq.getMachine()
        start_currents = tokamak.controlCurrents()
        
        end_currents, _ = optimize.leastsq(self.psinorm_difference, start_currents, args=(eq,))
        
        tokamak.setControlCurrents(end_currents)
        
        # Ensure that the last constraint used is set in the Equilibrium
        eq._constraints = self 

def _epd_residual(coeffs, mags, fsv, fdv, fkv, xcc, ycc, bgv, bge, iha, izd):
    '''
    This is the residual function to minimize using scipy.optimize.leastsq.

    '''

    f = _epd_function(coeffs, fsv, fdv, fkv, xcc, ycc, bgv, bge, iha, izd)
    residual = mags - f
    return residual 

def curve_Predicted(x, T):
    """Fill the missing point of a distillation curve

    Parameters
    ------------
    x : list
        Array with mole/volume/weight fraction, [-]
    T : list
        Array with boiling point temperatures, [K]

    Returns
    -------
    TBP : list
        Calculated distillation data

    Examples
    --------
    Example 4.7 from [54]_

    >>> xw = [0.261, 0.254, 0.183, 0.14, 0.01, 0.046, 0.042, 0.024, 0.015, \
        0.009, 0.007]
    >>> x = [sum(xw[:i+1]) for i, xi in enumerate(xw)]
    >>> T = [365, 390, 416, 440, 461, 482, 500, 520, 539, 556, 573]
    >>> parameter, r = curve_Predicted(x, T)
    >>> "%.0f %.4f %.4f" % tuple(parameter)
    '327 0.2028 1.3802'
    """
    x = array(x)
    T = array(T)
    p0 = [T[0], 0.1, 1]

    def errf(p, xw, Tb):
        return _Tb_Predicted(p, xw) - Tb

    p, cov, info, mesg, ier = leastsq(errf, p0, args=(x, T), full_output=True)

    ss_err = (info['fvec']**2).sum()
    ss_tot = ((T-T.mean())**2).sum()
    r = 1-(ss_err/ss_tot)
    return p, r