Fitting of Spectrum1D #210

brechmos-stsci · 2018-05-29T18:39:33Z

This PR enables the ability to fit a Spectrum1D object based on a list of astropy.modeling.models where each model in the list is set with a initial conditions. There is the ability to choose different astropy.modeling.fitters as well.

Not mergable. Just a WIP PR.

bsipocz · 2018-05-30T15:46:25Z

docs/fitting.rst

+
+   Original spectrum (blue) and fitted spectrum (orange).
+
+.. automodpai:: specuilts.fitting


argh.. darn fingers :). Thanks! Will fix.

oh, wait, not just 1 but 2 typos in one line! coffee, need more coffee

nmearl · 2018-05-31T04:45:29Z

I'm a bit confused by this approach. Why have the user provide a list of separate models? Why not have the user provide a single compound model? It seems that the fitmodel/_simple functions assume that the user only wants to add the list of models, but the astropy modeling machinery is super sophisticated, allowing arithmetic and binary operations to compose complicated compound models. It's not clear how this would be supported using this api abstraction. It also seems like the user has no control over the fitting initialization, since the fitmodel function only asks for an uninitialized class, and not a fitter instance (and the function doesn't have e.g. a fitter_kwargs dict).

Also, the docs seem to give the impression that the returned object is a list of models matching those that were put in, but that's a bit of an over simplification, yeah? The result is a model -- a compound model, which the user can query for a list of the composing sub models, or use directly as a simple model representation of the spectrum, or even use as an input into a more complex spectrum model.

From my understanding, this approach feels a little more complicated than simply asking the user to look at the astropy docs for how to compose a compound model, and feed that, along with some spectrum data, into a fitter. I think I might just be missing the bigger picture?

brechmos-stsci · 2018-05-31T12:24:44Z

@nmearl This is a work in progress. @eteq and I had a long discussion about the steps to fitting spectra and how to make it as modular as possible. This is the last of the fitting steps in that framework he and I discussed. He is firming up the framework and interface from our discussion. Once that is completed, I will come back to this work.

cc @eteq

nmearl · 2018-05-31T13:50:20Z

Great! Sorry to jump the gun.

brechmos-stsci · 2018-07-12T13:15:54Z

@eteq In the https://github.com/spacetelescope/dat_pyinthesky/blob/master/pyinthesky_specutils_fitting.ipynb you have the window parameter for the fitting. This essentially tells the fitter to do the fit only within that window/region for a model (if multiple models then would have multiple windows, one for each).

Given the analysis work and discussion of region, I wonder if it is better to rename window to region. They are effectively the same thing ("area on which to work"). I think generally they are similar enough that they should be named the same thing for consistency.

bsipocz · 2018-08-13T16:16:00Z

@brechmos-stsci - it seems that this will need a rebase as many unrelated changes has been merged in from master.

brechmos-stsci · 2018-08-14T13:03:43Z

Going down a rabbit hole and need some guidance. For fitting one could potentially pass in a compound astropy model, but astropy fitting currently can not fit a compound model that has sub-models that has parameters with astropy units.

So, my plan was to convert the compound model and spectral data to a common set of units, strip the units, do the fitting, add the units back onto the fitted model and then convert fitted model units back to the original units.

So, going down this rabbit hole:

Attempt 1) I originally was just going take each sub-model of the compound model, create a new sub-model of the same type, then convert and set its parameters. The problem is for models such as Chebyshev1D one needs to know the 'order' to create a new instance. So, not easy to do model creation in simple, general way given some models have input that define the number of parameters.

Attempt 2) So, rather than create a new sub-model I was going to do a deepcopy() of each sub-model and convert, strip and set the parameter values. But in doing this method on Gaussian1D, if the input sub-model amplitude has a Quantity parameter then I can't set a non-Quantity value (it errors out).

So, basically, I have gone down this rabbit hole and there has to be a better way to do this. It all comes down to the fact that astropy fitting does not like compound models that have units.

@eteq @stscicrawford @nmearl @bsipocz Any ideas on something a little cleaner I could do?

[Update: I figured out a way to do it with stripping units, but I still wonder if there is a better way to do it]

nmearl · 2018-08-14T17:57:49Z

This has been an on going issue, I have a couple of open issues with this. I tackled the same thing when developing spectacle. To avoid this issue, I had to ensure that every model used provided a _parameter_units_for_data_units method. This is because astropy fitters natively try to remove units before fitting, and then will attempt to add them back if it can.

Finally, for some reason with a long winded explanation, I had to create defined classes for the compound models. This may or may not be necessary, but it seemed to alleviate some issues with astropy dealing directly with fitting a compound model as opposed to just what it thinks is a custom model.

I'll certainly comb through your PR and see if I can't offer some extra help if I can.

bsipocz · 2018-08-14T18:10:56Z

compound models are also due to be refactored, but I don't exactly know the timescale and extent of it though, you may want to check it with Nadia and Perry.

brechmos-stsci · 2018-08-15T20:02:25Z

After numerous discussion with Nadia, Tom, etc etc, I now have a working set of code (ie it fits compound models with units). @bsipocz, understood, but not sure how long that will be to get into astropy.

This PR is not quite ready, but it passes tests on my machine and I plotted the data and fitted model/compound model and it looks decent:

I will need to clean up this PR a bit.

So, as a summary, astropy fitting will not fit a CompoundModel that has units. So, I stripped the units off the compound model and created a new model. After discussions with several people, it seems the only way to get a representation of the submodels with operators is to use model._tree.traverse_postorder() SO, this uses a private variable (tsk tsk). If someone has a better way I could do this, please tell me. Then, once I have the traversal I can re-create the compound model without units (after conversion to spectral units), fit, and then add the units back on (same type of procedure). When I add the units back on, the models in the tree of the compound model does not have updated parameters so I have to grab the fitted values from the compound models parameters list. There are many subtleties here.

Tomorrow, I might clean up this entry but I wanted to get all my thoughts down before I go for the day.

brechmos-stsci · 2018-08-17T15:15:37Z

Ok. after long last, many discussions with people about astropy modeling, compound models and the idiosyncrasies I finally have a version of spectra fitting that works (for at least the test cases :) ).

@stscicrawford @eteq @bsipocz and anyone else, please take a look through it.

Here are the caveats:

It uses a private variable ._tree in CompoundModel. After talking with Perry, Nadia and some others this appears to be the only way to get the CompoundModel operators.
- I looked at removing the units of each parameter of each sub-model in the CompoundModel but if the parameter of a Model is a Quantity it can only be replaced by a Quantity -- so no-go on that one.
- I also tried creating a new unit-less model and then just placing into the CompoundModel, but CompoundModel does not allow indexed assignment of a new model (e.g., compound_model[0] = gaussian_unitless_model does not work)
In the future, when astropy fitting has the ability to fit data to CompoundModels, where both have units, then much of the code in fitmodels.py can be removed.
Also in the future, when CompoundModels is re-written this will have to be re-addressed.

eteq

Lots of inline comments, most pretty straightforward typos and the like. A few that didn't fit inline:

Should have a section on continuum fitting in the docs. For now this could just be a quick stub that shows a call to the generic fitter (and a few words on what "generic" means here), to be later expanded if/when there are cleverer continuum-fitting functions.
I'm confused why the exciser is in utils rather than a method on Spectrum1D or SpectralRegion. Or is it that #289 is effectively replacing this anyway and it's just a parallel-PR's issue?
I really like the examples you showed in your comment #210 (comment) . Am I right that's from the tests? I think it would be neat to turn that into a "gallery" in the docs - you could use the .. plot:: directive to implement all of them by just importing the test functions and plotting the result?

docs/fitting.rst

eteq · 2018-08-31T17:54:45Z

docs/fitting.rst

+    >>> np.random.seed(42)
+
+
+For the purpose of the example, build a ``spectrum1d`` variable that will be used in the fitting:


"variable" -> "object"

eteq · 2018-08-31T17:57:03Z

docs/fitting.rst

+    >>> spectrum1d = Spectrum1D(spectral_axis=wave_um, flux=flux_e1)
+
+
+.. figure:: img/fitting_original.jpg


This is the result of the above, right? If you use the plot sphinx directive, it will automatically do this instead of you having to manage the image files.

See e.g. http://docs.astropy.org/en/stable/visualization/normalization.html#matplotlib-normalization
and the corresponding rst at
https://github.com/astropy/astropy/blob/master/docs/visualization/normalization.rst

I'll take a look at this.

eteq · 2018-08-31T17:58:40Z

docs/fitting.rst

+
+   Spectrum to be fit.
+
+Now that there is a ``spectrum1d`` to fit, the real fitting setup must happen.  First create


``spectrum1d``

->

`~specutils.Spectrum1D`

(to make it a live link)

specutils/fitting/fitmodels.py

eteq · 2018-08-31T21:24:29Z

specutils/fitting/fitmodels.py

+
+    # In this case the window defines the area around the center of each model
+    if window is not None and isinstance(window, (float, int)):
+        center = model.mean


Not all models have the same definition of "center". My instinct is to right a heuristic function to determine the center - it should work for Gaussian, Lorentzian, and Voigt, and also try some obvious ones like "mean" or "center", but other than that just give up and tell the user "sorry, you've used too weird of a model for me to know what you mean".

That said, a valid response is "won't do in this PR, make an issue to track it for later"!

eteq · 2018-08-31T21:33:28Z

specutils/fitting/fitmodels.py

+        # Make a new instance of the class.
+        #
+
+        if isinstance(sub_model, models.PolynomialModel):


Alternative solution to this that's a bit more general:

if any([p.default == p.empty for p in inspect.signature(sub_model.__class__).parameters.values() if p.kind==p.POSITIONAL_OR_KEYWORD])

You might want to unpack that into multiple lines (😉 ), but it looks to me like it does the job: if a model class has a keyword without a default, it will be treated like Polynomial. I think that's basically the key thing here, right? That there's no default?

You'd probably then have to programatically pick out the parameter name using the signature to replace the sub_model.degree in the next line, but that should work.

(Also falls under "you can assign this to @eteq since it's his crazy idea, though 😉)

Ok. I would like to suggest this goes into a new PR once this PR is merged.

specutils/fitting/fitmodels.py

brechmos-stsci · 2018-09-04T13:47:32Z

I'm confused why the exciser is in utils rather than a method on Spectrum1D or SpectralRegion. Or is it that #289 is effectively replacing this anyway and it's just a parallel-PR's issue?

It is the #289 / parallel PR thing... I'll modify to use SpectralRegions

brechmos-stsci · 2018-09-04T18:19:15Z

I really like the examples you showed in your comment ...

They are from the tests. I added them into the docs using ..plot.

drdavella · 2018-09-04T18:52:38Z

Current test failures are unrelated and need to be addressed upstream in https://github.com/astropy/ci-helpers. A workaround in the meantime would be setting CONDA_VERSION in .travis.yml to something sufficiently recent (like 4.5?).

eteq · 2018-09-06T21:07:20Z

I made some doc updates in 3273b62 as promised, but other than that (and a series of issues I'll create for enhancements), looks good, and will merge as soon as tests pass.

Bad data popup

brechmos-stsci changed the title ~~WIP: Fitting of Spectrum1D~~ Fitting of Spectrum1D May 30, 2018

bsipocz reviewed May 30, 2018

View reviewed changes

brechmos-stsci changed the title ~~Fitting of Spectrum1D~~ WIP: Fitting of Spectrum1D May 31, 2018

eteq mentioned this pull request Jul 17, 2018

Line Fitting #247

Closed

brechmos added 17 commits August 17, 2018 10:32

initial fitting structure

95f40d6

added more documentation

01e8977

fixed for rename

83a91bc

added first pass at fitting document

a215257

added documentation, small dir refactor

d975f2c

fixed *two* typos in the automodapi

98659a9

fixed change in directory structure

084dba7

added in region cutouts (excise) of spectra.

8cc476b

updated api and tests

b20ce6a

changed name to fit_lines to match erick's notebook

5a29bd5

udpated the continuum fitting and tests

6e237d5

set the tolerance

dcf981f

first pass at cleanup of fitting (model and compoundmodels)

f14f0f9

continuum fixes

fa0c5f6

fixed for Chebyshev1D model

52ef187

remove units, fit, add units back

a023b4b

works for just adding models

8bf3ecb

brechmos-stsci force-pushed the fitting branch from a212ac9 to c9b77eb Compare August 17, 2018 14:39

added fitting into index doctree

3b4597f

brechmos-stsci changed the title ~~WIP: Fitting of Spectrum1D~~ Fitting of Spectrum1D Aug 17, 2018

eteq mentioned this pull request Aug 24, 2018

Add a check for whether the spectrum is continuum subtracted for centroid #290

Closed

eteq requested changes Aug 31, 2018

View reviewed changes

brechmos added 2 commits September 4, 2018 09:15

added text re initial guess code will be added later.

ffa15f4

updates and mods based on PR comments

6c37563

brechmos added 3 commits September 4, 2018 12:53

updated to changed method name

b9fe3e4

added dynamically generated plots

a8d55a2

removed unnecesary static images

906bef9

added continuum fitting

cdc3972

nmearl added this to the v0.6.0 milestone Sep 5, 2018

eteq mentioned this pull request Sep 5, 2018

[WIP] Continuum normalization #291

Open

4 tasks

brechmos-stsci mentioned this pull request Sep 6, 2018

Centroid calculation #270

Merged

more contextual words about fitting in the docs

3273b62

eteq approved these changes Sep 6, 2018

View reviewed changes

eteq merged commit 9c96d21 into astropy:master Sep 6, 2018

brechmos-stsci deleted the fitting branch December 5, 2018 19:07

SaOgaz pushed a commit that referenced this pull request Mar 25, 2019

Merge pull request #210 from brechmos-stsci/bad-data-popup

c60a442

Bad data popup

nmearl mentioned this pull request Apr 11, 2019

Gwcs refactor #462

Merged

2 tasks


		Original spectrum (blue) and fitted spectrum (orange).

		.. automodpai:: specuilts.fitting

		>>> np.random.seed(42)


		For the purpose of the example, build a ``spectrum1d`` variable that will be used in the fitting:

		>>> spectrum1d = Spectrum1D(spectral_axis=wave_um, flux=flux_e1)


		.. figure:: img/fitting_original.jpg


		Spectrum to be fit.

		Now that there is a ``spectrum1d`` to fit, the real fitting setup must happen. First create

Uh oh!

Fitting of Spectrum1D #210

Fitting of Spectrum1D #210

Uh oh!

Conversation

brechmos-stsci commented May 29, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brechmos-stsci May 30, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nmearl commented May 31, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

brechmos-stsci commented May 31, 2018

Uh oh!

nmearl commented May 31, 2018

Uh oh!

brechmos-stsci commented Jul 12, 2018

Uh oh!

bsipocz commented Aug 13, 2018

Uh oh!

brechmos-stsci commented Aug 14, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nmearl commented Aug 14, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

bsipocz commented Aug 14, 2018

Uh oh!

brechmos-stsci commented Aug 15, 2018

Uh oh!

brechmos-stsci commented Aug 17, 2018

Uh oh!

eteq left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

brechmos-stsci commented Sep 4, 2018

Uh oh!

brechmos-stsci commented Sep 4, 2018

Uh oh!

drdavella commented Sep 4, 2018

Uh oh!

eteq commented Sep 6, 2018

Uh oh!

Reviewers

Assignees

Labels

Projects

brechmos-stsci May 30, 2018 •

edited

Loading

nmearl commented May 31, 2018 •

edited

Loading

brechmos-stsci commented Aug 14, 2018 •

edited

Loading

nmearl commented Aug 14, 2018 •

edited

Loading