How to use the tractor.fitpsf.em_init_params function in tractor

np.all(klpsf[:, -1] == 0.)):
                klpsf = klpsf[1:-1, 1:-1]
            else:
                break
        mypsf = PixelizedPSF(klpsf)

    elif psf == 'kl-gm':
        from tractor.emfit import em_fit_2d
        from tractor.fitpsf import em_init_params

        # Create Gaussian mixture model PSF approximation.
        klpsf = psfield.getPsfAtPoints(bandnum, x0 + W / 2, y0 + H / 2)
        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)

v3=False):
        '''
        optional P0 = (w,mu,var): initial parameter guess.

        w has shape (N,)
        mu has shape (N,2)
        var (variance) has shape (N,2,2)

        optional xy0 = int x0,y0 origin of stamp.
        '''
        from tractor.emfit import em_fit_2d_reg
        from tractor.fitpsf import em_init_params
        if P0 is not None:
            w, mu, var = P0
        else:
            w, mu, var = em_init_params(N, None, None, None)
        stamp = stamp.copy()

        if xy0 is None:
            xm, ym = -(stamp.shape[1] // 2), -(stamp.shape[0] // 2)
        else:
            xm, ym = xy0

        if v3:
            tpsf = GaussianMixturePSF(w, mu, var)
            tim = Image(data=stamp, invvar=1e6 * np.ones_like(stamp),
                        psf=tpsf)
            h, w = tim.shape
            src = PointSource(PixPos(w // 2, h // 2), Flux(1.))
            tr = Tractor([tim], [src])
            tr.freezeParam('catalog')
            tim.freezeAllBut('psf')

N = 9
	assert(((H % N) == 0) and ((W % N) == 0))
	# Select which of the NxN PSF images applies to our cutout.
	ix = int(N * float(xc) / OW)
	iy = int(N * float(yc) / OH)
	print('PSF image number', ix,iy)
	PW,PH = W/N, H/N
	print('PSF image shape', PW,PH)

	psfim = psfim[iy*PH: (iy+1)*PH, ix*PW: (ix+1)*PW]
	print('my PSF image shape', PW,PH)
	psfim = np.maximum(psfim, 0)
	psfim /= np.sum(psfim)

	K = psfK
	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)
	tpsf = GaussianMixturePSF(w, mu, sig)

	tsky = ConstantSky(0.)

	obj = phdr['OBJECT'].strip()

	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:

plt.clf()
    for i, y in enumerate([0, 500, 1000]):
        for j, x in enumerate([0, 500, 1000]):
            P = pyfits.open('psf-1-%i-%i.fits' % (x, y))[0].data
            P /= P.max()
            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

obj = phdr['OBJECT'].strip()

	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()

hdr = P[2].header
	badbits = [hdr.get('MP_%s' % nm) for nm in ['BAD', 'SAT', 'INTRP', 'CR']]
	print('Bad bits:', badbits)
	badmask = sum([1 << bit for bit in badbits])
	#print 'Bad mask:', badmask
	#print 'Mask dtype', mask.dtype
	invvar[(mask & int(badmask)) > 0] = 0.
	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]

cas.cut(np.logical_or(Ibest, Iprim))
    print('Cut to', len(cas), 'PRIMARY + BEST-not-near-PRIMARY')

    I, J, d = match_radec(cas.ra, cas.dec, cas.ra,
                          cas.dec, 2. / 3600., notself=True)
    plt.clf()
    loghist((cas.ra[I] - cas.ra[J]) * 3600.,
            (cas.dec[I] - cas.dec[J]) * 3600., 200)
    plt.title('CAS self-matches')
    ps.savefig()

    psf = pyfits.open('wise-psf-w1-500-500.fits')[0].data
    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)
    w1psf = GaussianMixturePSF(w, mu, sig)
    w1psf.computeRadius()

    print('PSF radius:', w1psf.getRadius(), 'pixels')

assert(all(np.isfinite(img.ravel())))
	assert(all(np.isfinite(invvar.ravel())))

	psf = pyfits.open(psffn)[0].data
	print('psf', psf.shape)
	psf /= psf.sum()

	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()