diff --git a/.github/workflows/cd.yaml b/.github/workflows/cd.yaml index 5d9cd60b..c951fd76 100644 --- a/.github/workflows/cd.yaml +++ b/.github/workflows/cd.yaml @@ -58,7 +58,7 @@ jobs: ls -ltrh ls -ltrh dist - name: Publish to Test PyPI - uses: pypa/gh-action-pypi-publish@v1.9.0 + uses: pypa/gh-action-pypi-publish@v1.12.2 with: repository-url: https://test.pypi.org/legacy/ verbose: true @@ -96,5 +96,5 @@ jobs: name: artifact path: dist - - uses: pypa/gh-action-pypi-publish@v1.9.0 + - uses: pypa/gh-action-pypi-publish@v1.12.2 if: startsWith(github.ref, 'refs/tags') diff --git a/.github/workflows/coverage.yaml b/.github/workflows/coverage.yaml index 1a39bdd9..22e0d18d 100644 --- a/.github/workflows/coverage.yaml +++ b/.github/workflows/coverage.yaml @@ -52,9 +52,10 @@ jobs: pytest -vvv --cov=${{ env.PROJECT_NAME }} --cov-report=xml --cov-report=term tests/ - name: Upload coverage reports to Codecov with GitHub Action - uses: codecov/codecov-action@v4 + uses: codecov/codecov-action@v5 + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} with: - token: ${{ secrets.CODECOV_TOKEN }} files: ./coverage.xml flags: unittests name: codecov-umbrella diff --git a/astrophot/fit/lm.py b/astrophot/fit/lm.py index b00ebca8..44d40460 100644 --- a/astrophot/fit/lm.py +++ b/astrophot/fit/lm.py @@ -307,7 +307,6 @@ def step(self, chi2) -> torch.Tensor: # Keep track of chi^2 improvement even if it fails curvature test if chi2 <= scarry_best[1]: scarry_best = (ha, chi2, self.L) - nostep = False # Check for high curvature, in which case linear approximation is not valid. avoid this step rho = torch.linalg.norm(a) / torch.linalg.norm(h) diff --git a/astrophot/models/__init__.py b/astrophot/models/__init__.py index 7edea0eb..81edb2c8 100644 --- a/astrophot/models/__init__.py +++ b/astrophot/models/__init__.py @@ -8,6 +8,7 @@ from .flatsky_model import * from .planesky_model import * from .gaussian_model import * +from .multi_gaussian_expansion_model import * from .spline_model import * from .relspline_model import * from .psf_model_object import * diff --git a/astrophot/models/_shared_methods.py b/astrophot/models/_shared_methods.py index ca442bce..c00b3b26 100644 --- a/astrophot/models/_shared_methods.py +++ b/astrophot/models/_shared_methods.py @@ -14,6 +14,9 @@ moffat_torch, nuker_torch, ) +from ..utils.conversions.coordinates import ( + Rotate_Cartesian, +) from ..utils.decorators import ignore_numpy_warnings, default_internal from ..image import ( Image_List, @@ -280,6 +283,30 @@ def radial_evaluate_model(self, X=None, Y=None, image=None, parameters=None): ) +@default_internal +def transformed_evaluate_model(self, X=None, Y=None, image=None, parameters=None, **kwargs): + if X is None or Y is None: + Coords = image.get_coordinate_meshgrid() + X, Y = Coords - parameters["center"].value[..., None, None] + X, Y = self.transform_coordinates(X, Y, image, parameters) + return self.radial_model( + self.radius_metric(X, Y, image=image, parameters=parameters), + image=image, + parameters=parameters, + ) + + +# Transform Coordinates +###################################################################### +@default_internal +def inclined_transform_coordinates(self, X, Y, image=None, parameters=None): + X, Y = Rotate_Cartesian(-(parameters["PA"].value - image.north), X, Y) + return ( + X, + Y / parameters["q"].value, + ) + + # Exponential ###################################################################### @default_internal diff --git a/astrophot/models/galaxy_model_object.py b/astrophot/models/galaxy_model_object.py index 4e422b5e..7bad13b8 100644 --- a/astrophot/models/galaxy_model_object.py +++ b/astrophot/models/galaxy_model_object.py @@ -112,24 +112,5 @@ def initialize(self, target=None, parameters: Optional[Parameter_Node] = None, * if parameters["q"].uncertainty is None: parameters["q"].uncertainty = parameters["q"].value * self.default_uncertainty - @default_internal - def transform_coordinates(self, X, Y, image=None, parameters=None): - X, Y = Rotate_Cartesian(-(parameters["PA"].value - image.north), X, Y) - return ( - X, - Y / parameters["q"].value, - ) - - @default_internal - def evaluate_model( - self, X=None, Y=None, image=None, parameters: Parameter_Node = None, **kwargs - ): - if X is None or Y is None: - Coords = image.get_coordinate_meshgrid() - X, Y = Coords - parameters["center"].value[..., None, None] - XX, YY = self.transform_coordinates(X, Y, image, parameters) - return self.radial_model( - self.radius_metric(XX, YY, image, parameters), - image=image, - parameters=parameters, - ) + from ._shared_methods import inclined_transform_coordinates as transform_coordinates + from ._shared_methods import transformed_evaluate_model as evaluate_model diff --git a/astrophot/models/moffat_model.py b/astrophot/models/moffat_model.py index 8f6701a5..06961c8c 100644 --- a/astrophot/models/moffat_model.py +++ b/astrophot/models/moffat_model.py @@ -1,4 +1,5 @@ import torch +import numpy as np from .galaxy_model_object import Galaxy_Model from .psf_model_object import PSF_Model @@ -6,6 +7,7 @@ from ..utils.decorators import ignore_numpy_warnings, default_internal from ..utils.parametric_profiles import moffat_np from ..utils.conversions.functions import moffat_I0_to_flux, general_uncertainty_prop +from ..param import Param_Unlock, Param_SoftLimits __all__ = ["Moffat_Galaxy", "Moffat_PSF"] @@ -155,3 +157,35 @@ def total_flux_uncertainty(self, parameters=None): ) from ._shared_methods import radial_evaluate_model as evaluate_model + + +class Moffat2D_PSF(Moffat_PSF): + + model_type = f"moffat2d {PSF_Model.model_type}" + parameter_specs = { + "q": {"units": "b/a", "limits": (0, 1), "uncertainty": 0.03}, + "PA": { + "units": "radians", + "limits": (0, np.pi), + "cyclic": True, + "uncertainty": 0.06, + }, + } + _parameter_order = Moffat_PSF._parameter_order + ("q", "PA") + usable = True + model_integrated = False + + @select_target + @default_internal + def initialize(self, target=None, parameters=None, **kwargs): + with Param_Unlock(parameters["q"]), Param_SoftLimits(parameters["q"]): + if parameters["q"].value is None: + parameters["q"].value = 0.9 + + with Param_Unlock(parameters["PA"]), Param_SoftLimits(parameters["PA"]): + if parameters["PA"].value is None: + parameters["PA"].value = 0.1 + super().initialize(target=target, parameters=parameters) + + from ._shared_methods import inclined_transform_coordinates as transform_coordinates + from ._shared_methods import transformed_evaluate_model as evaluate_model diff --git a/astrophot/models/multi_gaussian_expansion_model.py b/astrophot/models/multi_gaussian_expansion_model.py new file mode 100644 index 00000000..dd71726b --- /dev/null +++ b/astrophot/models/multi_gaussian_expansion_model.py @@ -0,0 +1,171 @@ +import torch +import numpy as np +from scipy.stats import iqr + +from .psf_model_object import PSF_Model +from .model_object import Component_Model +from ._shared_methods import ( + select_target, +) +from ..utils.initialize import isophotes +from ..utils.angle_operations import Angle_COM_PA +from ..utils.conversions.coordinates import ( + Rotate_Cartesian, +) +from ..param import Param_Unlock, Param_SoftLimits, Parameter_Node +from ..utils.decorators import ignore_numpy_warnings, default_internal + +__all__ = ["Multi_Gaussian_Expansion"] + + +class Multi_Gaussian_Expansion(Component_Model): + """Model that represents a galaxy as a sum of multiple Gaussian + profiles. The model is defined as: + + I(R) = sum_i flux_i * exp(-0.5*(R_i / sigma_i)^2) / (2 * pi * q_i * sigma_i^2) + + where $R_i$ is a radius computed using $q_i$ and $PA_i$ for that component. All components share the same center. + + Parameters: + q: axis ratio to scale minor axis from the ratio of the minor/major axis b/a, this parameter is unitless, it is restricted to the range (0,1) + PA: position angle of the semi-major axis relative to the image positive x-axis in radians, it is a cyclic parameter in the range [0,pi) + sigma: standard deviation of each Gaussian + flux: amplitude of each Gaussian + """ + + model_type = f"mge {Component_Model.model_type}" + parameter_specs = { + "q": {"units": "b/a", "limits": (0, 1)}, + "PA": {"units": "radians", "limits": (0, np.pi), "cyclic": True}, + "sigma": {"units": "arcsec", "limits": (0, None)}, + "flux": {"units": "log10(flux)"}, + } + _parameter_order = Component_Model._parameter_order + ("q", "PA", "sigma", "flux") + usable = True + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + # determine the number of components + for key in ("q", "sigma", "flux"): + if self[key].value is not None: + self.n_components = self[key].value.shape[0] + break + else: + self.n_components = kwargs.get("n_components", 3) + + @torch.no_grad() + @ignore_numpy_warnings + @select_target + @default_internal + def initialize(self, target=None, parameters=None, **kwargs): + super().initialize(target=target, parameters=parameters) + + target_area = target[self.window] + target_dat = target_area.data.detach().cpu().numpy() + if target_area.has_mask: + mask = target_area.mask.detach().cpu().numpy() + target_dat[mask] = np.median(target_dat[np.logical_not(mask)]) + if parameters["sigma"].value is None: + with Param_Unlock(parameters["sigma"]), Param_SoftLimits(parameters["sigma"]): + parameters["sigma"].value = np.logspace( + np.log10(target_area.pixel_length.item() * 3), + max(target_area.shape.detach().cpu().numpy()) * 0.7, + self.n_components, + ) + parameters["sigma"].uncertainty = ( + self.default_uncertainty * parameters["sigma"].value + ) + if parameters["flux"].value is None: + with Param_Unlock(parameters["flux"]), Param_SoftLimits(parameters["flux"]): + parameters["flux"].value = np.log10( + np.sum(target_dat[~mask]) / self.n_components + ) * np.ones(self.n_components) + parameters["flux"].uncertainty = 0.1 * parameters["flux"].value + + if not (parameters["PA"].value is None or parameters["q"].value is None): + return + edge = np.concatenate( + ( + target_dat[:, 0], + target_dat[:, -1], + target_dat[0, :], + target_dat[-1, :], + ) + ) + edge_average = np.nanmedian(edge) + edge_scatter = iqr(edge[np.isfinite(edge)], rng=(16, 84)) / 2 + icenter = target_area.plane_to_pixel(parameters["center"].value) + + if parameters["PA"].value is None: + weights = target_dat - edge_average + Coords = target_area.get_coordinate_meshgrid() + X, Y = Coords - parameters["center"].value[..., None, None] + X, Y = X.detach().cpu().numpy(), Y.detach().cpu().numpy() + if target_area.has_mask: + seg = np.logical_not(target_area.mask.detach().cpu().numpy()) + PA = Angle_COM_PA(weights[seg], X[seg], Y[seg]) + else: + PA = Angle_COM_PA(weights, X, Y) + + with Param_Unlock(parameters["PA"]), Param_SoftLimits(parameters["PA"]): + parameters["PA"].value = ((PA + target_area.north) % np.pi) * np.ones( + self.n_components + ) + if parameters["PA"].uncertainty is None: + parameters["PA"].uncertainty = (5 * np.pi / 180) * torch.ones_like( + parameters["PA"].value + ) # default uncertainty of 5 degrees is assumed + if parameters["q"].value is None: + q_samples = np.linspace(0.2, 0.9, 15) + try: + pa = parameters["PA"].value.item() + except: + pa = parameters["PA"].value[0].item() + iso_info = isophotes( + target_area.data.detach().cpu().numpy() - edge_average, + (icenter[1].detach().cpu().item(), icenter[0].detach().cpu().item()), + threshold=3 * edge_scatter, + pa=(pa - target.north), + q=q_samples, + ) + with Param_Unlock(parameters["q"]), Param_SoftLimits(parameters["q"]): + parameters["q"].value = q_samples[ + np.argmin(list(iso["amplitude2"] for iso in iso_info)) + ] * torch.ones(self.n_components) + if parameters["q"].uncertainty is None: + parameters["q"].uncertainty = parameters["q"].value * self.default_uncertainty + + @default_internal + def total_flux(self, parameters=None): + return torch.sum(10 ** parameters["flux"].value) + + @default_internal + def evaluate_model(self, X=None, Y=None, image=None, parameters=None, **kwargs): + if X is None or Y is None: + Coords = image.get_coordinate_meshgrid() + X, Y = Coords - parameters["center"].value[..., None, None] + + if parameters["PA"].value.numel() == 1: + X, Y = Rotate_Cartesian(-(parameters["PA"].value - image.north), X, Y) + X = X.repeat(parameters["q"].value.shape[0], *[1] * X.ndim) + Y = torch.vmap(lambda q: Y / q)(parameters["q"].value) + else: + X, Y = torch.vmap(lambda pa: Rotate_Cartesian(-(pa - image.north), X, Y))( + parameters["PA"].value + ) + Y = torch.vmap(lambda q, y: y / q)(parameters["q"].value, Y) + + R = self.radius_metric(X, Y, image, parameters) + return torch.sum( + torch.vmap( + lambda A, R, sigma, q: (A / (2 * np.pi * q * sigma**2)) + * torch.exp(-0.5 * (R / sigma) ** 2) + )( + image.pixel_area * 10 ** parameters["flux"].value, + R, + parameters["sigma"].value, + parameters["q"].value, + ), + dim=0, + ) diff --git a/docs/source/_toc.yml b/docs/source/_toc.yml index 4c0f8233..ebd38f80 100644 --- a/docs/source/_toc.yml +++ b/docs/source/_toc.yml @@ -6,6 +6,7 @@ root: index chapters: - file: getting_started - file: install + - file: astrophotdemo - file: fastfit - file: tutorials/index - file: coordinates diff --git a/docs/source/astrophotdemo.rst b/docs/source/astrophotdemo.rst new file mode 100644 index 00000000..4f56e886 --- /dev/null +++ b/docs/source/astrophotdemo.rst @@ -0,0 +1,11 @@ +AstroPhot Demo +============== + +Go to the `AstroPhot Demo `_ to see a live demo of AstroPhot in action. + +In the demo you can upload your own FITS file and try fitting it with AstroPhot. +Add models and tweak the parameters until they seem reasonably good by eye, then +run the ``Optimize`` button to fit the model to the data. The demo is hosted on +Streamlit sharing and is free to use. Since the demo runs on the free version of +Streamlit, it may be slow to run all the functions (especially plotting), but it +will give you a sense of how AstroPhot works. diff --git a/docs/source/tutorials/ModelZoo.ipynb b/docs/source/tutorials/ModelZoo.ipynb index 876089f7..cc8a5307 100644 --- a/docs/source/tutorials/ModelZoo.ipynb +++ b/docs/source/tutorials/ModelZoo.ipynb @@ -231,6 +231,37 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2D Moffat PSF\n", + "\n", + "Like a Moffat, but can have a axis ratio and position angle. This could be used to make parametric spikes, or account for very slight asymmetry in a PSF." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "M = ap.models.AstroPhot_Model(\n", + " model_type=\"moffat2d psf model\",\n", + " parameters={\"n\": 2.0, \"Rd\": 10.0, \"q\": 0.7, \"PA\": 3.14 / 3},\n", + " target=psf_target,\n", + ")\n", + "print(M.parameter_order)\n", + "print(tuple(P.units for P in M.parameters))\n", + "M.initialize()\n", + "\n", + "fig, ax = plt.subplots(1, 2, figsize=(14, 6))\n", + "ap.plots.psf_image(fig, ax[0], M)\n", + "ap.plots.radial_light_profile(fig, ax[1], M)\n", + "ax[0].set_title(M.name)\n", + "plt.show()" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -624,6 +655,44 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Multi Gaussian Expansion\n", + "\n", + "A multi gaussian expansion is essentially a model made of overlapping gaussian models that share the same center. However, they are combined into a single model for computational efficiency. Another advantage of the MGE is that it is possible to determine a deprojection of the model from 2D into a 3D shape since the projection of a 3D gaussian is a 2D gaussian. Note however, that in some configurations this deprojection is not unique. See Cappellari 2002 for more details.\n", + "\n", + "Note: The ``PA`` can be either a single value (same for all components) or an array with values for each component." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "M = ap.models.AstroPhot_Model(\n", + " model_type=\"mge model\",\n", + " parameters={\n", + " \"center\": [50, 50],\n", + " \"q\": [0.9, 0.8, 0.6, 0.5],\n", + " \"PA\": 30 * np.pi / 180,\n", + " \"sigma\": [4.0, 8.0, 16.0, 32.0],\n", + " \"flux\": np.ones(4) / 4,\n", + " },\n", + " target=basic_target,\n", + ")\n", + "print(M.parameter_order)\n", + "print(tuple(P.units for P in M.parameters))\n", + "M.initialize()\n", + "\n", + "fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n", + "ap.plots.model_image(fig, ax, M)\n", + "ax.set_title(M.name)\n", + "plt.show()" + ] + }, { "cell_type": "markdown", "metadata": {},