Piff 114 moments

[eacharles]

Probably need to iterate a bit to get full coverage.

[eacharles]

Hmm. The python 2.7 test is failing in installing the dependencies.
And I had to skip the test_gp_interp_anisotropic, which seems to be failing b/c of a change to the treegp code.
Finally, I'm not sure how to get the hsm to return a failure code without actually just raising an error elsewhere, so I don't know if that check is needed.

[rmjarvis]

copy_star = noise_star.copy()

[eacharles]

Star doesn't actually have a copy method. If I add one

def copy(self):
   return copy.deepcopy(self) 

This works fine

[rmjarvis]

Sorry. My mistake. I thought we had one. I'm pretty sure we don't want deepcopy here though. Just copy. We normally want to leave it as shallow as possible.

[eacharles]

That's doesn't work with the current code. It fails unless I use this line in test_moments.py:

noisy_star = copy.deepcopy(noiseless_star)

If I leave that line as a shallow copy the variance of the moments distributions are all zero. Offhand it isn't clear to my why that should fail.

[rmjarvis]

OK. Strange. I thought the line you had here would have been equivalent to a shallow copy. Just copying the pointers to data and fit. I thought that was basically what copy.copy did -- just shallow copy all the attributes.

[rmjarvis]

FYI, I decided to try to close out this issue that's been lingering for over a year now. I did all the things I had suggested in my review from back then, since I think Eric has moved on to other projects.

I also decided to tackle the derivations of the variances, since I'd never been happy with the fudge factors that were in there. So now the math is all included in inline comments, and the analytic variance estimates are all accurate, at least for small e1,e2. (I may try to add the ellipticity terms while the math is all still fresh, but they get pretty complicated for a lot of the moments, so I started with just getting things accurate for roundish stars.)

There is a script in the devel directory called calibrate_moment_errors.py, which I had written as a way to figure out the right fudge factors, but now, with the improved formulae, is rather a script to verify that no fudge factors are required anymore.

[aaronroodman]

I have been swamped by my new LSST responsibilities – but I will look at what you have. If you have an analytic version that includes the terms dependent on the kernel that would be awesome…

…

_________________________________________________________________________ Prof. Aaron Roodman Dept. of Particle Physics & Astrophysics Kavli Institute for Particle Astrophysics & Cosmology SLAC National Accelerator Laboratory Stanford University SLAC National Accelerator Laboratory E-mail: ***@***.******@***.***> 2575 Sand Hill Rd. Phone: 650-926-2705 MS 29 Menlo Park, CA 94025 URL: http://www.slac.stanford.edu/~roodman

_________________________________________________________________________ From: Mike Jarvis ***@***.***> Date: Tuesday, October 26, 2021 at 4:56 PM To: rmjarvis/Piff ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [rmjarvis/Piff] Piff 114 moments (#115) FYI, I decided to try to close out this issue that's been lingering for over a year now. I did all the things I had suggested in my review from back then, since I think Eric has moved on to other projects. I also decided to tackle the derivations of the variances, since I'd never been happy with the fudge factors that were in there. So now the math is all included in inline comments, and the analytic variance estimates are all accurate, at least for small e1,e2. (I may try to add the ellipticity terms while the math is all still fresh, but they get pretty complicated for a lot of the moments, so I started with just getting things accurate for roundish stars.) There is a script in the devel directory called calibrate_moment_errors.py, which I had written as a way to figure out the right fudge factors, but now, with the improved formulae, is rather a script to verify that no fudge factors are required anymore. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#115 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABARV4ZMSQYIX7J6XIBAFZDUI45ZPANCNFSM4OHUL3FA>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[rmjarvis]

Yep, that's the main thing that was required to make things work out. Figuring out dW/dI_k.

Basically, we have

dW/dI_k = dW/du0 du0/dI_k + dW/dv0 dv0/dI_k + dW/dssq dssq/dI_k + dW/de1 de1/dI_k + dW/de2 de2/dI_k.

The first factor in each is fairly straightforward, coming from the Gaussian form of W.
The second factors are trickier, but they come from the constraint equations that HSM imposes.

The math is all there in comments. Some of the algebra is left to the reader, but hopefully it's complete enough to be relatively followable.

[aaronroodman]

Two comments...

1. Should flag=True really be a RunTime error? IIRC this does occur every so often in running, and what I have in my test branch is this instead:

if flag: # just return -1's when HSM fails... logger.info("HSM failed") ret = (-1,-1,-1,-1,-1,-1) if third_order: ret += (-1,-1,-1,-1) if fourth_order: ret += (-1,-1,-1,-1,-1) if radial: ret += (-1,-1,-1,-1,-1,-1) if errors: ret = ret + ret return ret

2. I have added a parallel utility function to return variable names for the moments, given the same arguments. I make use of this to refer to the variables by name in OptAtmo. Something like the code below - although it needs a fix since you've improved the logic for the overlaps between fourth_order and radial. Can you add this too? (sorry, I'm still not competent enough with GitHub to do myself)

`def get_moment_names(third_order=False, fourth_order=False, radial=False, errors=False, kernel=None):
"""Fill vector with moment names, based on arguments to calculate_moments

:return: names         A list of moment variable names
"""

names = []

# 0th, 1st, 2nd order moments are always included
names = ['Flux', 'du', 'dv', 'e0', 'e1', 'e2']

# 3rd order
if third_order:
    names += ['M21', 'M12', 'M30', 'M03']

# 4th order
if fourth_order:
    names += ['M22', 'M31', 'M13', 'M40', 'M04']

# radial
if radial:
    names += ['M22', 'M33', 'M44', 'M22n', 'M33n', 'M44n']

# errors
if errors:

    # TODO: rename these _var instead of _err, since these are actually the variance...

    # 0th, 1st, 2nd order moments are always included
    names += ['Flux_err', 'du_err', 'dv_err', 'e0_err', 'e1_err', 'e2_err']

    # 3rd order
    if third_order:
        names += ['M21_err', 'M12_err', 'M30_err', 'M03_err']

    # 4th order
    if fourth_order:
        names += ['M22_err', 'M31_err','M13_err', 'M40_err', 'M04_err']

    # radial
    if radial:
        names += ['M22_err', 'M33_err', 'M44_err', 'M22n_err', 'M33n_err', 'M44n_err']

return names

`

[rmjarvis]

1. I would think we do want it to raise an error rather than return -1's. All -1 values seems like it would give bad results if used naively. So you're going to want to check for -1s and do something appropriate. The canonical way to do that in python is:

try:
    moments = calculate_moments(star)
except RuntimeError:
    logger.info("calculate_moments failed")
    # .. do something appropriate with the star if we don't have useful moments.
    # i.e. probably skip it in whatever context we're in here.

[rmjarvis]

2. Would it make more sense to have calculate_moments return a dict? Then all the names are directly attached to the value. So you would access moments['M20'] and errors['M20'], rather than trying to keep straight which order the parameters are returned and using the corresponding one.

[rmjarvis]

BTW, in the dict version, you could still reproduce the full list with moments.values(), and the names would be moments.keys() if you really just want them as a list. But if you have a use for the names, it seems like returning a dict that has both the keys and values together would be appropriate.

[aaronroodman]

Returning a dict would be ok. I wrote some code for OptAtmo to parse the moment names based on the yaml configuration, to pull out just the values used in the Chi2, and could modify that to use the keys and items of a dict. For the moment-based fitting this gets called a lot, but I guess there isn't that much overhead in forming a small dictionary each time.

As to the error checking, I would have to change my code that calls this to adjust to the error condition, as would any other user of it. Not a big modification of course, and you are right that it ensures that the failure is noticed by the calling code. I already have a method that calls this given a list of stars, and can put the checking there. So thats ok with me too.

[rmjarvis]

OK, I switched it to return a dict for the moments, and if requested a second dict with the same keys for the variances. I think that's a nicer API than returning a long opaque list of values. The values in the dict are ordered, so you can turn it back to a long list if you want, with the same ordering that we currently have. But I suspect user code will be significantly more legible using the key names for access.

[rmjarvis]

I also looked into getting the ellipticity terms correct, but they quickly got out of hand as lots of high order moments became relevant (M42, M24, M04, M40, etc.). So I think we should be content that the variances are pretty good for round stars, but maybe not so great when the ellipticity is more than a few percent.