How to use the lmfit.lineshapes.gaussian function in lmfit

def residual(pars, x, data=None):
        g1 = gaussian(x, pars['a1'], pars['c1'], pars['w1'])
        g2 = gaussian(x, pars['a2'], pars['c2'], pars['w2'])
        model = g1 + g2
        if data is None:
            return model
        return (model - data)

import numpy as np
import pytest

from lmfit.lineshapes import gaussian
from lmfit.models import GaussianModel, LinearModel

try:
    import numdifftools  # noqa: F401
    calc_covar_options = [False, True]
except ImportError:
    calc_covar_options = [False]


np.random.seed(7)
x = np.linspace(0, 20, 401)
y = gaussian(x, amplitude=24.56, center=7.6543, sigma=1.23)
y -= 0.20*x + 3.333 + np.random.normal(scale=0.23, size=len(x))
mod = GaussianModel(prefix='peak_') + LinearModel(prefix='bkg_')

pars = mod.make_params(peak_amplitude=21.0, peak_center=7.0,
                       peak_sigma=2.0, bkg_intercept=2, bkg_slope=0.0)

# set bounds for use with 'differential_evolution' and 'brute'
pars['bkg_intercept'].set(min=0, max=10)
pars['bkg_slope'].set(min=-5, max=5)
pars['peak_amplitude'].set(min=20, max=25)
pars['peak_center'].set(min=5, max=10)
pars['peak_sigma'].set(min=0.5, max=2)


def per_iteration(pars, iter, resid, *args, **kws):
    """Iteration callback, will abort at iteration 23."""

        f2 = lambda x, amp, cen, sig: gaussian(x, amplitude=amp, center=cen,
                                               sigma=sig)
        model2 = Model(f2)

def residual(pars, x, data=None):
        g1 = gaussian(x, pars['a1'], pars['c1'], pars['w1'])
        g2 = gaussian(x, pars['a2'], pars['c2'], pars['w2'])
        model = g1 + g2
        if data is None:
            return model
        return (model - data)

def test_param_set():
    np.random.seed(2015)
    x = np.arange(0, 20, 0.05)
    y = gaussian(x, amplitude=15.43, center=4.5, sigma=2.13)
    y = y + 0.05 - 0.01*x + np.random.normal(scale=0.03, size=len(x))

    model  = VoigtModel()
    params = model.guess(y, x=x)

    # test #1:  gamma is constrained to equal sigma
    assert(params['gamma'].expr == 'sigma')
    params.update_constraints()
    sigval = params['sigma'].value
    assert_allclose(params['gamma'].value, sigval, 1e-4, 1e-4, '', True)

    # test #2: explicitly setting a param value should work, even when
    #          it had been an expression.  The value will be left as fixed
    gamval = 0.87543
    params['gamma'].set(value=gamval)
    assert(params['gamma'].expr is None)

def get_gaussianmodel(amplitude=1.0, center=5.0, sigma=1.0, noise=0.1):
    # create data to be fitted
    np.random.seed(7392)
    x = np.linspace(-20, 20, 201)
    y = gaussian(x, amplitude, center=center, sigma=sigma)
    y = y + np.random.normal(size=len(x), scale=noise)

    model = GaussianModel()
    params = model.make_params(amplitude=amplitude/5.0,
                               center=center-1.0,
                               sigma=sigma*2.0)
    return x, y, model, params

def test_multidatasets():
    # create 5 datasets
    x  = np.linspace( -1, 2, 151)
    data = []
    for i in np.arange(5):
        amp  =  2.60 + 1.50*np.random.rand()
        cen  = -0.20 + 1.50*np.random.rand()
        sig  =  0.25 + 0.03*np.random.rand()
        dat  = gaussian(x, amp, cen, sig) + \
               np.random.normal(size=len(x), scale=0.1)
        data.append(dat)

    # data has shape (5, 151)
    data = np.array(data)
    assert(data.shape) == (5, 151)

    # create 5 sets of parameters, one per data set
    pars = Parameters()
    for iy, y in enumerate(data):
        pars.add( 'amp_%i' % (iy+1), value=0.5, min=0.0,  max=200)
        pars.add( 'cen_%i' % (iy+1), value=0.4, min=-2.0,  max=2.0)
        pars.add( 'sig_%i' % (iy+1), value=0.3, min=0.01, max=3.0)

    # but now constrain all values of sigma to have the same value
    # by assigning sig_2, sig_3, .. sig_5 to be equal to sig_1

def flexible_func(x, amplitude, center, sigma, **kwargs):
            return gaussian(x, amplitude, center, sigma)

def residual(pars, x, sigma=None, data=None):
    yg = gaussian(x, pars['amp_g'], pars['cen_g'], pars['wid_g'])
    yl = lorentzian(x, pars['amp_l'], pars['cen_l'], pars['wid_l'])

    slope = pars['line_slope']
    offset = pars['line_off']
    model = yg + yl + offset + x*slope

    if data is None:
        return model
    if sigma is None:
        return model - data
    return (model - data) / sigma