Rationalize/document error simulation functions

[mjuric]

A discussion with @rhiannonlynne and @ivezic revealed that we have a number of (each slightly different!) functions for computation of astrometric precision, SNR, photometric precision, etc. Here are just a few:

• https://github.com/lsst/sims_photUtils/blob/master/python/lsst/sims/photUtils/SignalToNoise.py
• https://github.com/lsst/sims_maf/blob/master/python/lsst/sims/maf/utils/astrometryUtils.py
• https://github.com/lsst/mops_daymops/blob/99af0c2817f9fd55e2d6187d0ff2f0b271847b6d/python/add_astrometric_noise.py

We should:

• Document that the functions to be used are in https://github.com/lsst/sims_photUtils/blob/master/python/lsst/sims/photUtils/SignalToNoise.py and others should be ignored
• Submit PRs to mark all other functions as "obsolete" and point users to the above.
• @ivezic needs to confirm whether the formulae in https://github.com/lsst/sims_photUtils/blob/557295365d7cbb4ba14068f51719138413954264/python/lsst/sims/photUtils/SignalToNoise.py#L502 are total or per-coordinate. That should then be documented.
• We should extend calcAstrometricError function to use FWHM (and document that the FWHM to use is the geometric FWHM)

Pinging @eggls6 to keep him in the loop (and discuss when we next meet on how/what we use for simulating the solar system catalog).

[ivezic]

Here is my contribution in the form of revised (and tested!) code.
For theoretical background, see comments in the code
calcRandomAstrometricErrorPerCoord. (I am having formatting
issues with the code - ask me directly for the code if need be).

def calcAstrometricError(mag, m5, nvisit=1, FWHMeff=700.0, error_sys = 10.0, astErrCoeff=0.60):
    """
    Calculate the astrometric error, for object catalog purposes.
    The effective FWHMeff MUST BE given in miliarcsec (NOT arcsec!).
    Systematic error, error_sys, must be given in miliarcsec. 
    Returns astrometric error in miliarcsec.
    *** This error corresponds to a single-coordinate error ***
        the total astrometric error (e.g. relevant when 
        matching two catalogs) will be sqrt(2) times larger! 
    """
    # The astrometric error can be applied to parallax or proper motion (for nvisit>1).
    # If applying to proper motion, should also divide by the # of years of the survey.
    # This is also referenced in the LSST overview paper (arXiv:0805.2366, ls.st/lop)
    # 
    # Notes:
    #  - assumes sqrt(Nvisit) scaling, which is the best-case scenario
    #  - calcRandomAstrometricError assumes maxiumm likelihood solution, 
    #     which is also the best-case scenario
    #  - the systematic error, error_sys = 10 mas, corresponds to the 
    #     design spec from the LSST Science Requirements Document (ls.st/srd)
    # 
    # first compute SNR 
    rgamma = 0.039
    xval = numpy.power(10, 0.4*(mag-m5))
    SNR = 1/numpy.sqrt((0.04-rgamma)*xval + rgamma*xval*xval)
    # random astrometric error for a single visit 
    error_rand = calcRandomAstrometricErrorPerCoord(FWHMeff, SNR, astErrCoeff)
    # random astrometric error for nvisit observations
    error_rand = error_rand / numpy.sqrt(nvisit)
    # add systematic error floor:
    astrom_error = numpy.sqrt(error_sys * error_sys + error_rand*error_rand)
    return astrom_error, SNR, error_rand 


def calcRandomAstrometricErrorPerCoord(FWHMeff, SNR, AstromErrCoeff = 0.60):
    """
    Calculate the random astrometric error, as a function of 
    effective FWHMeff and signal-to-noise ratio SNR
    Returns astrometric error in the same units as FWHM.
    ** This error corresponds to a single-coordinate error **
    the total astrometric error (e.g. relevant when matching
    two catalogs) will be sqrt(2) times larger. 
    
    The coefficient AstromErrCoeff for Maximum Likelihood
    solution is given by 
       AstromErrCoeff = <P^2> / <|dP/dx|^2> * 1/FWHMeff
    where P is the point spread function, P(x,y). 
    
    For a single-Gaussian PSF, AstromErrCoeff = 0.60
    For a double-Gaussian approximation to Kolmogorov 
    seeing, AstromErrCoeff = 0.55; however, given the 
    same core seeing (FWHMgeom) as for a single-Gaussian
    PSF, the resulting error will be 36% larger because
    FWHMeff is 1.22 times larger and SNR is 1.22 times
    smaller, compared to error for single-Gaussian PSF.
    Although Kolmogorov seeing is a much better approximation
    of the free atmospheric seeing than single Gaussian seeing, 
    the default value of AstromErrCoeff is set to the
    more conservative value.
    Note also that AstromErrCoeff = 1.0 is often used in
    practice to empirically account for other error sources.
    
    The <June2020 version of this code used AstromErrCoeff=1.0
    and FWHMgeom instead of FWHMeff. Given that 
    FWHMeff/FWHMgeom=1.22 for Kolmogorov seeing, the old 
    version was within 3% from the total astrometric error 
    (that is, not error per coordinate, but sqrt(2) larger)
    produced with this updated version.
    """
    return AstromErrCoeff * FWHMeff / SNR

[ivezic]

I agree with Mario that there should be only one instance of this code. Indeed, thanks to some bad memories of similar issues, I would suggest to consider intentionally braking all other instances (e.g. by returning errors and pointing to the "official" single instance). Nevertheless, I leave this decision to the "professionals". Please let me know if you need anything else from me.

[rhiannonlynne]

Quick question Zeljko -- did you mean FWHMeff or FWHMgeom? I would have thought this would relate to FWHMgeom, not FWHMeff.

…

On Sun, Jun 7, 2020 at 2:47 PM Zeljko Ivezic ***@***.***> wrote: I agree with Mario that there should be only one instance of this code. Indeed, thanks to some bad memories of similar issues, I would suggest to consider intentionally braking all other instances (e.g. by returning errors and pointing to the "official" single instance). Nevertheless, I leave this decision to the "professionals". Please let me know if you need anything else from me. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#191 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAOFBX6X2Q2BBPMVCC3Z3EDRVQDGLANCNFSM4NRAPYKQ> .

-- ========================= Dr R. Lynne Jones Research Scientist @ University of Washington / Dirac Institute / Rubin Observatory she/hers =========================

[rhiannonlynne]

Although .. now I properly page through your notes, I see you explicitly referred to using FWHMeff and have some calculations about what the difference should be -- is the change based on calculations you did or is there a reference somewhere else?

…

On Mon, Jun 8, 2020 at 12:51 PM Lynne Jones ***@***.***> wrote: Quick question Zeljko -- did you mean FWHMeff or FWHMgeom? I would have thought this would relate to FWHMgeom, not FWHMeff. On Sun, Jun 7, 2020 at 2:47 PM Zeljko Ivezic ***@***.***> wrote: > I agree with Mario that there should be only one instance of this code. > Indeed, thanks to some bad memories of similar issues, I would suggest to > consider intentionally braking all other instances (e.g. by returning > errors and pointing to the "official" single instance). Nevertheless, I > leave this decision to the "professionals". Please let me know if you need > anything else from me. > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#191 (comment)>, or > unsubscribe > <https://github.com/notifications/unsubscribe-auth/AAOFBX6X2Q2BBPMVCC3Z3EDRVQDGLANCNFSM4NRAPYKQ> > . > -- ========================= Dr R. Lynne Jones Research Scientist @ University of Washington / Dirac Institute / Rubin Observatory she/hers =========================

-- ========================= Dr R. Lynne Jones Research Scientist @ University of Washington / Dirac Institute / Rubin Observatory she/hers =========================

[ivezic]

It is true that astrometric error is controlled but the PSF core width rather than its wings. However, for a given profile, FWHMeff and FWHMgeom are monotonically related. One can calibrate the relation using either one and I chose FWHMeff because it's needed for SNR too, while FWHMgeom might be unavailable. The scaling factor can be easily obtained numerically from the ML-derived expression and this is how I estimated it.