No more coaddition. Use NN2 flux difference method.

[laldoroty]

I think I should rewrite to use Masao's suggested approach, which is: do single-epoch DIA on each image N times, where N is the number of template images. Then, coadd the difference images at the end. So, I would not coadd anything until the very end. It is more computationally expensive, but it builds on code that already works, and I think I'm wasting a lot of time trying to hunt this down. I could also simultaneously rewrite so that I rotate the reference to match the science instead of the other way around (like I'm currently doing).

To make the science image where I get the ZPT (i.e., the science image that's un-subtracted but has been cross-convolved with the reference and the decorrelation kernel applied), I could just make N of these, do photometry on the stars in all N images, and take the median exactly as I was doing before. This also may allow for SFFT on smaller stamp sizes (1K)--the issue with making the stamps smaller was there weren't enough stars in the stamps to get a ZPT. This would multiply the number of stars by N.

Thoughts?

Originally posted by @laldoroty in https://github.com/laldoroty/phrosty/issues/12#issuecomment-2246186857

Then, @wmwv responded:

If you're going to go done the N-combinations path, I might offer:

"The NN2 Flux Difference Method for Constructing Variable Object Light Curves" https://ui.adsabs.harvard.edu/abs/2005AJ....130.2272B/abstract Barris, Tonry, Novicki, Wood-Vasey """ We present a new method for optimally extracting point-source time variability information from a series of images. Differential photometry is generally best accomplished by subtracting two images separated in time, since this removes all constant objects in the field. By removing background sources such as the host galaxies of supernovae, such subtractions make possible the measurement of the proper flux of point-source objects superposed on extended sources. In traditional difference photometry, a single image is designated as the ``template'' image and is subtracted from all other observations. This procedure does not take all the available information into account and for suboptimal template images may produce poor results. Given N total observations of an object, we show how to obtain an estimate of the vector of fluxes from the individual images using the antisymmetric matrix of flux differences formed from the N(N-1)/2 distinct possible subtractions and provide a prescription for estimating the associated uncertainties. We then demonstrate how this method improves results over the standard procedure of designating one image as a template and differencing against only that image. """

I think going back to single-image subtractions make sense. Coadding the difference images doesn't seem necessary. I would photometry and then combine in catalog space.

So, going to rewrite code to use this method. Rewrote code structure plan: https://github.com/Roman-Supernova-PIT/diff-img/discussions/14

[laldoroty]

Test code, seems to work.

v1 = np.array([3.12])
v2 = np.array([3.22])
v3 = np.array([3.33])
v4 = np.array([3.1])
vecs = np.hstack([v1,v2,v3,v4])

n = 4
N = int(n*(n-1)/2)

print('N, n:')
print(N, n)

dv1 = v1-v2
dv2 = v1-v3
dv3 = v1-v4
dv4 = v2-v3
dv5 = v2-v4
dv6 = v3-v4

A = np.hstack([dv1,dv2,dv3,dv4,dv5,dv6])

print(A)
print('A shape:')
print(A.shape)

C = np.diagflat([0.1,0.12,0.2,0.14,0.16,0.09])

print('C shape:')
print(C.shape)

i, j = np.triu_indices(n, k=1)
X = np.zeros((N,n))

print('X shape:')
print(X.shape)
np.put_along_axis(X, i[:, None], 1, axis=1)
np.put_along_axis(X, j[:, None], -1, axis=1)

XTCXXT = np.linalg.lstsq(X.T @ np.linalg.inv(C) @ X, X.T, rcond=None)
print('XTCXXT shape')
print(XTCXXT[0].shape)

v = XTCXXT[0] @ np.linalg.inv(C) @ A

print('SOLUTION:')
print(vecs - v)

Output:

SOLUTION:
[3.1925 3.1925 3.1925 3.1925]

[laldoroty]

@wmwv when you're back from vacation, can you take a look at the above and confirm the right idea is there? Particularly with the solution being an n-length vector instead of a scalar. I expect that for real data, (vecs-v) will not contain all identical values.

[laldoroty]

I am going to close this issue. I am not working on it and the branch is littered with other stuff. I will re-open a new issue for sorting out the photometry.

[laldoroty]

Merged to new branch associated with issue #20. Closing.