Numerical instability in Mapping.rate

[TallJimbo]

I'm seeing an unexpected near-zero value in Mapping.rate for a transform that otherwise looks perfectly fine (with the correct rate many orders of magnitude larger). In particular, this is a 2-d composite mapping, and if I compute the first derivative matrices of the two component mappings via Mapping.rate, those matrices and their product looks fine, but the matrix derived from the composite mapping's rate is singular.

I've included below a how-to-reproduce Python script that uses the Python bindings in Rubin Observatory's astshim package, but I think it should be easy translatable into any other kind of script that calls AST under the hood. If you do need to install it, I'd recommend installing the full Rubin software stack.

#!/usr/bin/env python

import astshim
import numpy as np

patch_to_sky = astshim.Mapping.fromString("""
Begin CmpMap 	# Compound Mapping
    Nin = 2 	# Number of input coordinates
    IsSimp = 1 	# Mapping has been simplified
 IsA Mapping 	# Mapping between coordinate systems
    MapA = 	# First component Mapping
       Begin WinMap 	# Map one window on to another
          Nin = 2 	# Number of input coordinates
       IsA Mapping 	# Mapping between coordinate systems
          Sft1 = 0.014543440805923862 	# Shift for axis 1
          Scl1 = -9.696273622190721e-07 	# Scale factor for axis 1
          Sft2 = -0.014543440805923862 	# Shift for axis 2
          Scl2 = 9.696273622190721e-07 	# Scale factor for axis 2
       End WinMap
    MapB = 	# Second component Mapping
       Begin CmpMap 	# Compound Mapping
          Nin = 2 	# Number of input coordinates
       IsA Mapping 	# Mapping between coordinate systems
          InvA = 1 	# First Mapping used in inverse direction
          MapA = 	# First component Mapping
             Begin WcsMap 	# FITS-WCS sky projection
                Nin = 2 	# Number of input coordinates
                Invert = 1 	# Mapping inverted
             IsA Mapping 	# Mapping between coordinate systems
                Type = "TAN" 	# Gnomonic projection
             End WcsMap
          MapB = 	# Second component Mapping
             Begin CmpMap 	# Compound Mapping
                Nin = 2 	# Number of input coordinates
             IsA Mapping 	# Mapping between coordinate systems
                InvA = 1 	# First Mapping used in inverse direction
                MapA = 	# First component Mapping
                   Begin SphMap 	# Cartesian to Spherical mapping
                      Nin = 3 	# Number of input coordinates
                      Nout = 2 	# Number of output coordinates
                      Invert = 1 	# Mapping inverted
                   IsA Mapping 	# Mapping between coordinate systems
                      UntRd = 1 	# All input vectors have unit length
                      PlrLg = 0 	# Polar longitude (rad.s)
                   End SphMap
                MapB = 	# Second component Mapping
                   Begin CmpMap 	# Compound Mapping
                      Nin = 3 	# Number of input coordinates
                      Nout = 2 	# Number of output coordinates
                   IsA Mapping 	# Mapping between coordinate systems
                      MapA = 	# First component Mapping
                         Begin MatrixMap 	# Matrix transformation
                            Nin = 3 	# Number of input coordinates
                         IsA Mapping 	# Mapping between coordinate systems
                            M0 = -0.27751547555781247 	# Forward matrix value
                            M1 = -0.79955425190606422 	# Forward matrix value
                            M2 = 0.53263323129978957 	# Forward matrix value
                            M3 = -0.36944878902416273 	# Forward matrix value
                            M4 = 0.60059387131316444 	# Forward matrix value
                            M5 = 0.70908010409942102 	# Forward matrix value
                            M6 = -0.88684426655106752 	# Forward matrix value
                            M7 = 0 	# Forward matrix value
                            M8 = -0.46206844394039626 	# Forward matrix value
                            IM0 = -0.27751547555781247 	# Inverse matrix value
                            IM1 = -0.36944878902416273 	# Inverse matrix value
                            IM2 = -0.88684426655106752 	# Inverse matrix value
                            IM3 = -0.79955425190606422 	# Inverse matrix value
                            IM4 = 0.60059387131316444 	# Inverse matrix value
                            IM5 = 0 	# Inverse matrix value
                            IM6 = 0.53263323129978957 	# Inverse matrix value
                            IM7 = 0.70908010409942102 	# Inverse matrix value
                            IM8 = -0.46206844394039626 	# Inverse matrix value
                            Form = "Full" 	# Matrix storage form
                         End MatrixMap
                      MapB = 	# Second component Mapping
                         Begin SphMap 	# Cartesian to Spherical mapping
                            Nin = 3 	# Number of input coordinates
                            Nout = 2 	# Number of output coordinates
                         IsA Mapping 	# Mapping between coordinate systems
                            UntRd = 1 	# All input vectors have unit length
                            PlrLg = 0.92655267202648273 	# Polar longitude (rad.s)
                         End SphMap
                   End CmpMap
             End CmpMap
       End CmpMap
 End CmpMap
 """)

sky_to_ccd = astshim.Mapping.fromString("""
 Begin CmpMap 	# Compound Mapping
    Nin = 2 	# Number of input coordinates
    IsSimp = 1 	# Mapping has been simplified
 IsA Mapping 	# Mapping between coordinate systems
    InvA = 1 	# First Mapping used in inverse direction
    MapA = 	# First component Mapping
       Begin SphMap 	# Cartesian to Spherical mapping
          Nin = 3 	# Number of input coordinates
          Nout = 2 	# Number of output coordinates
       IsA Mapping 	# Mapping between coordinate systems
          UntRd = 1 	# All input vectors have unit length
          PlrLg = 0.92910062914125369 	# Polar longitude (rad.s)
       End SphMap
    MapB = 	# Second component Mapping
       Begin CmpMap 	# Compound Mapping
          Nin = 3 	# Number of input coordinates
          Nout = 2 	# Number of output coordinates
       IsA Mapping 	# Mapping between coordinate systems
          InvA = 1 	# First Mapping used in inverse direction
          MapA = 	# First component Mapping
             Begin MatrixMap 	# Matrix transformation
                Nin = 3 	# Number of input coordinates
             IsA Mapping 	# Mapping between coordinate systems
                M0 = -0.28337891056398906 	# Forward matrix value
                M1 = -0.80108194229191043 	# Forward matrix value
                M2 = 0.52722302186208359 	# Forward matrix value
                M3 = -0.37926313225344088 	# Forward matrix value
                M4 = 0.59855469402037131 	# Forward matrix value
                M5 = 0.70561445193497041 	# Forward matrix value
                M6 = -0.88082681019646292 	# Forward matrix value
                M7 = 0 	# Forward matrix value
                M8 = -0.47343862372975476 	# Forward matrix value
                IM0 = -0.28337891056398906 	# Inverse matrix value
                IM1 = -0.37926313225344088 	# Inverse matrix value
                IM2 = -0.88082681019646292 	# Inverse matrix value
                IM3 = -0.80108194229191043 	# Inverse matrix value
                IM4 = 0.59855469402037131 	# Inverse matrix value
                IM5 = 0 	# Inverse matrix value
                IM6 = 0.52722302186208359 	# Inverse matrix value
                IM7 = 0.70561445193497041 	# Inverse matrix value
                IM8 = -0.47343862372975476 	# Inverse matrix value
                Form = "Full" 	# Matrix storage form
             End MatrixMap
          MapB = 	# Second component Mapping
             Begin CmpMap 	# Compound Mapping
                Nin = 3 	# Number of input coordinates
                Nout = 2 	# Number of output coordinates
             IsA Mapping 	# Mapping between coordinate systems
                MapA = 	# First component Mapping
                   Begin SphMap 	# Cartesian to Spherical mapping
                      Nin = 3 	# Number of input coordinates
                      Nout = 2 	# Number of output coordinates
                      Invert = 1 	# Mapping inverted
                   IsA Mapping 	# Mapping between coordinate systems
                      UntRd = 1 	# All input vectors have unit length
                      PlrLg = 0 	# Polar longitude (rad.s)
                   End SphMap
                MapB = 	# Second component Mapping
                   Begin CmpMap 	# Compound Mapping
                      Nin = 2 	# Number of input coordinates
                   IsA Mapping 	# Mapping between coordinate systems
                      MapA = 	# First component Mapping
                         Begin WcsMap 	# FITS-WCS sky projection
                            Nin = 2 	# Number of input coordinates
                            Invert = 1 	# Mapping inverted
                         IsA Mapping 	# Mapping between coordinate systems
                            Type = "TAN" 	# Gnomonic projection
                         End WcsMap
                      MapB = 	# Second component Mapping
                         Begin CmpMap 	# Compound Mapping
                            Nin = 2 	# Number of input coordinates
                         IsA Mapping 	# Mapping between coordinate systems
                            InvA = 1 	# First Mapping used in inverse direction
                            MapA = 	# First component Mapping
                               Begin WinMap 	# Map one window on to another
                                  Nin = 2 	# Number of input coordinates
                               IsA Mapping 	# Mapping between coordinate systems
                                  Sft1 = -0.00064954183224266336 	# Shift for axis 1
                                  Scl1 = 0.017453292519943295 	# Scale factor for axis 1
                                  Sft2 = 0.00016792057273134152 	# Shift for axis 2
                                  Scl2 = 0.017453292519943295 	# Scale factor for axis 2
                               End WinMap
                            MapB = 	# Second component Mapping
                               Begin CmpMap 	# Compound Mapping
                                  Nin = 2 	# Number of input coordinates
                               IsA Mapping 	# Mapping between coordinate systems
                                  InvA = 1 	# First Mapping used in inverse direction
                                  MapA = 	# First component Mapping
                                     Begin MatrixMap 	# Matrix transformation
                                        Nin = 2 	# Number of input coordinates
                                     IsA Mapping 	# Mapping between coordinate systems
                                        M0 = -25.534438635766644 	# Forward matrix value
                                        M1 = 51.302461187649776 	# Forward matrix value
                                        M2 = 51.298795950170486 	# Forward matrix value
                                        M3 = 25.532559965718377 	# Forward matrix value
                                        IM0 = -0.0077755125420108477 	# Inverse matrix value
                                        IM1 = 0.015623303379535295 	# Inverse matrix value
                                        IM2 = 0.015622187192986517 	# Inverse matrix value
                                        IM3 = 0.0077760846594382505 	# Inverse matrix value
                                        Form = "Full" 	# Matrix storage form
                                     End MatrixMap
                                  MapB = 	# Second component Mapping
                                     Begin CmpMap 	# Compound Mapping
                                        Nin = 2 	# Number of input coordinates
                                     IsA Mapping 	# Mapping between coordinate systems
                                        MapA = 	# First component Mapping
                                           Begin UnitNormMap 	# Compute unit vector and norm relative to a centre
                                              Nin = 2 	# Number of input coordinates
                                              Nout = 3 	# Number of output coordinates
                                              Invert = 1 	# Mapping inverted
                                           IsA Mapping 	# Mapping between coordinate systems
                                           End UnitNormMap
                                        MapB = 	# Second component Mapping
                                           Begin CmpMap 	# Compound Mapping
                                              Nin = 3 	# Number of input coordinates
                                              Nout = 2 	# Number of output coordinates
                                           IsA Mapping 	# Mapping between coordinate systems
                                              MapA = 	# First component Mapping
                                                 Begin CmpMap 	# Compound Mapping
                                                    Nin = 3 	# Number of input coordinates
                                                 IsA Mapping 	# Mapping between coordinate systems
                                                    Series = 0 	# Component Mappings applied in parallel
                                                    InvB = 1 	# Second Mapping used in inverse direction
                                                    MapA = 	# First component Mapping
                                                       Begin UnitMap 	# Unit (null) Mapping
                                                          Nin = 2 	# Number of input coordinates
                                                       IsA Mapping 	# Mapping between coordinate systems
                                                       End UnitMap
                                                    MapB = 	# Second component Mapping
                                                       Begin PolyMap 	# Polynomial transformation
                                                          Nin = 1 	# Number of input coordinates
                                                       IsA Mapping 	# Mapping between coordinate systems
                                                          MPF1 = 5 	# Max. power of input 1 in any forward polynomial
                                                          NCF1 = 5 	# No. of coeff.s for forward polynomial 1
                                                          CF1 = 9.7050637265321229e-05 	# Coeff 1 of forward polynomial 1
                                                          CF2 = 0 	# Coeff 2 of forward polynomial 1
                                                          CF3 = 9.0056225011787966e-12 	# Coeff 3 of forward polynomial 1
                                                          CF4 = 0 	# Coeff 4 of forward polynomial 1
                                                          CF5 = 2.346396828627123e-17 	# Coeff 5 of forward polynomial 1
                                                          PF1 = 1 	# Power of i/p 1 for coeff 1 of fwd poly 1
                                                          PF2 = 2 	# Power of i/p 1 for coeff 2 of fwd poly 1
                                                          PF3 = 3 	# Power of i/p 1 for coeff 3 of fwd poly 1
                                                          PF4 = 4 	# Power of i/p 1 for coeff 4 of fwd poly 1
                                                          PF5 = 5 	# Power of i/p 1 for coeff 5 of fwd poly 1
                                                          IterInv = 1 	# Use an iterative inverse
                                                          NiterInv = 30 	# Max number of iterations for iterative inverse
                                                          TolInv = 1e-08 	# Target relative error for iterative inverse
                                                       End PolyMap
                                                 End CmpMap
                                              MapB = 	# Second component Mapping
                                                 Begin CmpMap 	# Compound Mapping
                                                    Nin = 3 	# Number of input coordinates
                                                    Nout = 2 	# Number of output coordinates
                                                 IsA Mapping 	# Mapping between coordinate systems
                                                    InvA = 1 	# First Mapping used in inverse direction
                                                    InvB = 1 	# Second Mapping used in inverse direction
                                                    MapA = 	# First component Mapping
                                                       Begin UnitNormMap 	# Compute unit vector and norm relative to a centre
                                                          Nin = 2 	# Number of input coordinates
                                                          Nout = 3 	# Number of output coordinates
                                                       IsA Mapping 	# Mapping between coordinate systems
                                                       End UnitNormMap
                                                    MapB = 	# Second component Mapping
                                                       Begin WinMap 	# Map one window on to another
                                                          Nin = 2 	# Number of input coordinates
                                                       IsA Mapping 	# Mapping between coordinate systems
                                                          Sft1 = -20.365000000000002 	# Shift for axis 1
                                                          Scl1 = 0.01 	# Scale factor for axis 1
                                                          Sft2 = -62.254999999999995 	# Shift for axis 2
                                                          Scl2 = 0.01 	# Scale factor for axis 2
                                                       End WinMap
                                                 End CmpMap
                                           End CmpMap
                                     End CmpMap
                               End CmpMap
                         End CmpMap
                   End CmpMap
             End CmpMap
       End CmpMap
 End CmpMap
""")

patch_position = [17365.72735521638, -0.030469189292489318]
patch_to_sky_matrix = np.array(
    [[patch_to_sky.rate(patch_position, 1, 1), patch_to_sky.rate(patch_position, 1, 2)],
     [patch_to_sky.rate(patch_position, 2, 1), patch_to_sky.rate(patch_position, 2, 2)]],
)
sky_position = patch_to_sky.applyForward(patch_position)
sky_to_ccd_matrix = np.array(
    [[sky_to_ccd.rate(sky_position, 1, 1), sky_to_ccd.rate(sky_position, 1, 2)],
     [sky_to_ccd.rate(sky_position, 2, 1), sky_to_ccd.rate(sky_position, 2, 2)]],
)
print("Expected matrix (from matrix product of separate linearizations):")
print(np.dot(sky_to_ccd_matrix, patch_to_sky_matrix))
patch_to_ccd = patch_to_sky.then(sky_to_ccd)
patch_to_ccd_matrix = np.array(
    [[patch_to_ccd.rate(patch_position, 1, 1), patch_to_ccd.rate(patch_position, 1, 2)],
     [patch_to_ccd.rate(patch_position, 2, 1), patch_to_ccd.rate(patch_position, 2, 2)]],
)
print("Actual matrix from linearization of composite mapping:")
print(patch_to_ccd_matrix)

[timj]

I believe this is with version 9.2.11 (but it doesn't look like 9.2.12 has had any changes in this area).

[timj]

@dsberry this is the output of the test script:

Expected matrix (from matrix product of separate linearizations):
[[ 4.43891233e-01 -5.33461802e-04]
 [-8.94236561e-01  1.07803087e-03]]
Actual matrix from linearization of composite mapping:
[[ 4.43891233e-01  1.05236724e-17]
 [-8.94236560e-01  5.11801974e-18]]

[dsberry]

I've been able to reproduce this. I'll look into a fix some time over the next week or so. Sorry it can't be a bit sooner but I'm in the thick of things at the moment.

[dsberry]

Hopefully the commit I've just made will fix the issue.

[timj]

82444f6

[timj]

@dsberry would it be possible to make a new release with this fix?

[dsberry]

Yes, but it will be next week before I can get round to it.

[dsberry]

@timj New release 9.2.13 now available on github.