search

dstndstn / tractor

The Tractor: measuring astronomical sources via probabilistic inference
Other
86 stars 24 forks source link

Support for pixel-to-pixel covariance? #85

Closed dlakaplan closed 6 years ago

dlakaplan commented 6 years ago

Is there support for pixel-to-pixel covariance (such as in radio synthesis imaging, or even in drizzled HST images) in tractor? It looks like the interface to Image() supports inverse error/variance maps, but no covariance between pixels. Or can the Image.psf implement covariance?

If I want to implement that myself (say by subclassing Image()) can you give any pointers on what needs changing? Is there some modification of getChiImage() that would be sufficient?

dstndstn commented 6 years ago

There's no existing support for that. Are you thinking about changing the error model: rather than having pixelwise-independent errors, you want there to be a correlation matrix between pixels? If so, then I think you want to look in tractor/engine.py for getLogLikelihood(), which currently calls getChiImages() and then just takes the sum of pixelwise chi-squared values. You instead want to call getChiImage() and push it through the covariance matrix.

That will work for making the getLogProb() function work. The other part of the puzzle (if you're using the existing optimizer rather than doing MCMC) is to get the derivatives right. The code in lsqr_optimizer.py is brutal, but it does have one or two useful comments, including this on the matrix equation being solved:

    # We want to minimize:
    #   || chi + (d(chi)/d(params)) * dparams ||^2
    # So  b = chi
    #     A = -d(chi)/d(params)
    #     x = dparams
    #
    # chi = (data - model) / std = (data - model) * inverr
    # derivs = d(model)/d(param)
    # A matrix = -d(chi)/d(param)
    #          = + (derivs) * inverr

changing that to handle correlated errors will I guess require replacing the thing to be minimized from chi^2 to your correlation matrix equation. If that matrix is compact then maybe there's an easy-ish modification of it...

(And as an aside, drizzle is totally unjustified and shouldn't be used for science; one should model the individual images)

dlakaplan commented 6 years ago

Thank you. Generally agree about drizzle, but I was just trying to make the request for this a bit more general than just radio astronomy.

The covariance matrix is pretty compact. So maybe this is workable, or we could use MCMC as you suggest.

Thanks,

David