Repository: spacetelescope/tweakwcs (language: python)

Issue:
Add multivariate output fitting tool
Since many transforms between coordinate systems involve two input coordinates and two output coordinates, there is a need to adapt astropy modeling to be able to handle more than one output value to be fit. This provides such a tool that will accept an astropy model and astropy fitter instance along with the input values and output values to fit.
Add multivariate output fitting tool
Since many transforms between coordinate systems involve two input coordinates and two output coordinates, there is a need to adapt astropy modeling to be able to handle more than one output value to be fit. This provides such a tool that will accept an astropy model and astropy fitter instance along with the input values and output values to fit.


Discussion / hints:
# [Codecov](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
> Merging [#137](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) (5035ae7) into [master](https://codecov.io/gh/spacetelescope/tweakwcs/commit/a89b390139f0931a1d02d93209c0bb53e6f4f1e3?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) (a89b390) will **decrease** coverage by `0.53%`.
> The diff coverage is `67.85%`.

> :exclamation: Current head 5035ae7 differs from pull request most recent head fed6786. Consider uploading reports for the commit fed6786 to get more accurate results
[![Impacted file tree graph](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137/graphs/tree.svg?width=650&height=150&src=pr&token=YjZkswRHyv&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)

```diff
@@            Coverage Diff             @@
##           master     #137      +/-   ##
==========================================
- Coverage   85.93%   85.39%   -0.54%     
==========================================
  Files          10       11       +1     
  Lines        1827     1883      +56     
  Branches      270      277       +7     
==========================================
+ Hits         1570     1608      +38     
- Misses        252      269      +17     
- Partials        5        6       +1     
```


| [Impacted Files](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | |
|---|---|---|
| [tweakwcs/multivariate.py](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-dHdlYWt3Y3MvbXVsdGl2YXJpYXRlLnB5) | `67.85% <67.85%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope). Last update [a89b390...fed6786](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

@perrygreenfield @nden @stsci-hack 

Could you please describe the main use for this functionality that prompted this PR? That is, what kind of model (number of inputs and outputs) do you plan to use with this code? Is this like the test function in this PR: a model with 2 inputs and 2 outputs?
The main driver for this PR was the inability of the current astropy.modeling support for effectively replicating what 'tweakwcs' does; specifically, performing an 'rscale' fit between an input set of coordinates and a reference catalog of coordinates.  The astropy.modeling notebook used for JWSTUserTraining introduced the use of 2D Polynomials for fitting, but did not actually use them with a fitter because they did not work like the 1D case used earlier in the notebook.  The changes in this PR allow the type of fitting of 2D coordinates from the notebook to actually be fit using a fitter. The test implemented in this PR effectively demonstrates the basic reason for this PR.  

We need this functionality in astropy.modeling so that we can expand on the models we can fit beyond what 'tweakwcs' can do as 'tweakwcs' can only fit for linear global transformations between catalogs (coordinate frames). The ultimate goal is to be able to fit proper motions for each source in addition to the global transformations all at the same time which astropy.modeling can support once the problems with fitting addressed by this PR are implemented.  
> ... use of 2D Polynomials for fitting ...

I believe this particular fitting _is supported_: a 2D polynomial maps 2D domain to a 1D output. I believe the issue is when you want to fit **vector-valued** functions to, e.g., 2D data points.

For example, @perrygreenfield code in `gwcs.WCS.to_fits_sip()` already fits 2D Polynomials (SIP) to data points.
If you need to fit 2D-vector-models to 2D points, I have a different suggestion that works with existing `astropy` modeling infrastructure and is significantly simpler than proposed approach. This is an example based on your unit test. The main code is in the third block `In [82]`:

```python
In [80]: import numpy as np
    ...: from astropy.modeling.models import Shift, Rotation2D
    ...: from astropy.modeling.fitting import LevMarLSQFitter, _fitter_to_model_params

In [81]: inputs = [np.array([10., 10., 20., 20.]), np.array([10., 20., 20., 10.])]
    ...: outputs = [np.array([8.06101731, 0.98994949, 8.06101731, 15.13208512]),
    ...:            np.array([12.16223664, 19.23330445, 26.30437226, 19.23330445])]
    ...: r = Rotation2D()
    ...: r.fittable = True
    ...: model = (Shift() & Shift()) | r
    ...:

In [82]: class LMF2x2(LevMarLSQFitter):
    ...:     def objective_function(self, fps, *args):
    ...:         model, weights, inputs, meas = args
    ...:         _fitter_to_model_params(model, fps)
    ...:         return np.ravel(np.subtract(model(*inputs), meas))
    ...:

In [83]: f = LMF2x2()

In [84]: f(model, inputs, outputs)
Out[84]: <CompoundModel(offset_0=4.3, offset_1=-7.1, angle_2=45.)>
```

If you need weights - the code can be modified to work with weights.
I must admit this is a much easier way to do it. So let's dump this PR. Warren's eventual use case is more involved, but that is not relevant to this approach. I missed the fact that the objective function can be used to map multiple outputs into a single array.
# [Codecov](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
> Merging [#137](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) (5035ae7) into [master](https://codecov.io/gh/spacetelescope/tweakwcs/commit/a89b390139f0931a1d02d93209c0bb53e6f4f1e3?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) (a89b390) will **decrease** coverage by `0.53%`.
> The diff coverage is `67.85%`.

> :exclamation: Current head 5035ae7 differs from pull request most recent head fed6786. Consider uploading reports for the commit fed6786 to get more accurate results
[![Impacted file tree graph](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137/graphs/tree.svg?width=650&height=150&src=pr&token=YjZkswRHyv&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)

```diff
@@            Coverage Diff             @@
##           master     #137      +/-   ##
==========================================
- Coverage   85.93%   85.39%   -0.54%     
==========================================
  Files          10       11       +1     
  Lines        1827     1883      +56     
  Branches      270      277       +7     
==========================================
+ Hits         1570     1608      +38     
- Misses        252      269      +17     
- Partials        5        6       +1     
```


| [Impacted Files](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | |
|---|---|---|
| [tweakwcs/multivariate.py](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-dHdlYWt3Y3MvbXVsdGl2YXJpYXRlLnB5) | `67.85% <67.85%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope). Last update [a89b390...fed6786](https://codecov.io/gh/spacetelescope/tweakwcs/pull/137?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

@perrygreenfield @nden @stsci-hack 

Could you please describe the main use for this functionality that prompted this PR? That is, what kind of model (number of inputs and outputs) do you plan to use with this code? Is this like the test function in this PR: a model with 2 inputs and 2 outputs?
The main driver for this PR was the inability of the current astropy.modeling support for effectively replicating what 'tweakwcs' does; specifically, performing an 'rscale' fit between an input set of coordinates and a reference catalog of coordinates.  The astropy.modeling notebook used for JWSTUserTraining introduced the use of 2D Polynomials for fitting, but did not actually use them with a fitter because they did not work like the 1D case used earlier in the notebook.  The changes in this PR allow the type of fitting of 2D coordinates from the notebook to actually be fit using a fitter. The test implemented in this PR effectively demonstrates the basic reason for this PR.  

We need this functionality in astropy.modeling so that we can expand on the models we can fit beyond what 'tweakwcs' can do as 'tweakwcs' can only fit for linear global transformations between catalogs (coordinate frames). The ultimate goal is to be able to fit proper motions for each source in addition to the global transformations all at the same time which astropy.modeling can support once the problems with fitting addressed by this PR are implemented.  
> ... use of 2D Polynomials for fitting ...

I believe this particular fitting _is supported_: a 2D polynomial maps 2D domain to a 1D output. I believe the issue is when you want to fit **vector-valued** functions to, e.g., 2D data points.

For example, @perrygreenfield code in `gwcs.WCS.to_fits_sip()` already fits 2D Polynomials (SIP) to data points.
If you need to fit 2D-vector-models to 2D points, I have a different suggestion that works with existing `astropy` modeling infrastructure and is significantly simpler than proposed approach. This is an example based on your unit test. The main code is in the third block `In [82]`:

```python
In [80]: import numpy as np
    ...: from astropy.modeling.models import Shift, Rotation2D
    ...: from astropy.modeling.fitting import LevMarLSQFitter, _fitter_to_model_params

In [81]: inputs = [np.array([10., 10., 20., 20.]), np.array([10., 20., 20., 10.])]
    ...: outputs = [np.array([8.06101731, 0.98994949, 8.06101731, 15.13208512]),
    ...:            np.array([12.16223664, 19.23330445, 26.30437226, 19.23330445])]
    ...: r = Rotation2D()
    ...: r.fittable = True
    ...: model = (Shift() & Shift()) | r
    ...:

In [82]: class LMF2x2(LevMarLSQFitter):
    ...:     def objective_function(self, fps, *args):
    ...:         model, weights, inputs, meas = args
    ...:         _fitter_to_model_params(model, fps)
    ...:         return np.ravel(np.subtract(model(*inputs), meas))
    ...:

In [83]: f = LMF2x2()

In [84]: f(model, inputs, outputs)
Out[84]: <CompoundModel(offset_0=4.3, offset_1=-7.1, angle_2=45.)>
```

If you need weights - the code can be modified to work with weights.
I must admit this is a much easier way to do it. So let's dump this PR. Warren's eventual use case is more involved, but that is not relevant to this approach. I missed the fact that the objective function can be used to map multiple outputs into a single array.

Fix patch:
diff --git a/setup.py b/setup.py
index 0495502..74d5f21 100755
--- a/setup.py
+++ b/setup.py
@@ -77,6 +77,7 @@ def get_transforms_data():
     'pytest',
     'pytest-cov',
     'codecov',
+    'scipy'
 ]
 
 DOCS_REQUIRE = [
diff --git a/tweakwcs/linearfit.py b/tweakwcs/linearfit.py
index 77a04f0..7a898ad 100644
--- a/tweakwcs/linearfit.py
+++ b/tweakwcs/linearfit.py
@@ -11,14 +11,17 @@
 import logging
 import numbers
 import numpy as np
+from astropy.modeling.fitting import LevMarLSQFitter, _fitter_to_model_params
 
 from .linalg import inv
 from . import __version__  # noqa: F401
 
 __author__ = 'Mihai Cara, Warren Hack'
 
-__all__ = ['iter_linear_fit', 'build_fit_matrix', 'SUPPORTED_FITGEOM_MODES']
-
+__all__ = [
+    'iter_linear_fit', 'build_fit_matrix', 'SUPPORTED_FITGEOM_MODES',
+    '_LevMarLSQFitter2x2'
+]
 
 # Supported fitgeom modes and corresponding minobj
 SUPPORTED_FITGEOM_MODES = {
@@ -832,3 +835,14 @@ def build_fit_matrix(rot, scale=1):
                        [-sx * np.sin(rx), sy * np.cos(ry)]])
 
     return matrix
+
+
+class _LevMarLSQFitter2x2(LevMarLSQFitter):
+    """ Performs fits of 2D vector-models to 2D reference points. """
+    def objective_function(self, fps, *args):
+        model, weights, inputs, meas = args
+        _fitter_to_model_params(model, fps)
+        if weights is None:
+            return np.ravel(np.subtract(model(*inputs), meas))
+        else:
+            return np.ravel(weights * np.subtract(model(*inputs), meas))


Test patch:
diff --git a/tweakwcs/tests/test_linearfit.py b/tweakwcs/tests/test_linearfit.py
index bbdd234..64f90a2 100644
--- a/tweakwcs/tests/test_linearfit.py
+++ b/tweakwcs/tests/test_linearfit.py
@@ -8,6 +8,7 @@
 import math
 import pytest
 import numpy as np
+from astropy.modeling.models import Shift, Rotation2D
 from tweakwcs import linearfit, linalg
 
 
@@ -420,3 +421,16 @@ def test_iter_rscale_invalid_scale():
     # incorrect coordinate array dimensionality:
     with pytest.raises(ValueError):
         linearfit.fit_rscale(np.zeros((4, 2)), np.zeros((4, 2)), scale=0)
+
+
+def test_levmar2x2_multivariate():
+    inputs = [np.array([10., 10., 20., 20.]), np.array([10., 20., 20., 10.])]
+    outputs = [np.array([8.06101731, 0.98994949, 8.06101731, 15.13208512]),
+               np.array([12.16223664, 19.23330445, 26.30437226, 19.23330445])]
+    rot = Rotation2D()
+    rot.fittable = True
+    model = (Shift() & Shift()) | rot
+    fitter = linearfit._LevMarLSQFitter2x2()
+    finfo = fitter(model, inputs, outputs)
+    assert np.allclose(finfo.parameters, np.array([4.3, -7.1, 45.]),
+                       rtol=1e-5, atol=1e-5)