From 52339816b1ef9b3404fa7d91473d6e8368a9949d Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 11:36:27 -0400 Subject: [PATCH 01/20] move polynomial fitting to mapcat --- mapcat/database/pointing_residual.py | 3 +- mapcat/pointing/const.py | 4 +- mapcat/pointing/poly.py | 182 +++++++++++++++++++++++++++ 3 files changed, 185 insertions(+), 4 deletions(-) create mode 100644 mapcat/pointing/poly.py diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index 79ba2ec..dd3c6fa 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -40,7 +40,6 @@ class PointingResidualTable(SQLModel, table=True): residual_model: ConstantPointingModel = Field( discriminator="model_type", sa_type=JSONEncodedPydantic(ConstantPointingModel) ) - residual_stats: PointingModelStats = Field( - default_factory=PointingModelStats, sa_type=JSONEncodedPydantic(PointingModelStats) + residual_stats: PointingModelStats | None = Field( nullable=True, sa_type=JSONEncodedPydantic(PointingModelStats) ) map: DepthOneMapTable = Relationship(back_populates="pointing_residual") diff --git a/mapcat/pointing/const.py b/mapcat/pointing/const.py index 81010d5..2a9ca91 100644 --- a/mapcat/pointing/const.py +++ b/mapcat/pointing/const.py @@ -17,7 +17,7 @@ class ConstantPointingModel(PointingModelProtocol): dec_offset: AstroPydanticQuantity[u.deg] def predict(self, pos: SkyCoord) -> SkyCoord: - ra = pos.ra + self.ra_offset - dec = pos.dec + self.dec_offset + ra = pos.ra - self.ra_offset + dec = pos.dec - self.dec_offset return SkyCoord(ra=ra, dec=dec, frame=pos.frame) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py new file mode 100644 index 0000000..dab6e39 --- /dev/null +++ b/mapcat/pointing/poly.py @@ -0,0 +1,182 @@ +""" +Polynomial pointing model. +""" +import numpy as np +from typing import Literal + +from astropy.coordinates import SkyCoord +from astropy import units as u +from astropydantic import AstroPydanticQuantity +from mapcat.pointing.base import PointingModelProtocol, PointingModelStats + + +class PolynomialCoefficients(AstroPydanticQuantity): + """ + Coefficients for a polynomial pointing model. + for example a 2D polynomial of order 2, + coeffs={'x^2':1, + 'y^2':1, + 'xy':1, + 'y': 2, + 'x':2, + 'constant':3 + } + labels = {'x':'ra', 'y':'dec'} + """ + coeffs: dict[str,float] + labels: dict[str, str] + poly_order: int + + + +class PolynomialPointingModel(PointingModelProtocol): + model_type: Literal["polynomial"] = "polynomial" + + poly_order: int + ra_model_coefficients: PolynomialCoefficients | None = None + dec_model_coefficients: PolynomialCoefficients | None = None + + ## Basis terms for 2D polynomial fit + def _poly_terms(self, x, y): + terms = [] + for i in range(self.poly_order + 1): + for j in range(self.poly_order + 1 - i): + terms.append((x**i) * (y**j)) + return np.vstack(terms).T + + def _poly_keys(self): + keys = [] + for i in range(self.poly_order + 1): + for j in range(self.poly_order + 1 - i): + keys.append(f"x^{i}y^{j}") + return keys + + def calculate_model(self, + measured_positions: SkyCoord, + expected_positions: SkyCoord, + position_uncertainty: tuple[u.Quantity, u.Quantity] | u.Quantity | None = None, + weights: tuple[list[float], list[float]] | list[float] | None = None + ) -> SkyCoord: + """ + Calculate the polynomial coefficients for the pointing model + using the measured and expected positions. + + Optionally accept uncertainties in the measured positions + and measurement weights to perform a weighted fit. + + position_uncertainty can be provided as a tuple of (ra_uncertainties, dec_uncertainties) + or a single quantity that applies to both. + similarly, weights can be provided as a tuple of (ra_weights, dec_weights) + or a single list that applies to both. + + weights are resolved from uncertainties if not provided, using inverse variance weighting. + """ + # Calculate offsets + ra_offsets = measured_positions.ra - expected_positions.ra + dec_offsets = measured_positions.dec - expected_positions.dec + n = len(ra_offsets) + if n == 0: + raise ValueError("No positions provided for model calculation.") + + # Unpack position_uncertainty into ra_uncerts, dec_uncerts + if isinstance(position_uncertainty, tuple): + ra_uncerts, dec_uncerts = position_uncertainty + else: + ra_uncerts = dec_uncerts = position_uncertainty # None or single Quantity applied to both + + # Unpack weights into ra_weights, dec_weights + if isinstance(weights, tuple): + ra_weights, dec_weights = weights + else: + ra_weights = dec_weights = weights # None or single list applied to both + + # Lots of logic to check if weights exist, etc. + # Resolve weights from uncertainties if not provided + if ra_weights is None and dec_weights is None: + if ra_uncerts is not None and dec_uncerts is not None: + ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts**2) + elif ra_uncerts is not None: + dec_uncerts = np.ones(n) * np.mean(ra_uncerts) + ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts**2) + elif dec_uncerts is not None: + ra_uncerts = np.ones(n) * np.mean(dec_uncerts) + ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts**2) + else: + # No uncertainty or weights provided — uniform + ra_weights = dec_weights = np.ones(n) + elif ra_weights is None: + ra_weights = dec_weights + elif dec_weights is None: + dec_weights = ra_weights + + ras = measured_positions.ra.to_value(u.deg) + decs = measured_positions.dec.to_value(u.deg) + ## RA polynomial fit + A_ra = self._poly_terms(ras, decs) + y_ra = ra_offsets.to_value(u.deg) + w_ra = ra_weights + + ## Apply weights + Aw = A_ra * w_ra[:, None] + yw = y_ra * w_ra + coeffs_ra, *_ = np.linalg.lstsq(Aw, yw, rcond=None) + + ## Dec polynomial fit + A_dec = self._poly_terms(ras, decs) + y_dec = dec_offsets.to_value(u.deg) + w_dec = dec_weights + + Aw = A_dec * w_dec[:, None] + yw = y_dec * w_dec + coeffs_dec, *_ = np.linalg.lstsq(Aw, yw, rcond=None) + + ra_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_ra)} + dec_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_dec)} + + self.ra_model_coefficients=PolynomialCoefficients(coeffs=ra_coeff_dict, labels={'x':'ra', 'y':'dec'}, poly_order=self.poly_order) + self.dec_model_coefficients=PolynomialCoefficients(coeffs=dec_coeff_dict, labels={'x':'ra', 'y':'dec'}, poly_order=self.poly_order) + + + def model_fn( + self, x: u.Quantity, y: u.Quantity, coeffs: np.ndarray + ) -> u.Quantity: + x = np.atleast_1d(x.to_value(u.deg)) + y = np.atleast_1d(y.to_value(u.deg)) + T = self._poly_terms(x, y) + return (T @ coeffs) * u.deg + + def extract_coefficients(self) -> np.ndarray: + ra_coeff_array = np.zeros(len(self.ra_model_coefficients.coeffs)) + for i, key in enumerate(self._poly_keys()): + ra_coeff_array[i] = self.ra_model_coefficients.coeffs.get(key, 0) + + dec_coeff_array = np.zeros(len(self.dec_model_coefficients.coeffs)) + for i, key in enumerate(self._poly_keys()): + dec_coeff_array[i] = self.dec_model_coefficients.coeffs.get(key, 0) + + + return ra_coeff_array, dec_coeff_array + + def predict(self, pos: SkyCoord) -> SkyCoord: + racoeffs,deccoeffs = self.extract_coefficients() + ra_offset = self.model_fn(pos.ra, pos.dec, racoeffs)[0] + dec_offset = self.model_fn(pos.ra, pos.dec, deccoeffs)[0] + ra = pos.ra - ra_offset*u.deg + dec = pos.dec - dec_offset*u.deg + + return SkyCoord(ra=ra, dec=dec, frame=pos.frame) + + def calculate_statistics(self, positions: SkyCoord): + new_positions = self.predict(positions) + ra_residuals = (new_positions.ra - positions.ra).to(u.arcsec) + dec_residuals = (new_positions.dec - positions.dec).to(u.arcsec) + mean_ra = np.mean(ra_residuals) + mean_dec = np.mean(dec_residuals) + std_ra = np.std(ra_residuals) + std_dec = np.std(dec_residuals) + return PointingModelStats( + mean_ra_offset=mean_ra, + mean_dec_offset=mean_dec, + stddev_ra_offset=std_ra, + stddev_dec_offset=std_dec, + ) \ No newline at end of file From 2676d8bf459e8732d22b38bcd24ef0088d0f5fa7 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 12:12:02 -0400 Subject: [PATCH 02/20] polynomialcoefficients is basemodel. not sure why I had astropydantic quantity there --- mapcat/pointing/poly.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index dab6e39..dc9c460 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -2,15 +2,15 @@ Polynomial pointing model. """ import numpy as np +from pydantic import BaseModel from typing import Literal from astropy.coordinates import SkyCoord from astropy import units as u -from astropydantic import AstroPydanticQuantity from mapcat.pointing.base import PointingModelProtocol, PointingModelStats -class PolynomialCoefficients(AstroPydanticQuantity): +class PolynomialCoefficients(BaseModel): """ Coefficients for a polynomial pointing model. for example a 2D polynomial of order 2, From f63815fa82730970db151f4172e599d2f490cb04 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 12:23:49 -0400 Subject: [PATCH 03/20] unit fix --- mapcat/pointing/poly.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index dc9c460..7a102f9 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -161,8 +161,8 @@ def predict(self, pos: SkyCoord) -> SkyCoord: racoeffs,deccoeffs = self.extract_coefficients() ra_offset = self.model_fn(pos.ra, pos.dec, racoeffs)[0] dec_offset = self.model_fn(pos.ra, pos.dec, deccoeffs)[0] - ra = pos.ra - ra_offset*u.deg - dec = pos.dec - dec_offset*u.deg + ra = pos.ra - ra_offset + dec = pos.dec - dec_offset return SkyCoord(ra=ra, dec=dec, frame=pos.frame) From 1a79df4b78f56f68559696a4c7ea0afee0aacdb1 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 13:25:17 -0400 Subject: [PATCH 04/20] allow residual_model to be constant pointing or polynomial, descriminated by model_type --- mapcat/database/pointing_residual.py | 10 ++++++++-- mapcat/pointing/poly.py | 5 +++-- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index dd3c6fa..c2fc261 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -2,15 +2,22 @@ Table containing pointing residuals. """ +from typing import Annotated + from sqlmodel import Field, Relationship, SQLModel from mapcat.pointing.const import ConstantPointingModel +from mapcat.pointing.poly import PolynomialPointingModel from mapcat.pointing.base import PointingModelStats from .depth_one_map import DepthOneMapTable from .json import JSONEncodedPydantic +PointingModel = Annotated[ + ConstantPointingModel | PolynomialPointingModel, + Field(discriminator="model_type") +] class PointingResidualTable(SQLModel, table=True): """ @@ -37,8 +44,7 @@ class PointingResidualTable(SQLModel, table=True): foreign_key="depth_one_maps.map_id", ondelete="CASCADE", ) - residual_model: ConstantPointingModel = Field( - discriminator="model_type", sa_type=JSONEncodedPydantic(ConstantPointingModel) + residual_model: PointingModel = Field(sa_type=JSONEncodedPydantic(PointingModel) ) residual_stats: PointingModelStats | None = Field( nullable=True, sa_type=JSONEncodedPydantic(PointingModelStats) ) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 7a102f9..72d9861 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -25,6 +25,7 @@ class PolynomialCoefficients(BaseModel): """ coeffs: dict[str,float] labels: dict[str, str] + unit: u.Unit = u.deg poly_order: int @@ -133,8 +134,8 @@ def calculate_model(self, ra_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_ra)} dec_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_dec)} - self.ra_model_coefficients=PolynomialCoefficients(coeffs=ra_coeff_dict, labels={'x':'ra', 'y':'dec'}, poly_order=self.poly_order) - self.dec_model_coefficients=PolynomialCoefficients(coeffs=dec_coeff_dict, labels={'x':'ra', 'y':'dec'}, poly_order=self.poly_order) + self.ra_model_coefficients=PolynomialCoefficients(coeffs=ra_coeff_dict, labels={'x':'ra', 'y':'dec'}, unit=u.deg, poly_order=self.poly_order) + self.dec_model_coefficients=PolynomialCoefficients(coeffs=dec_coeff_dict, labels={'x':'ra', 'y':'dec'}, unit=u.deg, poly_order=self.poly_order) def model_fn( From b969befdcfc720fc39add4eca698bf072a9d4d94 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 13:30:48 -0400 Subject: [PATCH 05/20] astropy units not kosher for pydantic --- mapcat/pointing/poly.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 72d9861..82f6965 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -1,6 +1,7 @@ """ Polynomial pointing model. """ +from astropydantic import AstropydanticUnit import numpy as np from pydantic import BaseModel from typing import Literal @@ -25,7 +26,7 @@ class PolynomialCoefficients(BaseModel): """ coeffs: dict[str,float] labels: dict[str, str] - unit: u.Unit = u.deg + unit: AstropydanticUnit = AstropydanticUnit(u.deg) poly_order: int From 4f68fab71f3d0600fafb76cf7e5401336d76444b Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 13:31:39 -0400 Subject: [PATCH 06/20] typo --- mapcat/pointing/poly.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 82f6965..1d7d64f 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -1,7 +1,7 @@ """ Polynomial pointing model. """ -from astropydantic import AstropydanticUnit +from astropydantic import AstroPydanticUnit import numpy as np from pydantic import BaseModel from typing import Literal @@ -26,7 +26,7 @@ class PolynomialCoefficients(BaseModel): """ coeffs: dict[str,float] labels: dict[str, str] - unit: AstropydanticUnit = AstropydanticUnit(u.deg) + unit: AstroPydanticUnit = AstroPydanticUnit(u.deg) poly_order: int From a2d05a75e714b11db36b24df8b7db6bef0535b43 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 13:33:16 -0400 Subject: [PATCH 07/20] default deg unit --- mapcat/pointing/poly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 1d7d64f..5a04e15 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -26,7 +26,7 @@ class PolynomialCoefficients(BaseModel): """ coeffs: dict[str,float] labels: dict[str, str] - unit: AstroPydanticUnit = AstroPydanticUnit(u.deg) + unit: AstroPydanticUnit = u.deg poly_order: int From d53ee74870b71e7a1026aa1b1574260b332f235b Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 13:57:32 -0400 Subject: [PATCH 08/20] fix the union of the two pointing models --- mapcat/database/pointing_residual.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index c2fc261..dda3079 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -2,7 +2,7 @@ Table containing pointing residuals. """ -from typing import Annotated +from typing import Annotated, Union from sqlmodel import Field, Relationship, SQLModel @@ -14,8 +14,7 @@ from .depth_one_map import DepthOneMapTable from .json import JSONEncodedPydantic -PointingModel = Annotated[ - ConstantPointingModel | PolynomialPointingModel, +PointingModel = Annotated[Union[ConstantPointingModel, PolynomialPointingModel], Field(discriminator="model_type") ] From a1cbca883f1d4c9a1da9c90b5eeae4467e6fd124 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 14:09:28 -0400 Subject: [PATCH 09/20] allow for union which wraps the actual pydantic classes --- mapcat/database/json.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/mapcat/database/json.py b/mapcat/database/json.py index bcf9cee..f99fdad 100644 --- a/mapcat/database/json.py +++ b/mapcat/database/json.py @@ -4,6 +4,8 @@ See: https://github.com/fastapi/sqlmodel/pull/1324 - this can be removed at some point. """ +import typing + from sqlalchemy.types import JSON, TypeDecorator @@ -23,4 +25,8 @@ def process_bind_param(self, value, dialect): def process_result_value(self, value, dialect): if value is None: return None - return self.pydantic_class.model_validate(value) + cls = self.pydantic_class + args = typing.get_args(cls) + if args: + cls = next(a for a in args if a is not type(None)) + return cls.model_validate(value) From 0719fdcdc804a22a232c1e912429da53608fb915 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 14:15:37 -0400 Subject: [PATCH 10/20] remove annotate --- mapcat/database/pointing_residual.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index dda3079..749a18b 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -2,7 +2,7 @@ Table containing pointing residuals. """ -from typing import Annotated, Union +from typing import Union from sqlmodel import Field, Relationship, SQLModel @@ -14,9 +14,7 @@ from .depth_one_map import DepthOneMapTable from .json import JSONEncodedPydantic -PointingModel = Annotated[Union[ConstantPointingModel, PolynomialPointingModel], - Field(discriminator="model_type") -] +PointingModel = Union[ConstantPointingModel, PolynomialPointingModel] class PointingResidualTable(SQLModel, table=True): """ From b3a9a796393fe39d8f73dc4797dbac684b33ad9f Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Tue, 28 Apr 2026 14:23:03 -0400 Subject: [PATCH 11/20] ok, it needs the type adaptor --- mapcat/database/json.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/mapcat/database/json.py b/mapcat/database/json.py index f99fdad..0c6f2b8 100644 --- a/mapcat/database/json.py +++ b/mapcat/database/json.py @@ -4,8 +4,7 @@ See: https://github.com/fastapi/sqlmodel/pull/1324 - this can be removed at some point. """ -import typing - +from pydantic import TypeAdapter from sqlalchemy.types import JSON, TypeDecorator @@ -16,6 +15,7 @@ class JSONEncodedPydantic(TypeDecorator): def __init__(self, pydantic_class, *args, **kwargs): super().__init__(*args, **kwargs) self.pydantic_class = pydantic_class + self._adapter = TypeAdapter(pydantic_class) def process_bind_param(self, value, dialect): if value is None: @@ -25,8 +25,4 @@ def process_bind_param(self, value, dialect): def process_result_value(self, value, dialect): if value is None: return None - cls = self.pydantic_class - args = typing.get_args(cls) - if args: - cls = next(a for a in args if a is not type(None)) - return cls.model_validate(value) + return self._adapter.validate_python(value) From 99d825a51b7d4221a098f3d149067add59c71f4a Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Wed, 29 Apr 2026 10:58:57 -0400 Subject: [PATCH 12/20] add simple tests for pointing models --- tests/test_pointing.py | 47 ++++++++++++++++++++++++++++++++++++++---- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/tests/test_pointing.py b/tests/test_pointing.py index 0768218..44fb0f3 100644 --- a/tests/test_pointing.py +++ b/tests/test_pointing.py @@ -2,15 +2,19 @@ Tests for the pointing residual models. """ -from astropy.units import deg +from astropy import units as u +from astropy.coordinates import SkyCoord +import numpy as np from sqlmodel import select from mapcat.database import DepthOneMapTable, PointingResidualTable from mapcat.pointing.const import ConstantPointingModel +from mapcat.pointing.poly import PolynomialPointingModel + def test_add_retrieve_pointing(database_sessionmaker): - model = ConstantPointingModel(ra_offset=0.5 * deg, dec_offset=0.5 * deg) + model = ConstantPointingModel(ra_offset=0.5 * u.deg, dec_offset=0.5 * u.deg) with database_sessionmaker() as session: sample_map = DepthOneMapTable( @@ -42,5 +46,40 @@ def test_add_retrieve_pointing(database_sessionmaker): recovered_pointing = recovered_map.pointing_residual[0] residual_model = recovered_pointing.residual_model - assert residual_model.ra_offset > 0.4 * deg - assert residual_model.dec_offset > 0.4 * deg + assert residual_model.ra_offset > 0.4 * u.deg + assert residual_model.dec_offset > 0.4 * u.deg + + +def test_make_constant_pointing_model(): + model = ConstantPointingModel(ra_offset=0.5 * u.deg, dec_offset=0.5 * u.deg) + + assert model.ra_offset == 0.5 * u.deg + assert model.dec_offset == 0.5 * u.deg + + og_pos = SkyCoord(ra=10 * u.deg, dec=20 * u.deg + ) + offset_pos = SkyCoord(og_pos.ra + model.ra_offset, og_pos.dec + model.dec_offset) + new_pos = model.predict(offset_pos) + assert new_pos.ra == og_pos.ra + assert new_pos.dec == og_pos.dec + + +def test_make_polynomial_pointing_model(): + model = PolynomialPointingModel(poly_order=2) + ras = np.arange(0, 10, 1) * u.deg + decs = np.arange(0, 10, 1) * u.deg + offset = 1.0 * u.arcmin + slope = 0.1 * u.arcmin / u.deg + + offset_positions = SkyCoord( + ra=ras + offset + slope * ras, dec=decs + offset + slope * decs + ) + model.calculate_model(measured_positions=offset_positions, expected_positions=SkyCoord(ras, decs)) + + assert model.ra_model_coefficients is not None + assert model.dec_model_coefficients is not None + + for i, offset_pos in enumerate(offset_positions): + predicted_pos = model.predict(offset_pos) + assert np.isclose(predicted_pos.ra.to_value(u.arcmin), ras[i].to_value(u.arcmin), atol=0.1) + assert np.isclose(predicted_pos.dec.to_value(u.arcmin), decs[i].to_value(u.arcmin), atol=0.1) From 0692598597f1bd23a502f4d090dedeb9e8759779 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Wed, 29 Apr 2026 15:06:20 -0400 Subject: [PATCH 13/20] feng shui --- mapcat/database/pointing_residual.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index 749a18b..9bb1f0f 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -43,6 +43,6 @@ class PointingResidualTable(SQLModel, table=True): ) residual_model: PointingModel = Field(sa_type=JSONEncodedPydantic(PointingModel) ) - residual_stats: PointingModelStats | None = Field( nullable=True, sa_type=JSONEncodedPydantic(PointingModelStats) + residual_stats: PointingModelStats | None = Field(nullable=True, default=None, sa_type=JSONEncodedPydantic(PointingModelStats) ) map: DepthOneMapTable = Relationship(back_populates="pointing_residual") From 823219b73f4dbb6ddd36f2bbef2d754bc7834a9a Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Wed, 29 Apr 2026 16:24:51 -0400 Subject: [PATCH 14/20] copilot suggestions, allow predict to work on array of positions --- mapcat/database/pointing_residual.py | 2 +- mapcat/pointing/poly.py | 38 ++++++++++++++++++---------- tests/test_pointing.py | 4 +++ 3 files changed, 30 insertions(+), 14 deletions(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index 9bb1f0f..6ce9082 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -26,7 +26,7 @@ class PointingResidualTable(SQLModel, table=True): ---------- map_id : int Internal ID of the depth one map - residual_model: ConstantPointingModel + residual_model: ConstantPointingModel | PolynomialPointingModel The pointing model to actually store in the database. residual_stats: PointingModelStats Statistics about the pointing residuals, such as mean and stddev of RA and Dec offsets diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 5a04e15..6d5bb67 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -1,13 +1,14 @@ """ Polynomial pointing model. """ + +from astropy.coordinates import SkyCoord +from astropy import units as u from astropydantic import AstroPydanticUnit import numpy as np from pydantic import BaseModel from typing import Literal -from astropy.coordinates import SkyCoord -from astropy import units as u from mapcat.pointing.base import PointingModelProtocol, PointingModelStats @@ -58,9 +59,9 @@ def calculate_model(self, expected_positions: SkyCoord, position_uncertainty: tuple[u.Quantity, u.Quantity] | u.Quantity | None = None, weights: tuple[list[float], list[float]] | list[float] | None = None - ) -> SkyCoord: + ): """ - Calculate the polynomial coefficients for the pointing model + Calculate and set the polynomial coefficients for the pointing model using the measured and expected positions. Optionally accept uncertainties in the measured positions @@ -96,13 +97,13 @@ def calculate_model(self, # Resolve weights from uncertainties if not provided if ra_weights is None and dec_weights is None: if ra_uncerts is not None and dec_uncerts is not None: - ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts**2) + ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts.to_value(u.deg)**2 + dec_uncerts.to_value(u.deg)**2) elif ra_uncerts is not None: - dec_uncerts = np.ones(n) * np.mean(ra_uncerts) - ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts**2) + dec_uncerts = np.ones(n) * np.mean(ra_uncerts.to_value(u.deg)) + ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts.to_value(u.deg)**2 + dec_uncerts**2) elif dec_uncerts is not None: - ra_uncerts = np.ones(n) * np.mean(dec_uncerts) - ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts**2) + ra_uncerts = np.ones(n) * np.mean(dec_uncerts.to_value(u.deg)) + ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts.to_value(u.deg)**2) else: # No uncertainty or weights provided — uniform ra_weights = dec_weights = np.ones(n) @@ -111,6 +112,11 @@ def calculate_model(self, elif dec_weights is None: dec_weights = ra_weights + ra_weights = np.asarray(ra_weights) + dec_weights = np.asarray(dec_weights) + assert len(ra_weights) == n, "Length of ra_weights must match number of positions" + assert len(dec_weights) == n, "Length of dec_weights must match number of positions" + ras = measured_positions.ra.to_value(u.deg) decs = measured_positions.dec.to_value(u.deg) ## RA polynomial fit @@ -147,7 +153,14 @@ def model_fn( T = self._poly_terms(x, y) return (T @ coeffs) * u.deg - def extract_coefficients(self) -> np.ndarray: + def extract_coefficients(self) -> tuple[np.ndarray, np.ndarray]: + """ + Extract the coefficients from the PolynomialCoefficients dataclasss and + return them as arrays in the correct order for the model function. + """ + if self.ra_model_coefficients is None or self.dec_model_coefficients is None: + raise ValueError("Model coefficients have not been calculated yet.") + ra_coeff_array = np.zeros(len(self.ra_model_coefficients.coeffs)) for i, key in enumerate(self._poly_keys()): ra_coeff_array[i] = self.ra_model_coefficients.coeffs.get(key, 0) @@ -156,13 +169,12 @@ def extract_coefficients(self) -> np.ndarray: for i, key in enumerate(self._poly_keys()): dec_coeff_array[i] = self.dec_model_coefficients.coeffs.get(key, 0) - return ra_coeff_array, dec_coeff_array def predict(self, pos: SkyCoord) -> SkyCoord: racoeffs,deccoeffs = self.extract_coefficients() - ra_offset = self.model_fn(pos.ra, pos.dec, racoeffs)[0] - dec_offset = self.model_fn(pos.ra, pos.dec, deccoeffs)[0] + ra_offset = self.model_fn(pos.ra, pos.dec, racoeffs) + dec_offset = self.model_fn(pos.ra, pos.dec, deccoeffs) ra = pos.ra - ra_offset dec = pos.dec - dec_offset diff --git a/tests/test_pointing.py b/tests/test_pointing.py index 44fb0f3..5d55ec1 100644 --- a/tests/test_pointing.py +++ b/tests/test_pointing.py @@ -83,3 +83,7 @@ def test_make_polynomial_pointing_model(): predicted_pos = model.predict(offset_pos) assert np.isclose(predicted_pos.ra.to_value(u.arcmin), ras[i].to_value(u.arcmin), atol=0.1) assert np.isclose(predicted_pos.dec.to_value(u.arcmin), decs[i].to_value(u.arcmin), atol=0.1) + + predicted_pos = model.predict(offset_positions) + assert np.all(np.isclose(predicted_pos.ra.to_value(u.arcmin), ras.to_value(u.arcmin), atol=0.1)) + assert np.all(np.isclose(predicted_pos.dec.to_value(u.arcmin), decs.to_value(u.arcmin), atol=0.1)) From b5eb5f99f2a3110d4742efdbff4f713b37a02915 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Wed, 29 Apr 2026 16:37:24 -0400 Subject: [PATCH 15/20] ruff --- mapcat/database/pointing_residual.py | 3 +-- mapcat/pointing/poly.py | 22 +++++++++++++++++++--- tests/test_pointing.py | 3 +-- 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index 6ce9082..32d01cd 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -6,11 +6,10 @@ from sqlmodel import Field, Relationship, SQLModel +from mapcat.pointing.base import PointingModelStats from mapcat.pointing.const import ConstantPointingModel from mapcat.pointing.poly import PolynomialPointingModel -from mapcat.pointing.base import PointingModelStats - from .depth_one_map import DepthOneMapTable from .json import JSONEncodedPydantic diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 6d5bb67..f6f9d2c 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -2,12 +2,13 @@ Polynomial pointing model. """ -from astropy.coordinates import SkyCoord +from typing import Literal + +import numpy as np from astropy import units as u +from astropy.coordinates import SkyCoord from astropydantic import AstroPydanticUnit -import numpy as np from pydantic import BaseModel -from typing import Literal from mapcat.pointing.base import PointingModelProtocol, PointingModelStats @@ -64,6 +65,7 @@ def calculate_model(self, Calculate and set the polynomial coefficients for the pointing model using the measured and expected positions. + Optionally accept uncertainties in the measured positions and measurement weights to perform a weighted fit. @@ -73,6 +75,15 @@ def calculate_model(self, or a single list that applies to both. weights are resolved from uncertainties if not provided, using inverse variance weighting. + + Raises + ------ + ValueError + If no positions are provided for model calculation. + ValueError + If the lengths of weights do not match the number of positions. + ValueError + If model coefficients have not been calculated yet when extracting coefficients. """ # Calculate offsets ra_offsets = measured_positions.ra - expected_positions.ra @@ -157,6 +168,11 @@ def extract_coefficients(self) -> tuple[np.ndarray, np.ndarray]: """ Extract the coefficients from the PolynomialCoefficients dataclasss and return them as arrays in the correct order for the model function. + + Raises + ------ + ValueError + If model coefficients have not been calculated yet. """ if self.ra_model_coefficients is None or self.dec_model_coefficients is None: raise ValueError("Model coefficients have not been calculated yet.") diff --git a/tests/test_pointing.py b/tests/test_pointing.py index 5d55ec1..a133dc0 100644 --- a/tests/test_pointing.py +++ b/tests/test_pointing.py @@ -2,9 +2,9 @@ Tests for the pointing residual models. """ +import numpy as np from astropy import units as u from astropy.coordinates import SkyCoord -import numpy as np from sqlmodel import select from mapcat.database import DepthOneMapTable, PointingResidualTable @@ -12,7 +12,6 @@ from mapcat.pointing.poly import PolynomialPointingModel - def test_add_retrieve_pointing(database_sessionmaker): model = ConstantPointingModel(ra_offset=0.5 * u.deg, dec_offset=0.5 * u.deg) From ce8e58ea11662ee6e8544cd27df320f49493c972 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Wed, 29 Apr 2026 16:52:25 -0400 Subject: [PATCH 16/20] ruff again.. now updated to 15.12 --- ...b2c3d4e5f6_add_residual_stats_to_pointing_residuals.py | 8 ++++---- ...d9bc4ba5bc0_change_pointing_model_to_pydantic_model.py | 8 ++++---- mapcat/database/pointing_residual.py | 4 +--- 3 files changed, 9 insertions(+), 11 deletions(-) diff --git a/mapcat/alembic/versions/a1b2c3d4e5f6_add_residual_stats_to_pointing_residuals.py b/mapcat/alembic/versions/a1b2c3d4e5f6_add_residual_stats_to_pointing_residuals.py index 1aaa0a1..225cd64 100644 --- a/mapcat/alembic/versions/a1b2c3d4e5f6_add_residual_stats_to_pointing_residuals.py +++ b/mapcat/alembic/versions/a1b2c3d4e5f6_add_residual_stats_to_pointing_residuals.py @@ -6,16 +6,16 @@ """ -from typing import Sequence, Union +from collections.abc import Sequence import sqlalchemy as sa from alembic import op # revision identifiers, used by Alembic. revision: str = "a1b2c3d4e5f6" -down_revision: Union[str, None] = "cd9bc4ba5bc0" -branch_labels: Union[str, Sequence[str], None] = None -depends_on: Union[str, Sequence[str], None] = None +down_revision: str | None = "cd9bc4ba5bc0" +branch_labels: str | Sequence[str] | None = None +depends_on: str | Sequence[str] | None = None def upgrade() -> None: diff --git a/mapcat/alembic/versions/cd9bc4ba5bc0_change_pointing_model_to_pydantic_model.py b/mapcat/alembic/versions/cd9bc4ba5bc0_change_pointing_model_to_pydantic_model.py index 1a1ed8c..1723b40 100644 --- a/mapcat/alembic/versions/cd9bc4ba5bc0_change_pointing_model_to_pydantic_model.py +++ b/mapcat/alembic/versions/cd9bc4ba5bc0_change_pointing_model_to_pydantic_model.py @@ -6,16 +6,16 @@ """ -from typing import Sequence, Union +from collections.abc import Sequence import sqlalchemy as sa from alembic import op # revision identifiers, used by Alembic. revision: str = "cd9bc4ba5bc0" -down_revision: Union[str, None] = "6ce7e94dfd2d" -branch_labels: Union[str, Sequence[str], None] = None -depends_on: Union[str, Sequence[str], None] = None +down_revision: str | None = "6ce7e94dfd2d" +branch_labels: str | Sequence[str] | None = None +depends_on: str | Sequence[str] | None = None def upgrade() -> None: diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index 32d01cd..b58ab6b 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -2,8 +2,6 @@ Table containing pointing residuals. """ -from typing import Union - from sqlmodel import Field, Relationship, SQLModel from mapcat.pointing.base import PointingModelStats @@ -13,7 +11,7 @@ from .depth_one_map import DepthOneMapTable from .json import JSONEncodedPydantic -PointingModel = Union[ConstantPointingModel, PolynomialPointingModel] +PointingModel = ConstantPointingModel | PolynomialPointingModel class PointingResidualTable(SQLModel, table=True): """ From 19b6de39384b6013e951b8f03fba6db9882a3035 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Thu, 30 Apr 2026 08:55:33 -0400 Subject: [PATCH 17/20] astropydantic requirement --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 1163952..73a6fe8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,6 +10,7 @@ name = "mapcat" version = "0.2.2" requires-python = ">=3.10" dependencies = [ + "astropydantic", "sqlmodel", "sqlalchemy[asyncio]", "aiosqlite", From b5a76acc46eb97967602088a8033ca3860da5e60 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Thu, 30 Apr 2026 11:10:05 -0400 Subject: [PATCH 18/20] rather than have both ra,dec uncertainties and weights, just have weights since they were redundant --- mapcat/pointing/poly.py | 31 ++++--------------------------- 1 file changed, 4 insertions(+), 27 deletions(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index f6f9d2c..b4d8788 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -58,23 +58,15 @@ def _poly_keys(self): def calculate_model(self, measured_positions: SkyCoord, expected_positions: SkyCoord, - position_uncertainty: tuple[u.Quantity, u.Quantity] | u.Quantity | None = None, weights: tuple[list[float], list[float]] | list[float] | None = None ): """ Calculate and set the polynomial coefficients for the pointing model using the measured and expected positions. - - Optionally accept uncertainties in the measured positions - and measurement weights to perform a weighted fit. - - position_uncertainty can be provided as a tuple of (ra_uncertainties, dec_uncertainties) - or a single quantity that applies to both. - similarly, weights can be provided as a tuple of (ra_weights, dec_weights) + weights can be provided as a tuple of (ra_weights, dec_weights) or a single list that applies to both. - - weights are resolved from uncertainties if not provided, using inverse variance weighting. + Raises ------ @@ -92,12 +84,6 @@ def calculate_model(self, if n == 0: raise ValueError("No positions provided for model calculation.") - # Unpack position_uncertainty into ra_uncerts, dec_uncerts - if isinstance(position_uncertainty, tuple): - ra_uncerts, dec_uncerts = position_uncertainty - else: - ra_uncerts = dec_uncerts = position_uncertainty # None or single Quantity applied to both - # Unpack weights into ra_weights, dec_weights if isinstance(weights, tuple): ra_weights, dec_weights = weights @@ -107,17 +93,8 @@ def calculate_model(self, # Lots of logic to check if weights exist, etc. # Resolve weights from uncertainties if not provided if ra_weights is None and dec_weights is None: - if ra_uncerts is not None and dec_uncerts is not None: - ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts.to_value(u.deg)**2 + dec_uncerts.to_value(u.deg)**2) - elif ra_uncerts is not None: - dec_uncerts = np.ones(n) * np.mean(ra_uncerts.to_value(u.deg)) - ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts.to_value(u.deg)**2 + dec_uncerts**2) - elif dec_uncerts is not None: - ra_uncerts = np.ones(n) * np.mean(dec_uncerts.to_value(u.deg)) - ra_weights = dec_weights = 1 / np.sqrt(ra_uncerts**2 + dec_uncerts.to_value(u.deg)**2) - else: - # No uncertainty or weights provided — uniform - ra_weights = dec_weights = np.ones(n) + # No weights provided — uniform + ra_weights = dec_weights = np.ones(n) elif ra_weights is None: ra_weights = dec_weights elif dec_weights is None: From b2126e1ca05ffc5e826bd3062c8a0e5171209505 Mon Sep 17 00:00:00 2001 From: "Allen M. Foster" Date: Thu, 30 Apr 2026 11:41:29 -0400 Subject: [PATCH 19/20] build the model, rather than calculate themodel --- mapcat/pointing/poly.py | 2 +- tests/test_pointing.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index b4d8788..104aa4f 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -55,7 +55,7 @@ def _poly_keys(self): keys.append(f"x^{i}y^{j}") return keys - def calculate_model(self, + def build_model(self, measured_positions: SkyCoord, expected_positions: SkyCoord, weights: tuple[list[float], list[float]] | list[float] | None = None diff --git a/tests/test_pointing.py b/tests/test_pointing.py index a133dc0..c11e5b9 100644 --- a/tests/test_pointing.py +++ b/tests/test_pointing.py @@ -73,7 +73,7 @@ def test_make_polynomial_pointing_model(): offset_positions = SkyCoord( ra=ras + offset + slope * ras, dec=decs + offset + slope * decs ) - model.calculate_model(measured_positions=offset_positions, expected_positions=SkyCoord(ras, decs)) + model.build_model(measured_positions=offset_positions, expected_positions=SkyCoord(ras, decs)) assert model.ra_model_coefficients is not None assert model.dec_model_coefficients is not None From bafb7cf87aa49e153c57b2661d49efc3451f5ca7 Mon Sep 17 00:00:00 2001 From: Josh Borrow Date: Fri, 1 May 2026 14:12:01 -0400 Subject: [PATCH 20/20] Formatting --- mapcat/database/pointing_residual.py | 7 +-- mapcat/pointing/poly.py | 78 ++++++++++++++++------------ tests/test_pointing.py | 31 +++++++---- 3 files changed, 72 insertions(+), 44 deletions(-) diff --git a/mapcat/database/pointing_residual.py b/mapcat/database/pointing_residual.py index b58ab6b..7f9a14e 100644 --- a/mapcat/database/pointing_residual.py +++ b/mapcat/database/pointing_residual.py @@ -13,6 +13,7 @@ PointingModel = ConstantPointingModel | PolynomialPointingModel + class PointingResidualTable(SQLModel, table=True): """ Table for tracking Pointing error for a depth one map, @@ -38,8 +39,8 @@ class PointingResidualTable(SQLModel, table=True): foreign_key="depth_one_maps.map_id", ondelete="CASCADE", ) - residual_model: PointingModel = Field(sa_type=JSONEncodedPydantic(PointingModel) - ) - residual_stats: PointingModelStats | None = Field(nullable=True, default=None, sa_type=JSONEncodedPydantic(PointingModelStats) + residual_model: PointingModel = Field(sa_type=JSONEncodedPydantic(PointingModel)) + residual_stats: PointingModelStats | None = Field( + nullable=True, default=None, sa_type=JSONEncodedPydantic(PointingModelStats) ) map: DepthOneMapTable = Relationship(back_populates="pointing_residual") diff --git a/mapcat/pointing/poly.py b/mapcat/pointing/poly.py index 104aa4f..38d200a 100644 --- a/mapcat/pointing/poly.py +++ b/mapcat/pointing/poly.py @@ -26,13 +26,13 @@ class PolynomialCoefficients(BaseModel): } labels = {'x':'ra', 'y':'dec'} """ - coeffs: dict[str,float] + + coeffs: dict[str, float] labels: dict[str, str] unit: AstroPydanticUnit = u.deg poly_order: int - class PolynomialPointingModel(PointingModelProtocol): model_type: Literal["polynomial"] = "polynomial" @@ -55,22 +55,23 @@ def _poly_keys(self): keys.append(f"x^{i}y^{j}") return keys - def build_model(self, - measured_positions: SkyCoord, - expected_positions: SkyCoord, - weights: tuple[list[float], list[float]] | list[float] | None = None - ): + def build_model( + self, + measured_positions: SkyCoord, + expected_positions: SkyCoord, + weights: tuple[list[float], list[float]] | list[float] | None = None, + ): """ Calculate and set the polynomial coefficients for the pointing model using the measured and expected positions. - weights can be provided as a tuple of (ra_weights, dec_weights) + weights can be provided as a tuple of (ra_weights, dec_weights) or a single list that applies to both. - + Raises ------ - ValueError + ValueError If no positions are provided for model calculation. ValueError If the lengths of weights do not match the number of positions. @@ -83,7 +84,7 @@ def build_model(self, n = len(ra_offsets) if n == 0: raise ValueError("No positions provided for model calculation.") - + # Unpack weights into ra_weights, dec_weights if isinstance(weights, tuple): ra_weights, dec_weights = weights @@ -102,8 +103,12 @@ def build_model(self, ra_weights = np.asarray(ra_weights) dec_weights = np.asarray(dec_weights) - assert len(ra_weights) == n, "Length of ra_weights must match number of positions" - assert len(dec_weights) == n, "Length of dec_weights must match number of positions" + assert len(ra_weights) == n, ( + "Length of ra_weights must match number of positions" + ) + assert len(dec_weights) == n, ( + "Length of dec_weights must match number of positions" + ) ras = measured_positions.ra.to_value(u.deg) decs = measured_positions.dec.to_value(u.deg) @@ -125,25 +130,34 @@ def build_model(self, Aw = A_dec * w_dec[:, None] yw = y_dec * w_dec coeffs_dec, *_ = np.linalg.lstsq(Aw, yw, rcond=None) - - ra_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_ra)} - dec_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_dec)} - - self.ra_model_coefficients=PolynomialCoefficients(coeffs=ra_coeff_dict, labels={'x':'ra', 'y':'dec'}, unit=u.deg, poly_order=self.poly_order) - self.dec_model_coefficients=PolynomialCoefficients(coeffs=dec_coeff_dict, labels={'x':'ra', 'y':'dec'}, unit=u.deg, poly_order=self.poly_order) - - def model_fn( - self, x: u.Quantity, y: u.Quantity, coeffs: np.ndarray - ) -> u.Quantity: + ra_coeff_dict = {key: coeff for key, coeff in zip(self._poly_keys(), coeffs_ra)} + dec_coeff_dict = { + key: coeff for key, coeff in zip(self._poly_keys(), coeffs_dec) + } + + self.ra_model_coefficients = PolynomialCoefficients( + coeffs=ra_coeff_dict, + labels={"x": "ra", "y": "dec"}, + unit=u.deg, + poly_order=self.poly_order, + ) + self.dec_model_coefficients = PolynomialCoefficients( + coeffs=dec_coeff_dict, + labels={"x": "ra", "y": "dec"}, + unit=u.deg, + poly_order=self.poly_order, + ) + + def model_fn(self, x: u.Quantity, y: u.Quantity, coeffs: np.ndarray) -> u.Quantity: x = np.atleast_1d(x.to_value(u.deg)) y = np.atleast_1d(y.to_value(u.deg)) T = self._poly_terms(x, y) return (T @ coeffs) * u.deg - + def extract_coefficients(self) -> tuple[np.ndarray, np.ndarray]: """ - Extract the coefficients from the PolynomialCoefficients dataclasss and + Extract the coefficients from the PolynomialCoefficients dataclasss and return them as arrays in the correct order for the model function. Raises @@ -153,26 +167,26 @@ def extract_coefficients(self) -> tuple[np.ndarray, np.ndarray]: """ if self.ra_model_coefficients is None or self.dec_model_coefficients is None: raise ValueError("Model coefficients have not been calculated yet.") - + ra_coeff_array = np.zeros(len(self.ra_model_coefficients.coeffs)) for i, key in enumerate(self._poly_keys()): ra_coeff_array[i] = self.ra_model_coefficients.coeffs.get(key, 0) - + dec_coeff_array = np.zeros(len(self.dec_model_coefficients.coeffs)) for i, key in enumerate(self._poly_keys()): dec_coeff_array[i] = self.dec_model_coefficients.coeffs.get(key, 0) - + return ra_coeff_array, dec_coeff_array - + def predict(self, pos: SkyCoord) -> SkyCoord: - racoeffs,deccoeffs = self.extract_coefficients() + racoeffs, deccoeffs = self.extract_coefficients() ra_offset = self.model_fn(pos.ra, pos.dec, racoeffs) dec_offset = self.model_fn(pos.ra, pos.dec, deccoeffs) ra = pos.ra - ra_offset dec = pos.dec - dec_offset return SkyCoord(ra=ra, dec=dec, frame=pos.frame) - + def calculate_statistics(self, positions: SkyCoord): new_positions = self.predict(positions) ra_residuals = (new_positions.ra - positions.ra).to(u.arcsec) @@ -186,4 +200,4 @@ def calculate_statistics(self, positions: SkyCoord): mean_dec_offset=mean_dec, stddev_ra_offset=std_ra, stddev_dec_offset=std_dec, - ) \ No newline at end of file + ) diff --git a/tests/test_pointing.py b/tests/test_pointing.py index c11e5b9..44fdff7 100644 --- a/tests/test_pointing.py +++ b/tests/test_pointing.py @@ -55,8 +55,7 @@ def test_make_constant_pointing_model(): assert model.ra_offset == 0.5 * u.deg assert model.dec_offset == 0.5 * u.deg - og_pos = SkyCoord(ra=10 * u.deg, dec=20 * u.deg - ) + og_pos = SkyCoord(ra=10 * u.deg, dec=20 * u.deg) offset_pos = SkyCoord(og_pos.ra + model.ra_offset, og_pos.dec + model.dec_offset) new_pos = model.predict(offset_pos) assert new_pos.ra == og_pos.ra @@ -71,18 +70,32 @@ def test_make_polynomial_pointing_model(): slope = 0.1 * u.arcmin / u.deg offset_positions = SkyCoord( - ra=ras + offset + slope * ras, dec=decs + offset + slope * decs - ) - model.build_model(measured_positions=offset_positions, expected_positions=SkyCoord(ras, decs)) + ra=ras + offset + slope * ras, dec=decs + offset + slope * decs + ) + model.build_model( + measured_positions=offset_positions, expected_positions=SkyCoord(ras, decs) + ) assert model.ra_model_coefficients is not None assert model.dec_model_coefficients is not None for i, offset_pos in enumerate(offset_positions): predicted_pos = model.predict(offset_pos) - assert np.isclose(predicted_pos.ra.to_value(u.arcmin), ras[i].to_value(u.arcmin), atol=0.1) - assert np.isclose(predicted_pos.dec.to_value(u.arcmin), decs[i].to_value(u.arcmin), atol=0.1) + assert np.isclose( + predicted_pos.ra.to_value(u.arcmin), ras[i].to_value(u.arcmin), atol=0.1 + ) + assert np.isclose( + predicted_pos.dec.to_value(u.arcmin), decs[i].to_value(u.arcmin), atol=0.1 + ) predicted_pos = model.predict(offset_positions) - assert np.all(np.isclose(predicted_pos.ra.to_value(u.arcmin), ras.to_value(u.arcmin), atol=0.1)) - assert np.all(np.isclose(predicted_pos.dec.to_value(u.arcmin), decs.to_value(u.arcmin), atol=0.1)) + assert np.all( + np.isclose( + predicted_pos.ra.to_value(u.arcmin), ras.to_value(u.arcmin), atol=0.1 + ) + ) + assert np.all( + np.isclose( + predicted_pos.dec.to_value(u.arcmin), decs.to_value(u.arcmin), atol=0.1 + ) + )