How to use the tractor.mixture_profiles.MixtureOfGaussians function in tractor

x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

mg = mp.MixtureOfGaussians([1.], [0.,1.], np.array([2.]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

mg = mp.MixtureOfGaussians([1.], [0.,1.], np.array([ [[1.3, 0.1],[0.1,3.1]], ]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

mg = mp.MixtureOfGaussians([1., 0.5], np.array([ [0.,1.], [-0.3,2] ]),
						   np.array([ [[1.3, 0.1],[0.1,3.1]], [[1.2, -0.8],[-0.8, 2.4]], ]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

# via test_hogg_galaxy.py

mg = mp.MixtureOfGaussians(
	np.array([  4.31155865e-05,   1.34300460e-03,   1.62488556e-02,
				1.13537806e-01,   4.19327122e-01,   4.49500096e-01]),
	np.array([[ 50., 66.],
			  [ 50., 66.],
			  [ 50., 66.],
			  [ 50., 66.],
			  [ 50., 66.],

mg = mp.MixtureOfGaussians([1.], [0.,0.], np.array([1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

mg = mp.MixtureOfGaussians([1.], [0.,1.], np.array([2.]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

mg = mp.MixtureOfGaussians([1.], [0.,1.], np.array([ [[1.3, 0.1],[0.1,3.1]], ]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

mg = mp.MixtureOfGaussians([1., 0.5], np.array([ [0.,1.], [-0.3,2] ]),
						   np.array([ [[1.3, 0.1],[0.1,3.1]], [[1.2, -0.8],[-0.8, 2.4]], ]))
x = mg.evaluate(np.array([0.,1.]))
x = mg.evaluate(np.array([0,0]))
x = mg.evaluate(np.array([-3.,2.]))
x = mg.evaluate(np.array([[-3.,2.], [-17,4], [4,-2]]))

# via test_hogg_galaxy.py

mg = mp.MixtureOfGaussians(
	np.array([  4.31155865e-05,   1.34300460e-03,   1.62488556e-02,

amix = self._getShearedProfile(img, px, py)
        fftmix = amix
        mogmix = None

        #print('Galaxy affine profile:', amix)

        if hybrid:
            # Split "amix" into terms that we will evaluate using MoG
            # vs FFT.
            vv = amix.var[:, 0, 0] + amix.var[:, 1, 1]
            nsigma = 3.
            # Terms that will wrap-around significantly...
            I = ((vv * nsigma**2) > (pW**2))
            if np.any(I):
                # print('Evaluating', np.sum(I), 'terms as MoGs')
                mogmix = mp.MixtureOfGaussians(amix.amp[I],
                                               amix.mean[I, :] + np.array([px, py])[np.newaxis, :],
                                               amix.var[I, :, :], quick=True)
            I = np.logical_not(I)
            if np.any(I):
                # print('Evaluating', np.sum(I), 'terms with FFT')
                fftmix = mp.MixtureOfGaussians(amix.amp[I], amix.mean[I, :], amix.var[I, :, :],
                                                quick=True)
            else:
                fftmix = None

        if fftmix is not None:
            Fsum = fftmix.getFourierTransform(v, w, zero_mean=True)
            # In Intel's mkl_fft library, the irfftn code path is faster than irfft2
            # (the irfft2 version sets args (to their default values) which triggers padding
            #  behavior, changing the FFT size and copying behavior)
            #G = np.fft.irfft2(Fsum * P, s=(pH, pW))

def getCircularMog(amps, sigmas):
    K = len(amps)
    amps = np.array(amps).astype(np.float32)
    means = np.zeros((K, 2), np.float32)
    vars = np.zeros((K, 2, 2), np.float32)
    for k in range(K):
        vars[k, 0, 0] = vars[k, 1, 1] = sigmas[k]**2
    return mp.MixtureOfGaussians(amps, means, vars, quick=True)

def apply_affine(self, shift, scale):
        '''
        NOTE, "scale" is transposed vs earlier versions of this code!

        NOTE, does NOT make a copy of amplitude!!

        shift: D-vector offset
        scale: DxD-matrix transformation
        '''
        assert(shift.shape == (self.D,))
        assert(scale.shape == (self.D, self.D))
        newmean = self.mean + shift
        newvar = np.zeros_like(self.var)
        for k in range(self.K):
            newvar[k, :, :] = np.linalg.multi_dot((scale, self.var[k,:,:], scale.T))
        return MixtureOfGaussians(self.amp, newmean, newvar, quick=True)

if hybrid:
            # Split "amix" into terms that we will evaluate using MoG
            # vs FFT.
            vv = amix.var[:, 0, 0] + amix.var[:, 1, 1]
            nsigma = 3.
            # Terms that will wrap-around significantly...
            I = ((vv * nsigma**2) > (pW**2))
            if np.any(I):
                # print('Evaluating', np.sum(I), 'terms as MoGs')
                mogmix = mp.MixtureOfGaussians(amix.amp[I],
                                               amix.mean[I, :] + np.array([px, py])[np.newaxis, :],
                                               amix.var[I, :, :], quick=True)
            I = np.logical_not(I)
            if np.any(I):
                # print('Evaluating', np.sum(I), 'terms with FFT')
                fftmix = mp.MixtureOfGaussians(amix.amp[I], amix.mean[I, :], amix.var[I, :, :],
                                                quick=True)
            else:
                fftmix = None

        if fftmix is not None:
            Fsum = fftmix.getFourierTransform(v, w, zero_mean=True)
            # In Intel's mkl_fft library, the irfftn code path is faster than irfft2
            # (the irfft2 version sets args (to their default values) which triggers padding
            #  behavior, changing the FFT size and copying behavior)
            #G = np.fft.irfft2(Fsum * P, s=(pH, pW))
            G = np.fft.irfftn(Fsum * P)

            assert(G.shape == (pH,pW))
            # FIXME -- we could try to be sneaky and Lanczos-interp
            # after cutting G down to nearly its final size... tricky
            # tho

def get_dev_mixture():
    return MixtureOfGaussians(dev_amp, np.zeros((dev_amp.size, 2)), dev_var)

varr = np.exp(np.array([f(sindex) for f in v0 + v1]))
        else:
            assert(len(matches) == 1)
            lo,hi,amp_funcs,logvar_funcs = matches[0]
            amps = np.array([f(sindex) for f in amp_funcs])
            amps /= amps.sum()
            varr = np.exp(np.array([f(sindex) for f in logvar_funcs]))

        # Core
        if sindex > 1.:
            core = self.core_func(sindex)
            amps *= (1. - core) / amps.sum()
            amps = np.append(amps, core)
            varr = np.append(varr, 0.)

        return mp.MixtureOfGaussians(amps, np.zeros((len(amps), 2)), varr)

def copy(self):
        return MixtureOfGaussians(self.amp, self.mean, self.var, quick=True)

S /= np.sqrt(np.linalg.det(S))
    print('Det', np.linalg.det(S))
    r = np.deg2rad(30.)
    cr = np.cos(r)
    sr = np.sin(r)
    S = np.dot(S, np.array([[cr, sr], [-sr, cr]]))
    print('Det', np.linalg.det(S))
    exp_mixture = exp_mixture.apply_affine(np.array([10, -30]), S)

    if psf is not None:
        exp_mixture = exp_mixture.convolve(psf)
    pa = [int(x) for x in [posmin[0], posmax[0] + 1, posmin[1], posmax[1] + 1]]
    exp_mix_patch = mixture_to_patch(exp_mixture, *pa)
    exp_mix_patch = exp_mix_patch.patch
    exp_patch = model_to_patch('exp', scale, posmin, posmax)
    dev_mixture = MixtureOfGaussians(
        dev_amp * scale * scale, np.zeros((dev_amp.size, 2)), dev_var * scale * scale)
    if psf is not None:
        dev_mixture = dev_mixture.convolve(psf)
    dev_mix_patch = mixture_to_patch(dev_mixture, *pa)
    dev_mix_patch = dev_mix_patch.patch
    dev_patch = model_to_patch('dev', scale, posmin, posmax)
    cmap = cm.gray
    vmin = -0.5
    vmax = 1.0
    factor = 100.
    plt.clf()
    plt.subplot(231)
    plt.imshow(exp_mix_patch, interpolation='nearest',
               origin='lower', cmap=cmap, vmin=vmin, vmax=vmax)
    plt.colorbar()
    plt.subplot(232)