How to use the tractor.emfit.em_fit_2d function in tractor

S = klpsf.shape[0]
        # number of Gaussian components
        K = 3
        w, mu, sig = em_init_params(K, None, None, None)
        II = klpsf.copy()
        II /= II.sum()
        # HIDEOUS HACK
        II = np.maximum(II, 0)
        #print('Multi-Gaussian PSF fit...')
        xm, ym = -(S // 2), -(S // 2)
        if savepsfimg is not None:
            plt.clf()
            plt.imshow(II, interpolation='nearest', origin='lower')
            plt.title('PSF image to fit with EM')
            plt.savefig(savepsfimg)
        res = em_fit_2d(II, xm, ym, w, mu, sig)
        #print('em_fit_2d result:', res)
        if res == 0:
            # print('w,mu,sig', w,mu,sig)
            mypsf = GaussianMixturePSF(w, mu, sig)
            mypsf.computeRadius()
        else:
            # Failed!  Return 'dg' model instead?
            print('PSF model fit', psf, 'failed!  Returning DG model instead')
            psf = 'dg'
    if psf == 'dg':
        print('Creating double-Gaussian PSF approximation')
        (a, s1, b, s2) = dgpsf
        mypsf = NCircularGaussianPSF([s1, s2], [a, b])

    if brightpsf:
        print('Wrapping PSF in SdssBrightPSF')

del I
	del var

	psfimg = pyfits.open(psffn)[0].data
	print('PSF image shape', psfimg.shape)
	# number of Gaussian components
	K = 3
	PS = psfimg.shape[0]
	w,mu,sig = em_init_params(K, None, None, None)
	II = psfimg.copy()
	II /= II.sum()
	# HACK
	II = np.maximum(II, 0)
	print('Multi-Gaussian PSF fit...')
	xm,ym = -(PS/2), -(PS/2)
	em_fit_2d(II, xm, ym, w, mu, sig)
	print('w,mu,sig', w,mu,sig)
	psf = GaussianMixturePSF(w, mu, sig)

	if bandname is None:
		# try looking up in filtermap.
		filt = phdr['FILTER']
		if filt in filtermap:
			print('Mapping filter', filt, 'to', filtermap[filt])
			bandname = filtermap[filt]
		else:
			print('No mapping found for filter', filt)
			bandname = flit

	photocal = cf.CfhtPhotoCal(hdr=phdr, bandname=bandname)

	filename = phdr['FILENAME'].strip()

from tractor.emfit import em_fit_2d
	from tractor.fitpsf import em_init_params

	# Create Gaussian mixture model PSF approximation.
	S = psf.shape[0]
	# number of Gaussian components
	K = 3
	w,mu,sig = em_init_params(K, None, None, None)
	II = psf.copy()
	II /= II.sum()
	# HIDEOUS HACK
	II = np.maximum(II, 0)
	print('Multi-Gaussian PSF fit...')
	xm,ym = -(S/2), -(S/2)
	em_fit_2d(II, xm, ym, w, mu, sig)
	print('w,mu,sig', w,mu,sig)
	mypsf = tractor.GaussianMixturePSF(w, mu, sig)


	P = mypsf.getPointSourcePatch(S/2, S/2)
	mn,mx = psf.min(), psf.max()
	ima = dict(interpolation='nearest', origin='lower',
			   vmin=mn, vmax=mx)
	plt.clf()
	plt.subplot(1,2,1)
	plt.imshow(psf, **ima)
	plt.subplot(1,2,2)
	pimg = np.zeros_like(psf)
	P.addTo(pimg)
	plt.imshow(pimg, **ima)
	plt.savefig('psf.png')

tim = Image(data=image, invvar=invvar, psf=tpsf, wcs=twcs, photocal=photocal,
				sky=tsky, name='CFHT coadd %s %s' % (obj, bandname))

	# set "zr" for plots
	sig = 1./np.median(tim.inverr)
	tim.zr = np.array([-1., +20.]) * sig

	if not doplots:
		return tim

	psfimpatch = Patch(-(PW/2), -(PH/2), psfim)
	# number of Gaussian components
	for K in range(1, 4):
		w,mu,sig = em_init_params(K, None, None, None)
		xm,ym = -(PW/2), -(PH/2)
		em_fit_2d(psfim, xm, ym, w, mu, sig)
		#print 'w,mu,sig', w,mu,sig
		psf = GaussianMixturePSF(w, mu, sig)
		patch = psf.getPointSourcePatch(0, 0)

		plt.clf()
		plt.subplot(1,2,1)
		plt.imshow(patch.getImage(), interpolation='nearest', origin='lower')
		plt.colorbar()
		plt.subplot(1,2,2)
		plt.imshow((patch - psfimpatch).getImage(), interpolation='nearest', origin='lower')
		plt.colorbar()
		plt.savefig('copsf-%i.png' % K)

	plt.clf()
	plt.imshow(psfim, interpolation='nearest', origin='lower')
	plt.colorbar()

# print 'Read PSF image:', psf.shape, 'range', psf.min(), psf.max()
        if positive:
            psf = np.maximum(psf, 0.)
        psf = PixelizedPSF(psf)
        cache[key] = psf
        return psf

    S = psf.shape[0]
    # number of Gaussian components
    K = 3
    w, mu, sig = em_init_params(K, None, None, None)
    II = psf.copy()
    II = np.maximum(II, 0)
    II /= II.sum()
    xm, ym = -(S / 2), -(S / 2)
    res = em_fit_2d(II, xm, ym, w, mu, sig)
    if res != 0:
        raise RuntimeError('Failed to fit PSF')
    print('W1 PSF:')
    print('  w', w)
    print('  mu', mu)
    print('  sigma', sig)
    psf = GaussianMixturePSF(w, mu, sig)
    psf.computeRadius()
    cache[key] = psf
    return psf

invvar=table[1].data
    skyobj = ba.ConstantSky(header['skyval'])
    psffn = 'J082742.02+212844.7-r-bpsf.fits.gz'
    psfimg = pyfits.open(psffn)[0].data
    print('PSF image shape', psfimg.shape)
    # number of Gaussian components
    PS = psfimg.shape[0]
    K = 3
    w,mu,sig = em_init_params(K, None, None, None)
    II = psfimg.copy()
    II /= II.sum()
    # HACK
    II = np.maximum(II, 0)
    print('Multi-Gaussian PSF fit...')
    xm,ym = -(PS/2), -(PS/2)
    em_fit_2d(II, xm, ym, w, mu, sig)
    print('w,mu,sig', w,mu,sig)
    psf = GaussianMixturePSF(w, mu, sig)

    sources = []
#    for run,camcol,field in zip(run,camcol,field):
#        sources.append(st.get_tractor_sources(run,camcol,field,bandname,bands=bands))
    wcs = Tan("J082742.02+212844.7-r.fits",0)
    wcs = FitsWcs(wcs)
    wcs.setX0Y0(1.,1.)
    photocal = ba.NasaSloanPhotoCal(bandname) #Also probably not right
    TI.append(en.Image(data=data,invvar=invvar,sky=skyobj,psf=psf,wcs=wcs,photocal=photocal,name = "NASA-Sloan Test"))
    
    lvl = logging.DEBUG
    logging.basicConfig(level=lvl,format='%(message)s',stream=sys.stdout)
    tims = [TI[0]]
    tractor = st.SDSSTractor(tims)

psfimg = pyfits.open(psffn)[0].data
	print('PSF image shape', psfimg.shape)
	from tractor.emfit import em_fit_2d
	from tractor.fitpsf import em_init_params
	# number of Gaussian components
	S = psfimg.shape[0]
	K = 3
	w,mu,sig = em_init_params(K, None, None, None)
	II = psfimg.copy()
	II /= II.sum()
	# HACK
	II = np.maximum(II, 0)
	print('Multi-Gaussian PSF fit...')
	xm,ym = -(S/2), -(S/2)
	em_fit_2d(II, xm, ym, w, mu, sig)
	print('w,mu,sig', w,mu,sig)
	psf = GaussianMixturePSF(w, mu, sig)

	photocal = cf.CfhtPhotoCal(hdr=phdr,
							   bandname='r')

	cftimg = Image(data=image, invvar=invvar, psf=psf, wcs=wcs,
				   sky=skyobj, photocal=photocal,
				   name='CFHT')
	return cftimg, cfsky, cfstd

plt.subplot(3, 3, 3 * i + j + 1)
            plt.imshow(np.log10(np.maximum(P, 1e-8)),
                       interpolation='nearest', origin='lower', vmax=0.01)
            # plt.colorbar()
    plt.savefig('psf-w1-xy.png')

    psf = pyfits.open('psf%i.fits' % 1)[0].data
    S = psf.shape[0]
    # number of Gaussian components
    for K in range(1, 6):
        w, mu, sig = em_init_params(K, None, None, None)
        II = psf.copy()
        II /= II.sum()
        II = np.maximum(II, 0)
        xm, ym = -(S / 2), -(S / 2)
        res = em_fit_2d(II, xm, ym, w, mu, sig)
        print('em_fit_2d result:', res)
        if res != 0:
            raise RuntimeError('Failed to fit PSF')
        print('w,mu,sig', w, mu, sig)
        mypsf = GaussianMixturePSF(w, mu, sig)
        mypsf.computeRadius()

        #
        mypsf.radius = S / 2
        mod = mypsf.getPointSourcePatch(0., 0.)
        mod = mod.patch
        mod /= mod.sum()

        plt.clf()
        plt.subplot(1, 2, 1)
        ima = dict(interpolation='nearest', origin='lower',