WIP: Support spectra & noise propagation #32

jehturner · 2017-09-21T20:39:25Z

This is the same PR as cmccully/lacosmicx#1, rebased on astroscrappy, so we can close that old one. It doesn't yet implement the intended behaviour summarized in that PR, so do not merge yet.

jehturner · 2017-09-21T20:53:05Z

Regarding this:

You'd like to get rid of pssl and my bgsub and replace both of them with bkg, which can accept either a fixed background level (like pssl) or an array of model values to be subtracted (like bgsub, mainly for spectroscopy).

There's actually no reason to support subtracting a constant background internally, is there? You were probably just proposing to get rid of pssl altogether for imaging and require the user to supply indat with its background still included, unless also passing a noise image.

coveralls · 2017-09-21T21:39:12Z

Coverage remained the same at 100.0% when pulling 110f162 on jehturner:spectroscopy into acaa02e on astropy:master.

coveralls · 2017-09-23T00:27:22Z

Coverage remained the same at 100.0% when pulling d628fc8 on jehturner:spectroscopy into 750cc76 on astropy:master.

jehturner · 2017-09-23T00:37:58Z

Does this look like what you had in mind, @cmccully? I have not added a noise argument yet, but this removes the pssl argument -- instead expecting indat to have its sky background included -- and renames my bgsub to bkg, an optional array that gets subtracted temporarily from the data during detection (I don't see an obvious reason to subtract a constant for imaging but one could).

jehturner · 2018-04-16T19:02:41Z

For completeness, I've had a go at adding a var parameter that accepts a user-supplied variance array, though perhaps it should go its own PR. Anyway, on reflection, it may be less useful than I had been thinking, given the following considerations:

If bkg contains all the signal that has been subtracted prior to running astroscrappy then there isn't really any additional information about the noise in var.
Like the main data array, I think var needs its own median filtering and cleaning, to determine what the noise would be in the absence of CRs (which is obviously a bit more processing).
If (unlike the main data array) var is filtered without subtracting & restoring bkg, I think that can lead to less accurate noise structure by blurring out fine background details (such as sky lines or IFU fibres).
In the case I've tested so far, there was very little difference in the cleaned result between using var or just indat, with only ~0.01% of pixels having significant differences. Those cleaned pixels appear similar but a little noisier when using var than they do with the original noise algorithm.

jehturner · 2018-04-16T21:16:18Z

I'll push what I did so you can see it. We can always revert that commit if needed...

I'm not returning a noise image here, but I suppose cleanvar could be added to the return tuple? I tend to re-clean my data after running LA Cosmic anyway, since I find it does a better job of identifying the CRs than removing them (eg. compared with local interpolation), so crmask seems like the most useful return value.

cmccully

This overall, this looks really helpful. I made a few minor comments about the implementation. I'm going to merge #41 before this so you can rebase to that if you would like.

cmccully · 2018-12-18T18:21:49Z

astroscrappy/astroscrappy.pyx

+        Input data array that will be used for cosmic ray detection. This
+        should include the sky background (or a mean background level, added
+        back in after sky subtraction), so that noise can be estimated
+        correctly from the data values.


maybe say unless the var array is also provided.

cmccully · 2018-12-18T18:22:48Z

astroscrappy/astroscrappy.pyx

+
+    var : float numpy array, optional
+        A pre-determined estimate of the data variance (ie. noise squared) in
+        each pixel, generated by previous processing of ``indat``. If provided,


I might remove "generated by previous processing of indat".

cmccully · 2018-12-18T18:23:36Z

astroscrappy/astroscrappy.pyx

+        each pixel, generated by previous processing of ``indat``. If provided,
+        this is used in place of an internal noise model based on ``indat``,
+        ``gain`` and ``readnoise``. This still gets median filtered and cleaned
+        internally, to estimate what the noise in each pixel *would* be in the


Is the median filtering done in place? We probably shouldn't. We should probably make a copy of the array if we are not already.

cmccully · 2018-12-18T18:25:46Z

astroscrappy/astroscrappy.pyx

+        internally, to estimate what the noise in each pixel *would* be in the
+        absence of cosmic rays, but without removing ``bkg`` temporarily (a
+        difference that could lead to less accurate results than the default
+        noise model close to fine structure). This argument should be provided


Would it just make sense to take the masked median of the variance array? Or just replace the masked regions with the median filter of the var array? This way you wouldn't have to worry about the bkg stuff.

cmccully · 2018-12-18T18:26:19Z

astroscrappy/astroscrappy.pyx

+        noise model close to fine structure). This argument should be provided
+        for correct results if ``bkg`` does not include all signal that has
+        previously been subtracted from ``indat`` (other than electronic bias).
+


You might want to note that var is expected to be in adu like bkg and indat (I really should have named that variable better).

cmccully · 2018-12-18T18:30:11Z

astroscrappy/astroscrappy.pyx

+    if var is not None:
+        goodvar = np.empty_like(gooddata, order='c')
+        goodvar[:] = var[np.logical_not(mask)]
+        var_level = median(goodvar, len(goodvar))


do you want an else var_level = background_level?

cmccully · 2018-12-18T18:32:49Z

astroscrappy/astroscrappy.pyx

+            # order to estimate the noise accurately (first saving a copy if
+            # using the same median to replace bad values):
+            if cleantype == 'median':
+                m5_nobkg = m5 if bkg is None else m5.copy()


Can you just replace this by working on the cleanvar below? See comment below.

cmccully · 2018-12-18T18:34:04Z

astroscrappy/astroscrappy.pyx

-
-        if cleantype != 'median':
+        if var is None:
+            noise = np.sqrt(m5 + readnoise * readnoise)


Instead of using m5 in the noise, why don't you use the 5x5 median of cleanvar. Basically calculate one m5 for the cleanarr for the cr detection and one m5 for the noise calculation. Then you don't have to worry about adding back in bkg etc.

cmccully · 2018-12-18T18:34:59Z

astroscrappy/astroscrappy.pyx

-            del m5
+            cleanarr[crinds] = m5_nobkg[crinds]
+            del m5_nobkg
+            if var is not None:


Similarly here, I think you can have a noise and just median of the cleanarr which would simplify this whole block considerably.

griffin-h · 2019-01-11T18:17:13Z

astroscrappy/astroscrappy.pyx

+    # Subtract the input sky model, if applicable.
+    if bkg is not None:
+        cleanarr -= bkg
+


I think this should happen much later in the function. The background needs to be included when you set cleanvar = cleanarr at line 221.

The background for the noise estimate gets restored again at L303, which needs to happen after the median filtering; see recent discussion at jehturner#1 for further explanation.

jehturner · 2019-03-06T15:07:25Z

Rebased as requested.

jehturner · 2019-03-06T15:40:04Z

Not sure what is going on with the first, Python 2.7 build test? It seems to install Python 3.7 with Conda and then fail, but I don't think it's relevant to this PR.

jehturner · 2019-03-08T00:47:53Z

Before answering your detailed suggestions (thanks), I'd like to point out that I was leaning towards reverting my last commit (see comment of 16 April). Input variance was something we discussed on Slack that seemed like a good idea at the time, but once I actually implemented var, I realized that it doesn't work as well as the original LACosmic noise estimate, for 2 reasons: First, var has no equivalent of the bkg image for indat, so continuum/sky structure gets left in the variance when filtering for cosmic rays, producing noise estimates that are less, rather than more, accurate. Second, bkg should help ensure that all the signal gets accounted for anyway when producing LACosmic's usual noise estimate, because it includes any "missing" flux that a separate variance array would otherwise tell you about. So I would probably just get rid of var to avoid overcomplicating things, but obviously if you want to keep it then we can. Admittedly, I'm mostly interested in spectroscopy, where bkg is going to be essential, whereas for imaging bkg is more likely to be a scalar and then you might get slightly more accurate noise structure using a var image as well (off the top of my head). Does that make sense?

jehturner · 2020-07-04T02:45:01Z

In our meeting earlier this week, @cmccully mentioned that he's keen to retain variance input as an option because he has a use case (I think for echelle data) where prior processing produces a non-trivial correspondence betweeen indat & var. I've therefore had another look at what's needed to address the limitations I mentioned above on 16 April 2018, but I've had quite a hard time remembering what I was doing here after 3 years...

The good news is that I think the code here already does most of the duplicate processing necessary to use input variance accurately and is "only" missing subtraction/restoration of the background structure (as for indat).

The problem is simply how to go about subtracting the background from the variance. Unless the only processing done prior to cosmic ray removal is bias subtraction and maybe stacking (in which case you don't really need var, because the original algorithm can derive the same thing from indat), the levels of background structure in var and indat may differ, even if just by a scaling factor or similar. AstroScrappy doesn't know anything about that, so the caller would be required to provide, say, bkg and varbkg separately, which is already getting a bit messy. Moreover, at present the run time for spectroscopy is already dominated by background fitting rather than by detect_cosmics itself, so if the calling function (equivalent to lacos_spec.cl) has to fit the sky+continuum for both indat & var separately, that's quite a large overhead. Alternatively, the caller might know enough about previous processing to convert bkg into varbkg directly, but that might sacrifice generality and/or produce a convoluted API -- I think it would be good to understand your use case a bit better here. And if the caller doesn't get this exactly right, things are liable to go (at least slightly) wrong without it being obvious.

So what are your thoughts on how exactly var would be used? I think I'm fairly clear on how detect_cosmics would need to be modified if we know how we want it to look.

cmccully · 2020-07-31T20:24:55Z

Can you rebase this into master @jehturner ?

jehturner · 2020-08-04T22:10:39Z

Rebase from master? OK.

…rappy's predecessor, lacosmicx (producing identical results in a quick test).

…d background internally instead of taking pre-subtracted input, as discussed with cmucully on the astropyspectroscopy slack channel in April.

…noise estimates, in place of the internal noise model. (May want some tweaking, eg. to protect bkg structure for noise estimates too??)

jehturner · 2020-08-04T22:26:11Z

That was suspiciously easy... Anyway, the diff is just 1 file now.

codecov · 2020-08-04T22:30:57Z

Codecov Report

Merging #32 into master will decrease coverage by 1.93%.
The diff coverage is 47.50%.

@@            Coverage Diff             @@
##           master      #32      +/-   ##
==========================================
- Coverage   94.88%   92.94%   -1.94%     
==========================================
  Files           7        7              
  Lines        1016     1049      +33     
  Branches       53       53              
==========================================
+ Hits          964      975      +11     
- Misses         52       74      +22

Impacted Files	Coverage Δ
astroscrappy/astroscrappy.pyx	`69.29% <47.50%> (-5.83%)`	⬇️
astroscrappy/utils/medutils.c	`100.00% <0.00%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update fd4acfe...0d45f31. Read the comment docs.

jehturner mentioned this pull request Sep 21, 2017

Add option bgsub, to pass a previously-subtracted background model cmccully/lacosmicx#1

Closed

crawfordsm mentioned this pull request Jun 15, 2018

CR rejection astropy/specreduce#21

Open

cmccully reviewed Dec 18, 2018

View reviewed changes

griffin-h reviewed Jan 11, 2019

View reviewed changes

jehturner force-pushed the spectroscopy branch from c43b73c to d628fc8 Compare March 6, 2019 15:05

jehturner added 4 commits August 4, 2020 18:17

Start off by adding the same bgsub parameter as in my fork of astrosc…

ba8f689

…rappy's predecessor, lacosmicx (producing identical results in a quick test).

Rename bgsub to bkg, remove pssl and temporarily subtract the supplie…

1bf914d

…d background internally instead of taking pre-subtracted input, as discussed with cmucully on the astropyspectroscopy slack channel in April.

Get the logic right for making a temporary m5 reference this time.

772482f

Add var parameter, to allow using a user-supplied variance array for …

0d45f31

…noise estimates, in place of the internal noise model. (May want some tweaking, eg. to protect bkg structure for noise estimates too??)

jehturner force-pushed the spectroscopy branch from d628fc8 to 0d45f31 Compare August 4, 2020 22:21

jehturner mentioned this pull request Sep 1, 2020

Feature/var and spectroscopy jehturner/astroscrappy#1

Closed

cmccully mentioned this pull request Nov 20, 2020

Feature/var and spectroscopy #53

Merged

cmccully closed this in #53 Feb 26, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

WIP: Support spectra & noise propagation #32

WIP: Support spectra & noise propagation #32

jehturner commented Sep 21, 2017

jehturner commented Sep 21, 2017

coveralls commented Sep 21, 2017

coveralls commented Sep 23, 2017 •

edited

Loading

jehturner commented Sep 23, 2017

jehturner commented Apr 16, 2018

jehturner commented Apr 16, 2018

cmccully left a comment

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

cmccully Dec 18, 2018

griffin-h Jan 11, 2019

jehturner Aug 7, 2020

jehturner commented Mar 6, 2019

jehturner commented Mar 6, 2019

jehturner commented Mar 8, 2019

jehturner commented Jul 4, 2020

cmccully commented Jul 31, 2020

jehturner commented Aug 4, 2020

jehturner commented Aug 4, 2020 •

edited

Loading

codecov bot commented Aug 4, 2020 •

edited

Loading

WIP: Support spectra & noise propagation #32

WIP: Support spectra & noise propagation #32

Conversation

jehturner commented Sep 21, 2017

jehturner commented Sep 21, 2017

coveralls commented Sep 21, 2017

coveralls commented Sep 23, 2017 • edited Loading

jehturner commented Sep 23, 2017

jehturner commented Apr 16, 2018

jehturner commented Apr 16, 2018

cmccully left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jehturner commented Mar 6, 2019

jehturner commented Mar 6, 2019

jehturner commented Mar 8, 2019

jehturner commented Jul 4, 2020

cmccully commented Jul 31, 2020

jehturner commented Aug 4, 2020

jehturner commented Aug 4, 2020 • edited Loading

codecov bot commented Aug 4, 2020 • edited Loading

Codecov Report

coveralls commented Sep 23, 2017 •

edited

Loading

jehturner commented Aug 4, 2020 •

edited

Loading

codecov bot commented Aug 4, 2020 •

edited

Loading