Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed window function math for tapering function #84

Merged
merged 7 commits into from
May 8, 2018
Merged

Conversation

acliu
Copy link
Contributor

@acliu acliu commented Apr 27, 2018

This fixes the window functions to account for tapering. Currently it assumes a frequency-independent primary beam, but fixing this will require some changes to pyuvdata that haven't been implemented yet.

@acliu acliu requested review from philbull and nkern April 27, 2018 22:43
@nkern
Copy link
Member

nkern commented Apr 27, 2018

the get_G call in PSpecData.pspec should have the taper kwarg propagated, but currently does not get it

@bhazelton
Copy link
Member

Out of curiosity, what do you need from pyuvdata that doesn't currently exist?

@acliu
Copy link
Contributor Author

acliu commented Apr 27, 2018

Right now we have something that returns Omega_pp(\nu_i) ~ \int A^2 (\nu_i) dOmega. Turns out we'll need Omega_pp(\nu_i, \nu_j) ~ \int A(\nu_i) A(\nu_j) dOmega as well.

@bhazelton
Copy link
Member

What does the i & j index? different antennas, frequencies, or something else?

@acliu
Copy link
Contributor Author

acliu commented Apr 27, 2018

Different frequencies.

@coveralls
Copy link

coveralls commented Apr 27, 2018

Coverage Status

Coverage increased (+0.1%) to 92.472% when pulling fe9a8c5 on window_fix into 813145b on master.

@acliu
Copy link
Contributor Author

acliu commented May 3, 2018

Ok, just pushed. Hopefully this time it's a little less buggy.

I mentioned to @nkern today that in working through this, I realized that the pspec scalar normalization infrastructure is not currently able to accommodate non-trivial tapers + non-trivial weighting matrices simultaneously. Details in the notes that I've attached, which will be posted to the HERA memo series.

There is a hacky way to do this, involving assigning half (i.e., the square root) of the weighting matrix to Q and having another square root of the weight act on the data. But this is beginning to feel like epicycles upon epicycles, so if there is a need to push further along this direction I'm ripping out the old stuff and implementing Q "properly" (i.e., I'll subclass and make a version where Q is truly equal to dC/dp).

Notes here: main.pdf

@philbull
Copy link
Collaborator

philbull commented May 7, 2018

@acliu On p6 of your note, should I read an "and" between the bullet points below Eq. 31, or an "or"?

In light of this, should we raise a warning when scalar is used with non-trivial weights R and/or a non-trivial taper?

Copy link
Collaborator

@philbull philbull left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Approved, pending some small changes to a few exception strings.

@@ -681,6 +685,75 @@ def get_G(self, key1, key2):
G[i,j] += np.einsum('ab,ba', iR1Q[i], iR2Q[j])

return G / 2.

def get_H(self, key1, key2, taper='none'):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see H_ab defined anywhere in the note?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, it's there. Equation 33. Unfortunate page break location makes it hard to see.

(See HERA memo #44). As currently implemented, this approximates the
primary beam as frequency independent. Under this approximation, the
our H_ab is defined using the equation above *except* we have
Q^{yet_another} rather than Q_b, where
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yet_another is not a great name for a variable...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changed it to Q^{tapered}

M = np.linalg.pinv(G, rcond=1e-12)
#U,S,V = np.linalg.svd(F)
#M = np.einsum('ij,j,jk', V.T, 1./S, U.T)
raise NotImplementedError
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please say what's not implemented, e.g. NoteImplementedError("G^-1 mode not currently supported.")

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done!

@philbull
Copy link
Collaborator

philbull commented May 7, 2018

(OK, I now see that you do raise warnings about this.)

@acliu acliu merged commit a3040b2 into master May 8, 2018
@acliu acliu deleted the window_fix branch May 8, 2018 02:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants