This Claude-generated issue was created based on findings from #73 where it was discovered that astLinearApprox was reporting that a spline map was linear even though it had deviations that were larger than the requested tolerance.
Summary
astLinearApprox (src/mapping.c:6491) decides linearity from a small, fixed set of sample positions expressed as constant fractions of the supplied box. Any structure in the Mapping that falls between those positions is invisible, however large it is and however well chosen the box is. Because the positions are fixed fractions rather than randomised, the failure is systematic and reproducible rather than a matter of luck: the same relative geometry defeats the test every time.
For a 2-input Mapping the function uses 17 positions in total:
- 4 fit positions, at the centres of the box faces (
src/mapping.c:6690). The returned fit is exact at these by construction.
- 13 test positions (
npoint = 1 + 2*(ndim_in + 2^ndim_in), src/mapping.c:6795 onward): one near the centre, one half way from the centre to each face centre, the four vertices, and one half way from the centre to each vertex.
In normalised box coordinates those 17 positions are:
fit: (0,0.5) (1,0.5) (0.5,0) (0.5,1)
test: (0.51,0.51)
(0.25,0.52) (0.76,0.52) (0.52,0.25) (0.52,0.76)
(0,0) (1,0) (0,1) (1,1)
(0.245,0.245) (0.765,0.245) (0.245,0.765) (0.765,0.765)
The gaps around (0.35, 0.35) and (0.35, 0.65) are roughly 0.15 of the box wide, which is comfortably large enough to hide a real feature.
The comment at src/mapping.c:6845 says the offsets "introduce some small random offsets to break any regularity in the grid", but the offsets are the compile-time constants 0.48/0.52/0.49/0.51. Nothing is randomised, so a caller cannot work around a miss by retrying.
Reproducer
A bi-cubic SplineMap with 30 coefficients per axis — the resolution a real astrometric fit uses — set to the identity via the Greville abscissae, with a single coefficient moved. A single-coefficient bump has support k/(n-k+1) = 4/27 = 14.8% of the box, which fits in the gap.
#include <stdio.h>
#include <math.h>
#include "ast.h"
#define KB 4
#define NB 30
#define TB (NB+KB)
#define L 1000.0
int main( void ) {
double t[ TB ], cu[ NB*NB ], cv[ NB*NB ];
double lbnd[ 2 ] = { 0.0, 0.0 }, ubnd[ 2 ] = { L, L };
double fit[ 6 ], tol = 0.1, dmax = 0.0;
AstSplineMap *sm;
int i, j, ok;
astBegin;
/* Clamped uniform knots. */
for( i = 0; i < KB; i++ ) t[ i ] = 0.0;
for( i = KB; i < NB; i++ ) t[ i ] = L*( i - KB + 1 )/(double)( NB - KB + 1 );
for( i = NB; i < TB; i++ ) t[ i ] = L;
/* Greville abscissae reproduce u = x and v = y exactly. */
#define GREV(k) ( ( t[(k)+1] + t[(k)+2] + t[(k)+3] )/3.0 )
for( i = 0; i < NB; i++ ) {
for( j = 0; j < NB; j++ ) {
cu[ i*NB + j ] = GREV(i);
cv[ i*NB + j ] = GREV(j);
}
}
/* Move one coefficient, near (0.33, 0.67) of the box. */
cu[ 10*NB + 19 ] += 10.0;
sm = astSplineMap( KB, KB, NB, NB, t, t, cu, cv, "OutUnit=1" );
ok = astLinearApprox( sm, lbnd, ubnd, tol, fit );
printf( "astLinearApprox -> %d, fit = %g %g / %g %g %g %g\n",
ok, fit[0], fit[1], fit[2], fit[3], fit[4], fit[5] );
for( j = 0; j <= 200; j++ ) {
double gx[ 201 ], gy[ 201 ], ax[ 201 ], ay[ 201 ], d;
for( i = 0; i <= 200; i++ ) { gx[i] = L*i/200.0; gy[i] = L*j/200.0; }
astTran2( (AstMapping *) sm, 201, gx, gy, 1, ax, ay );
for( i = 0; i <= 200; i++ ) {
d = hypot( ax[i] - ( fit[0] + fit[2]*gx[i] + fit[3]*gy[i] ),
ay[i] - ( fit[1] + fit[4]*gx[i] + fit[5]*gy[i] ) );
if( d > dmax ) dmax = d;
}
}
printf( "worst departure from that fit = %.4g = %.4g x tol\n",
dmax, dmax/tol );
astEnd;
return 0;
}Output:
astLinearApprox -> 1, fit = 2.70894e-14 0 / 1 5.68434e-17 0 1
worst departure from that fit = 4.418 = 44.18 x tol
The returned fit is bit-identical to the un-bumped identity case. The perturbation is not merely under-weighted, it is not seen at all, so the error is unbounded — it scales linearly with the coefficient:
| coefficient offset |
worst departure |
multiple of tol |
| 0.5 |
0.2209 |
2.2 |
| 2 |
0.8836 |
8.8 |
| 10 |
4.418 |
44 |
| 50 |
22.09 |
221 |
OutUnit=1 is used so the SplineMap behaves as the identity outside its knots; it is not required to produce the effect, it just makes the reproducer's arithmetic clean.
How much sampling would be needed
Tested with the same bump at three positions, checking the residual against the returned fit on an N×N grid:
| grid |
spacing (fraction of box) |
caught |
| 4x4 |
0.333 |
2 of 3 (coincidence) |
| 8x8 |
0.143 |
1 of 3 |
| 16x16 |
0.067 |
3 of 3 |
| 32x32 |
0.032 |
3 of 3 |
The spacing has to be comfortably finer than the smallest structure the Mapping can hold. 8×8 is about the same size as the bump and behaves like a coin toss.
Cost is not the obstacle — sampling that SplineMap:
| grid |
time |
| 8x8 (64 points) |
0.015 ms |
| 16x16 (256 points) |
0.062 ms |
| 32x32 (1024 points) |
0.257 ms |
For scale, one complete astWrite of a real LSST pixels-to-sky FrameSet through a FITS-WCS FitsChan measures 0.708 ms on the same machine, and FitLine in fitschan.c already uses NP = 101 positions per axis.
Documentation
The parameter description for tol is technically accurate — it says a fit is rejected if it deviates "by more than this amount at any point which is tested". But the Description above it reads as a claim about the whole range:
This function tests the forward coordinate transformation implemented by a
Mapping over a given range of input coordinates. If the transformation is
found to be linear to a specified level of accuracy...
A caller reading the Description would reasonably assume the box is covered. At minimum the Description should state that a fixed, sparse set of positions is sampled and that structure smaller than the spacing between them will not be detected.
Why this is not hypothetical
This was found through SIPIntWorld in fitschan.c, which used astLinearApprox to decide whether a Mapping could be described by the FITS-WCS SIP conventions and took the CDi_j values from the returned fit. A real LSST difference-image WCS containing a SplineMap was accepted as linear, and the CD matrix was taken from a fit that had never seen the spline. Depending on a second coincidence the write then either segfaulted or emitted a header with SIP coefficients and no CTYPE, CRPIX or CD cards at all, reported as success.
fitschan.c has since been changed to do its own least-squares fit over a 16×16 grid of positions spanning the image, so AST's FITS writer no longer depends on this behaviour. The general issue in astLinearApprox remains for other callers.
Options
Listed with their blast radius, because that is the deciding factor here.
-
Documentation only. State in the Description that a fixed sparse sample is used and that sub-sample-scale structure is not detected. No behaviour change, no risk. Does not help callers who need the guarantee.
-
Increase the number of test positions, leaving the fit as it is. This errs toward reporting "not linear", so the adaptive callers subdivide more: more accurate, somewhat slower, and resampled output changes in the direction of correctness. Still a behaviour change for every caller.
-
Add a sample-count or density parameter, or a sibling entry point. Callers who need the guarantee opt in; existing callers are untouched. This is the option I would suggest — it is the only one that fixes the problem without moving anyone else's pixels.
-
Change the fit itself to least squares over a denser grid. Note that for a genuinely linear Mapping the current 4-point fit is already exact, so this changes nothing in the case where the function answers "yes". It only alters the coefficients returned for Mappings that are not linear, where the useful answer is the rejection rather than the plane. Against that limited benefit, the internal callers are ResampleAdaptively (src/mapping.c:14815), RebinAdaptively (src/mapping.c:10606) and TranGridAdaptively (src/mapping.c:21076) — every adaptive subdivision in AST — so resampled and rebinned pixel values would change throughout Starlink. It is also public API in both C and Fortran (src/fmapping.c). I would not recommend this without a deliberate decision to accept that.
Whichever route is taken, the fixed 0.48/0.52/0.49/0.51 offsets should either be genuinely randomised as their comment claims, or the comment should be corrected.
This Claude-generated issue was created based on findings from #73 where it was discovered that
astLinearApproxwas reporting that a spline map was linear even though it had deviations that were larger than the requested tolerance.Summary
astLinearApprox(src/mapping.c:6491) decides linearity from a small, fixed set of sample positions expressed as constant fractions of the supplied box. Any structure in the Mapping that falls between those positions is invisible, however large it is and however well chosen the box is. Because the positions are fixed fractions rather than randomised, the failure is systematic and reproducible rather than a matter of luck: the same relative geometry defeats the test every time.For a 2-input Mapping the function uses 17 positions in total:
src/mapping.c:6690). The returned fit is exact at these by construction.npoint = 1 + 2*(ndim_in + 2^ndim_in),src/mapping.c:6795onward): one near the centre, one half way from the centre to each face centre, the four vertices, and one half way from the centre to each vertex.In normalised box coordinates those 17 positions are:
The gaps around
(0.35, 0.35)and(0.35, 0.65)are roughly 0.15 of the box wide, which is comfortably large enough to hide a real feature.The comment at
src/mapping.c:6845says the offsets "introduce some small random offsets to break any regularity in the grid", but the offsets are the compile-time constants0.48/0.52/0.49/0.51. Nothing is randomised, so a caller cannot work around a miss by retrying.Reproducer
A bi-cubic
SplineMapwith 30 coefficients per axis — the resolution a real astrometric fit uses — set to the identity via the Greville abscissae, with a single coefficient moved. A single-coefficient bump has supportk/(n-k+1)=4/27= 14.8% of the box, which fits in the gap.Output:
The returned fit is bit-identical to the un-bumped identity case. The perturbation is not merely under-weighted, it is not seen at all, so the error is unbounded — it scales linearly with the coefficient:
tolOutUnit=1is used so the SplineMap behaves as the identity outside its knots; it is not required to produce the effect, it just makes the reproducer's arithmetic clean.How much sampling would be needed
Tested with the same bump at three positions, checking the residual against the returned fit on an N×N grid:
The spacing has to be comfortably finer than the smallest structure the Mapping can hold. 8×8 is about the same size as the bump and behaves like a coin toss.
Cost is not the obstacle — sampling that SplineMap:
For scale, one complete astWrite of a real LSST pixels-to-sky FrameSet through a FITS-WCS FitsChan measures 0.708 ms on the same machine, and FitLine in fitschan.c already uses NP = 101 positions per axis.
Documentation
The parameter description for tol is technically accurate — it says a fit is rejected if it deviates "by more than this amount at any point which is tested". But the Description above it reads as a claim about the whole range:
A caller reading the Description would reasonably assume the box is covered. At minimum the Description should state that a fixed, sparse set of positions is sampled and that structure smaller than the spacing between them will not be detected.
Why this is not hypothetical
This was found through SIPIntWorld in fitschan.c, which used astLinearApprox to decide whether a Mapping could be described by the FITS-WCS SIP conventions and took the CDi_j values from the returned fit. A real LSST difference-image WCS containing a SplineMap was accepted as linear, and the CD matrix was taken from a fit that had never seen the spline. Depending on a second coincidence the write then either segfaulted or emitted a header with SIP coefficients and no CTYPE, CRPIX or CD cards at all, reported as success.
fitschan.c has since been changed to do its own least-squares fit over a 16×16 grid of positions spanning the image, so AST's FITS writer no longer depends on this behaviour. The general issue in astLinearApprox remains for other callers.
Options
Listed with their blast radius, because that is the deciding factor here.
Documentation only. State in the Description that a fixed sparse sample is used and that sub-sample-scale structure is not detected. No behaviour change, no risk. Does not help callers who need the guarantee.
Increase the number of test positions, leaving the fit as it is. This errs toward reporting "not linear", so the adaptive callers subdivide more: more accurate, somewhat slower, and resampled output changes in the direction of correctness. Still a behaviour change for every caller.
Add a sample-count or density parameter, or a sibling entry point. Callers who need the guarantee opt in; existing callers are untouched. This is the option I would suggest — it is the only one that fixes the problem without moving anyone else's pixels.
Change the fit itself to least squares over a denser grid. Note that for a genuinely linear Mapping the current 4-point fit is already exact, so this changes nothing in the case where the function answers "yes". It only alters the coefficients returned for Mappings that are not linear, where the useful answer is the rejection rather than the plane. Against that limited benefit, the internal callers are ResampleAdaptively (src/mapping.c:14815), RebinAdaptively (src/mapping.c:10606) and TranGridAdaptively (src/mapping.c:21076) — every adaptive subdivision in AST — so resampled and rebinned pixel values would change throughout Starlink. It is also public API in both C and Fortran (src/fmapping.c). I would not recommend this without a deliberate decision to accept that.
Whichever route is taken, the fixed 0.48/0.52/0.49/0.51 offsets should either be genuinely randomised as their comment claims, or the comment should be corrected.