initial bones of resample code, just have flux conserving and interpolation right now.

[SaOgaz]

See related issues: #359, #254, #255, #266

[SaOgaz]

Questions:

I discussed this a little bit with @eteq, the initial suggestion for weighting (and or uncertainity) was to have a weight keyword that could have one of three values, None (meaning no weighting), provided array (use provided array as weights) or 'uncer' (use the uncertainty defined in the input Spectrum1D as the weights). If the user choses 'uncer' but no uncertainty is in the Spectrum1D, we throw an error. Right now this is unimplemented.

How do we want to handle propagating the uncertainties.... the paper I pulled the algorithm from has a prescription for calculating the new associated uncertainties (it's just using normal uncertainty propogation). I haven't looked into how easy it would be to adapt the current code to use any built in uncertainty handling, buy my initial guess is that it won't be too difficult, since the flux doesn't get multiplied through until the very end. The thing I'm not sure about is, if someone is using the uncertainties for weights, how does that effect the output uncertainty?

Another uncertainty related question. if they are not defined in the input spectrum, do we leave them undefined? I'm leaning towards yes.

[hcferguson]

Anyone familiar with the uncertainties package?

[eteq]

@SaOgaz

I would say for this PR we just follow the propagation the paper suggests. The key point is that we need to be sure it's the right kind of uncertainty, though. That is, I think the paper is built around gaussian statistics, so you do that propagation if the uncertainty in the spectrum is either SD, variance, or ivar uncertainties (which is probably like 99% of the use cases anyway). One thing about this, though: if the keyword argument is weights, the "use the unceratinty" option should describing the weighting. I.e., instead of "uncer" I'd say it should be ivar_weighted, which means "convert the uncertainty to an inverse variance and use that as the weights" (which is also how you'd implement it).

As for the question of what to do if the weights aren't uncertainties - aren't the weights just the c_ij from https://arxiv.org/abs/1705.05165 ? In that case you should just be able to use the expression exactly as given in section 5 (the covariance matrix), regardless of whether the weights come from the uncertainties or somewhere else.

Related thought: since we are following that algorithm, perhaps the uncertainty should record the whole covariance matrix not just the diagonals (that is, the ). That would require making a new uncertainty class, but it would be easy to make one that can be converted to VarianceUncertainty or StdDevUncertainty by just taking the diagonal of the matrix? This may well be better done as a follow-on PR, though.

@hcferguson

Anyone familiar with the uncertainties package?

Yes, several years back I had a GSoC student work on that to see how well it can be used with some of the Astropy machinery. At the time the answer was "it's very hard because they've done some numpy stuff that conflicts deeply with a lot of other machinery". It may be useful as an eventual thing to complement the (fairly new) existing machinery we built in astropy to do uncertainty propogation (http://docs.astropy.org/en/stable/uncertainty/index.html). But it's hard to just use as it stands.

[eteq]

One other very minor comment to just make sure I don't forget: the resample functions should be pulled out into specutils.manipulation via from .resample import * in the __init__.py in manipulation. We agreed we generally want only 2-level imports at the most for the user API.

[eteq]

@SaOgaz - as we discussed out-of-band, I've added a test here to check for flux conservation. Basically the test I added is checking that flux is conserved between an input and output spectrum in the way I understand us to be aiming for. When I ran locally this didn't work... But we shall see what travis says.

[SaOgaz]

@eteq, huh, I did that a bit locally (not sure why I didn't think to add a test for it), and it came out correctly, but maybe my calculations were off?

[SaOgaz]

@eteq ahh, I think I might know why those tests are failing. For your input data, the bin edges aren't aligned, just input center wavelengths. i.e. the left edge of the non linespace wavelength axis is 3,500, the linespace left edge is 3666.66. The way I implemented this is that the flux conservation will only truly apply when the bin edges line up. (although i'm trying to edit your test example to follow this case and running into problems, so cannot yet confirm that results come out right aside from the testing i did earlier).

[SaOgaz]

side note, meant to add to init, just hadn't gotten around to it yet.

[SaOgaz]

As for the question of what to do if the weights aren't uncertainties - aren't the weights just the c_ij from https://arxiv.org/abs/1705.05165 ? In that case you should just be able to use the expression exactly as given in section 5 (the covariance matrix), regardless of whether the weights come from the uncertainties or somewhere else.

edit: sorry I'm getting a bit confused, I thought we were just using quadrature in the case of no, or separate weights (where separate weights would be including in the quadrature calculation) as described in the "calculation of associated uncertainties" section...

[SaOgaz]

I would say for this PR we just follow the propagation the paper suggests. The key point is that we need to be sure it's the right kind of uncertainty, though. That is, I think the paper is built around gaussian statistics, so you do that propagation if the uncertainty in the spectrum is either SD, variance, or ivar uncertainties (which is probably like 99% of the use cases anyway). One thing about this, though: if the keyword argument is weights, the "use the unceratinty" option should describing the weighting. I.e., instead of "uncer" I'd say it should be ivar_weighted, which means "convert the uncertainty to an inverse variance and use that as the weights" (which is also how you'd implement it).

Is this in the situation where someone isn't using one of the standard NDUncertainty subclasses? I would think that just checking the type and assuming we want to return the same type would be the most straightforward approach, in which case using a bland "uncer" to mean, use whatever uncertainty is there as a weight would be the best way to go.

[eteq]

@SaOgaz

edit: sorry I'm getting a bit confused, I thought we were just using quadrature in the case of no, or separate weights (where separate weights would be including in the quadrature calculation) as described in the "calculation of associated uncertainties" section...

Ah right, I see your point here, and I missed some of the subtlety. I think for the purposes of this PR, the "weights" are basically 1/(sigma^2) in the notation of the "Calculation of associated uncertainties" section. That doesn't make a ton of sense really in the single-spectrum case, but it's relevant for the "list of spectra" case. And in fact this is making me realize maybe we don't want to call this weights but rather pixel_weights. We might later add a spectrum_weights or similar, but given the fact that I think usually the user will want to use the uncertainty rather than the weights array, this makes sense as a way to keep those concepts separate.

Is this in the situation where someone isn't using one of the standard NDUncertainty subclasses? I would think that just checking the type and assuming we want to return the same type would be the most straightforward approach, in which case using a bland "uncer" to mean, use whatever uncertainty is there as a weight would be the best way to go.

Yes, if there's an NDUncertainty class definitely you use that information. (maybe write a simple to_ivar function that converts them all? really we might want to pass that upstream to astropy too but for now it's probably fine to have it be a private utility in specutils) My comment here was about the possibility of the user providing a "weights" array, in which case you'd need to say what kind of weight it was. But given the above realization I had, I think that's no longer relevant. Instead we just say "if weights are arrays they are interpreted as inverse variance" and then it's the user's business to know how to convert their weights into that. So I think you can ignore this concern now.

[eteq]

Just to note here - after some out-of-band discussion with @SaOgaz (included in the commits above) we decided the best way to push forward is to not have the weights at all for this PR. Rather instead this implements the flux-conserving algorithm (which doesn't really have meaningful per-pixel weights as described above if you want it to be flux-conserving...), but does propogate the uncertainties self consistently.

The whole "should there be a separate pixel_weight argument" is still relevant, especially in the context of other resampling algorithms (e.g. interpolation), but we decided that should be dealt with after this is implemented so that we can get something in to build off of.

[eteq]

(also the test failure appears to be unrealted to this PR, as it's also in master)

[SaOgaz]

I've added two interpolation versions as well, linear and spline.

[eteq]

A few minor comments, and one larger one (about splitting the InterpolatedResample) although all are up for debate if others disagree.

Also, it would be good to add at least a brief example in the docs about this. The most logical place is probably the "Manipulating Spectra" section? (manipulation.rst)

[eteq]

Can you reference either a PR # or a commit on github to cite this just so that there's history?

[SaOgaz]

Sure, there's no PR though. I was referencing his branch off his fork: nmearl@0ff6ef1. What's the best way to reference that?

[eteq]

I think it's find just to give the URL directly to that SHA (e.g. nmearl@0ff6ef1 ) - this is a private method so no need to be too picky about it

[SaOgaz]

Oh shoot, @eteq I accidentally over-wrote the change of yours I merged in, what did you change "This resample algorithim conserves overall flux (as opposed to flux density)" to?

[stscicrawford]

Tagging @philhodge, @jdavies-st, and @stscirij to take a look at how useful this functionality would be for the calibration software or if there is something that could be added in the future

[eteq]

I few fairly straightforward suggestions in some of the stuff I haven't yet reviewed.

One "find-and-replace" suggestion I'd like to make additionally, though. Right now the class names are all XXXResample. I think these would be better as XXXResampler - that is, class names should be nouns, because it's an "object", whereas "resample" is a verb. So the class should be "a thing that does resampling" - i.e. a "resampler". (Same comment for the variable names in the examples, but that's probably still just a global find-and-replace of "resample" -> "resampler"?)

[eteq]

I think it's find just to give the URL directly to that SHA (e.g. nmearl@0ff6ef1 ) - this is a private method so no need to be too picky about it

[eteq]

Alright I think we're all good here now. Thanks for your work on this @SaOgaz !