How to use the lmfit.minimize function in lmfit

mi = lmfit.minimize(residual, p, method='nelder', nan_policy='omit')
lmfit.printfuncs.report_fit(mi.params, min_correl=0.5)
if HASPYLAB:
    plt.figure()
    plt.plot(x, y, 'b')
    plt.plot(x, residual(mi.params) + y, 'r', label='best fit')
    plt.legend(loc='best')
    plt.show()

# Place bounds on the ln(sigma) parameter that emcee will automatically add
# to estimate the true uncertainty in the data since is_weighted=False
mi.params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2))

res = lmfit.minimize(residual, method='emcee', nan_policy='omit', burn=300,
                     steps=1000, thin=20, params=mi.params, is_weighted=False,
                     progress=False)

if HASPYLAB and HASCORNER:
    emcee_corner = corner.corner(res.flatchain, labels=res.var_names,
                                 truths=list(res.params.valuesdict().values()))
    plt.show()

if HASPYLAB:
    plt.plot(res.acceptance_fraction)
    plt.xlabel('walker')
    plt.ylabel('acceptance fraction')
    plt.show()

if hasattr(res, "acor"):
    print("Autocorrelation time for the parameters:")

y = (3.0*np.exp(-x/2) - 5.0*np.exp(-(x-0.1) / 10.) +
     0.1*np.random.randn(x.size))
if HASPYLAB:
    plt.plot(x, y, 'b')
    plt.show()

p = lmfit.Parameters()
p.add_many(('a1', 4), ('a2', 4), ('t1', 3), ('t2', 3., True))


def residual(p):
    v = p.valuesdict()
    return v['a1']*np.exp(-x/v['t1']) + v['a2']*np.exp(-(x-0.1) / v['t2']) - y


mi = lmfit.minimize(residual, p, method='nelder', nan_policy='omit')
lmfit.printfuncs.report_fit(mi.params, min_correl=0.5)
if HASPYLAB:
    plt.figure()
    plt.plot(x, y, 'b')
    plt.plot(x, residual(mi.params) + y, 'r', label='best fit')
    plt.legend(loc='best')
    plt.show()

# Place bounds on the ln(sigma) parameter that emcee will automatically add
# to estimate the true uncertainty in the data since is_weighted=False
mi.params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2))

res = lmfit.minimize(residual, method='emcee', nan_policy='omit', burn=300,
                     steps=1000, thin=20, params=mi.params, is_weighted=False,
                     progress=False)

ai, aj, ak = mat2euler(solution.rotation_matrix)

    params = lmfit.Parameters()
    params.add("center_x", value=solution.center_x, vary=vary_center)
    params.add("center_y", value=solution.center_y, vary=vary_center)
    params.add("ai", value=ai, vary=vary_angles)
    params.add("aj", value=aj, vary=vary_angles)
    params.add("ak", value=ak, vary=vary_angles)
    params.add("scale", value=solution.scale, vary=vary_scale, min=0.8, max=1.2)

    wavelength = get_electron_wavelength(accelarating_voltage)
    camera_length = camera_length * 1e10
    args = k_xy, lattice_recip, wavelength, camera_length

    res = lmfit.minimize(objfunc, params, args=args, method=method)

    if verbose:  # pragma: no cover
        lmfit.report_fit(res)

    p = res.params

    ai, aj, ak = p["ai"].value, p["aj"].value, p["ak"].value
    scale = p["scale"].value
    center_x = params["center_x"].value
    center_y = params["center_y"].value

    rotation_matrix = euler2mat(ai, aj, ak)

    k_xy = (k_xy + np.array((center_x, center_y)) * scale)
    cart = detector_to_fourier(k_xy, wavelength=wavelength, camera_length=camera_length)

def time_minimize(self):
        np.random.seed(201)
        x = np.linspace(0, 15, 601)

        data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) +
                np.random.normal(size=len(x), scale=0.3) )
        params = Parameters()
        params.add('amp',   value= 1,  min=0, max=100)
        params.add('decay', value= 0.0, min=0, max=10)
        params.add('shift', value= 0.0, min=-np.pi/2., max=np.pi/2)
        params.add('omega', value= 1.0, min=0, max=10)

        return minimize(obj_func, params, args=(x, data))

('l1',   5,   True,  0,      50,  None),
                    ('e1',   0.5,   True,  0,      1,  None),
                    ('a2',   100,   True,  0,      None,  None),
                    ('f2',   1304,   True, 1300,    1330,  None),
                    ('l2',   5,   True,  0,      50,  None),
                    ('e2',   0.5,   True,  0,      1,  None),
                    ('a3',   2,   True,  0,      None,  None),
                    ('f3',   1882,   True, 1870,    1895,  None),
                    ('l3',   4,   True,  0,      50,  None),
                    ('e3',   0.5,   True,  0,      1,  None))
                        
     
    # Fitting the RT spectrum
    bir = np.array([(1370,1500),(1700,1820),(1920,1950)]) # BIR diamond        
    ycorrRT, baselineRT, coeffsDRT = linbaseline(inputRT[:,0:2],bir,'gcvspline',0.15)
    result = minimize(residual, params,method = "leastsq", args=(ycorrRT[:,0], ycorrRT[:,1])) # fit data with  model from scipy
    youRT, peak1RT, peak2RT, peak3RT = residual(result.params, inputRT[:,0]) # the different peaks
    
    # Fitting the HT spectrum
    ycorrHT, baselineHT, coeffsDHT = linbaseline(inputHT[:,0:2],bir,'gcvspline',0.15)
    resultHT = minimize(residual, paramsHT,method = "leastsq", args=(ycorrHT[:,0], ycorrHT[:,1]))
    youHT, peak1HT, peak2HT, peak3HT = residual(resultHT.params, inputHT[:,0]) # the different peaks
        
    ax1.plot(inputRT[:,0], youRT+baselineRT[:,1],'b-',inputRT[:,0],peak1RT,'g-',inputRT[:,0],peak2RT,'g-',inputRT[:,0],peak3RT,'g-',inputRT[:,0],baselineRT[:,1],'g-')
    ax2.plot(inputHT[:,0], youHT+baselineHT[:,1],'b-',inputHT[:,0],peak1HT,'g-',inputHT[:,0],peak2HT,'g-',inputHT[:,0],peak3HT,'g-',inputHT[:,0],baselineHT[:,1],'g-')
    
    Freq13C_RT = result.params['f1'].value - (result.params['f3'].value - 1883.25) #f1 and f3 are the frequency of 13C and neon line
    Freq13C_HT = resultHT.params['f1'].value - (resultHT.params['f3'].value - 1883.25) #f1 and f3 are the frequency of 13C and neon line
        

    # WE CALCULATE PRESSURE WITH EQ. FROM MYSEN YAMASHITA 2010
    outputs[lg,0] = temperature[lg]

def fit_data(self):
        """Run a scipy leastsq regression."""
        self.parameters = Parameters()
        for f in self.functions:
            for p in f.parameters:
                p.original_name = p.name
                self.parameters[f.name+p.name] = p
                if p.value is None:
                    p.value = 1.0
                p.init_value = p.value
        self.result = minimize(self.residuals, self.parameters)
        if __version__ > '0.8.3':
            for parameter in self.parameters:
                self.parameters[parameter].value = \
                    self.result.params[parameter].value
                self.parameters[parameter].stderr = \
                    self.result.params[parameter].stderr
                self.parameters[parameter].correl = \
                    self.result.params[parameter].correl
        for f in self.functions:
            for p in f.parameters:
                p.name = p.original_name

use_fft = False



    # Set up parameters to fit galaxy with 2-dimensional sersic profile
    params = lmfit.Parameters()
    params.add('sersic_amplitide', value=initial_sersic_amplitide, vary=True)
    params.add('sersic_r_eff', value=initial_sersic_r_eff, vary=True, min=0.0, max=pod['adj_semimaj_pix'])
    params.add('sersic_n', value=initial_sersic_n, vary=True, min=0.1, max=10)
    params.add('sersic_x_0', value=initial_sersic_x_0, vary=False)
    params.add('sersic_y_0', value=initial_sersic_y_0, vary=False)
    params.add('sersic_ellip', value=initial_sersic_ellip, vary=True, min=0.5*initial_sersic_ellip, max=0.5*(1.0-initial_sersic_ellip)+initial_sersic_ellip)
    params.add('sersic_theta', value=initial_sersic_theta, vary=False)

    # Solve with LMfit to find parameters of best-fit sersic profile
    result = lmfit.minimize(Sersic_LMfit, params, args=(pod, cutout, psf, mask, use_fft), method='leastsq', ftol=1E-5, xtol=1E-5, maxfev=200)

    # Extract best-fit results
    sersic_amplitide = result.params['sersic_amplitide'].value
    sersic_r_eff = result.params['sersic_r_eff'].value
    sersic_n = result.params['sersic_n'].value
    sersic_x_0 = result.params['sersic_x_0'].value
    sersic_y_0 = result.params['sersic_y_0'].value
    sersic_ellip = result.params['sersic_ellip'].value
    sersic_theta = result.params['sersic_theta'].value

    # Construct underlying sersic map and convolved sersic map, using best-fit parameters
    sersic_model = astropy.modeling.models.Sersic2D(amplitude=sersic_amplitide, r_eff=sersic_r_eff, n=sersic_n, x_0=sersic_x_0, y_0=sersic_y_0, ellip=sersic_ellip, theta=sersic_theta)
    sersic_map = sersic_model(sersic_x, sersic_y)
    if use_fft==True:
        conv_map = astropy.convolution.convolve_fft(sersic_map, psf, normalize_kernel=True)
    elif use_fft==False:

if (result.errorbars and
                    np.isnan(result.params['Jmax'].stderr) == False):

                    self.succes_count += 1
                else:

                    # Failed errobar fitting, going to try and mess with
                    # starting poisition...
                    for i in xrange(self.Niter):
                        (vcmax_guess, jmax_guess,
                         rd_guess) = self.pick_random_starting_point()

                        params = self.setup_model_params(jmax_guess=jmax_guess,
                                                     vcmax_guess=vcmax_guess,
                                                     rd_guess=rd_guess)
                        result = minimize(self.residual, params,
                                         args=(curve_data, curve_data["Photo"]))
                        if result.errorbars:
                            self.succes_count += 1
                            break

                if print_to_screen:
                    print fname, curve_num
                    self.print_fit_to_screen(result)

                # Need to run the Farquhar model with the fitted params for
                # plotting...
                (An, Anc, Anj) = self.forward_run(result, curve_data)
                self.report_fits(wr, result, os.path.basename(fname),
                                 curve_data, An)

                self.make_plot(curve_data, curve_num, An, Anc, Anj, result)

if print_to_screen:
                    self.print_fit_to_screen(result)
                
                # Did we resolve the error bars during the fit? If yes then
                # move onto the next A-Ci curve
                if result.errorbars:
                    self.succes_count += 1
                
                # Otherwise we did not fit the data. Is it because of our 
                # starting position?? Lets vary this a little and redo the fits
                else:
                    for i in xrange(self.nstartpos):
                        params = self.try_new_params(jmax=True, vcmax=True, 
                                                     rd=True)
                        result = minimize(self.residual, params, 
                                          engine="leastsq", 
                                          args=(curve_data, curve_data["A"]))
                         
                        if print_to_screen:
                            self.print_fit_to_screen(result)
                        if result.errorbars:
                            succes_count += 1
                            break
                # Need to run the Farquhar model with the fitted params for
                # plotting...
                (An, Anc, Anj) = self.forward_run(result, curve_data)
                self.report_fits(wr, result, os.path.basename(fname), 
                                 curve_data, An)
                # Figure out which of our curves (if any) were measured at
                # 25 degrees C
                Tavg = np.mean(curve_data["Tleaf"]) - self.deg2kelvin

for name, mod in fit_models.items():
        fit_variables[name] = sorted(mod.energy.atoms(v.StateVariable).union({v.T, v.P}), key=str)
        # Extra factor '1e-100...' is to work around an annoying broadcasting bug for zero gradient entries
        fit_models[name].models['_broadcaster'] = 1e-100 * sympy.Mul(*fit_variables[name]) ** 3

    callables = {name: make_callable(mod.ast,
                                     itertools.chain(param_names, fit_variables[name]))
                 for name, mod in fit_models.items()}
    grad_callables = {name: make_callable(sympy.Matrix([mod.ast]).jacobian(fit_variables[name]),
                                          itertools.chain(param_names, fit_variables[name]))
                      for name, mod in fit_models.items()}

    #out = leastsq(residual_equilibrium, param_values,
    #              args=(data, dbf, comps, fit_models, callables),
    #              full_output=True, **kwargs)
    out = lmfit.minimize(residual_thermochemical, fit_params,
                         args=(data, dbf, comps, fit_models, callables, grad_callables))
    #fit = residual_equilibrium(data, dbf, comps, fit_models, callables)
    #ssq = np.linalg.norm(residual_equilibrium(out[0], data, dbf, comps,
    #                                          fit_models, callables))
    #return dict(zip(param_names, out[0])), out[1:]
    return fit_params.valuesdict(), out