From c967d9558410589275ccdfe098e985f1f5b8f794 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 18 Jul 2025 14:03:28 +0100 Subject: [PATCH 1/2] ellipse fixes --- autogalaxy/config/visualize/plots.yaml | 1 + .../ellipse/ellipse/ellipse_multipole.py | 4 +- autogalaxy/ellipse/fit_ellipse.py | 28 +- autogalaxy/ellipse/model/plotter_interface.py | 12 +- .../ellipse/plot/fit_ellipse_plot_util.py | 295 +++++------------- .../ellipse/plot/fit_ellipse_plotters.py | 150 ++++++++- autogalaxy/plot/__init__.py | 1 + autogalaxy/plot/mat_plot/two_d.py | 4 + autogalaxy/plot/visuals/two_d.py | 2 + autogalaxy/util/error_util.py | 83 +++++ 10 files changed, 341 insertions(+), 239 deletions(-) diff --git a/autogalaxy/config/visualize/plots.yaml b/autogalaxy/config/visualize/plots.yaml index 1bcb37df8..b26d6a05f 100644 --- a/autogalaxy/config/visualize/plots.yaml +++ b/autogalaxy/config/visualize/plots.yaml @@ -51,5 +51,6 @@ fit_interferometer: # Settings for plots of fits to inter fit_ellipse: # Settings for plots of ellipse fitting fits (e.g. FitEllipse) data : true # Plot the data of the ellipse fit? data_no_ellipse: true # Plot the data without the black data ellipses, which obscure noisy data? + ellipse_residuals: true # Plot the residuals of the ellipse fit? fit_quantity: {} # Settings for plots of fit quantities (e.g. FitQuantityPlotter). \ No newline at end of file diff --git a/autogalaxy/ellipse/ellipse/ellipse_multipole.py b/autogalaxy/ellipse/ellipse/ellipse_multipole.py index 5991039ed..89ac08517 100644 --- a/autogalaxy/ellipse/ellipse/ellipse_multipole.py +++ b/autogalaxy/ellipse/ellipse/ellipse_multipole.py @@ -38,7 +38,7 @@ class representing the multipole of an ellispe with, which is used to perform el self.multipole_comps = multipole_comps def points_perturbed_from( - self, pixel_scale, points, ellipse: Ellipse + self, pixel_scale, points, ellipse: Ellipse, n_i: int = 0 ) -> np.ndarray: """ Returns the (y,x) coordinates of the input points, which are perturbed by the multipole of the ellipse. @@ -57,7 +57,7 @@ def points_perturbed_from( The (y,x) coordinates of the input points, which are perturbed by the multipole. """ - angles = ellipse.angles_from_x0_from(pixel_scale=pixel_scale) + angles = ellipse.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i) radial = np.add( self.multipole_comps[1] * np.cos(self.m * (angles - ellipse.angle_radians)), diff --git a/autogalaxy/ellipse/fit_ellipse.py b/autogalaxy/ellipse/fit_ellipse.py index fa5b3d7d0..f0a2b70f2 100644 --- a/autogalaxy/ellipse/fit_ellipse.py +++ b/autogalaxy/ellipse/fit_ellipse.py @@ -56,6 +56,14 @@ def points_from_major_axis_from(self) -> np.ndarray: pixel_scale=self.dataset.pixel_scales[0], ) + if self.multipole_list is not None: + for multipole in self.multipole_list: + points = multipole.points_perturbed_from( + pixel_scale=self.dataset.pixel_scales[0], + points=points, + ellipse=self.ellipse, + ) + if self.interp.mask_interp is not None: i_total = 300 @@ -74,6 +82,16 @@ def points_from_major_axis_from(self) -> np.ndarray: pixel_scale=self.dataset.pixel_scales[0], n_i=i ) + if self.multipole_list is not None: + for multipole in self.multipole_list: + + points = multipole.points_perturbed_from( + pixel_scale=self.dataset.pixel_scales[0], + points=points, + ellipse=self.ellipse, + n_i=i + ) + if i == i_total: raise ValueError( @@ -88,14 +106,6 @@ def points_from_major_axis_from(self) -> np.ndarray: points = points[self.interp.mask_interp(points) == 0] - if self.multipole_list is not None: - for multipole in self.multipole_list: - points = multipole.points_perturbed_from( - pixel_scale=self.dataset.pixel_scales[0], - points=points, - ellipse=self.ellipse, - ) - return points @cached_property @@ -279,4 +289,4 @@ def figure_of_merit(self) -> float: ------- The figure of merit of the fit. """ - return self.log_likelihood + return -0.5* (self.chi_squared + self.noise_normalization) diff --git a/autogalaxy/ellipse/model/plotter_interface.py b/autogalaxy/ellipse/model/plotter_interface.py index 3e5f6a1a9..61412bd03 100644 --- a/autogalaxy/ellipse/model/plotter_interface.py +++ b/autogalaxy/ellipse/model/plotter_interface.py @@ -93,16 +93,16 @@ def should_plot(name): fit_list=fit_list, mat_plot_2d=mat_plot_2d, include_2d=self.include_2d ) - try: - fit_plotter.figures_2d(data=should_plot("data")) - except ValueError: - print(fit_plotter.fit_list[0].ellipse.major_axis) - print(fit_plotter.fit_list[0].ellipse.ell_comps) + fit_plotter.figures_2d( + data=should_plot("data"), + ellipse_residuals=should_plot("ellipse_residuals"), + ) if should_plot("data_no_ellipse"): fit_plotter.figures_2d( data=True, disable_data_contours=True, + ) fit_plotter.mat_plot_2d.use_log10 = True @@ -113,4 +113,4 @@ def should_plot(name): fit_plotter.figures_2d( data=True, disable_data_contours=True, - ) + ) \ No newline at end of file diff --git a/autogalaxy/ellipse/plot/fit_ellipse_plot_util.py b/autogalaxy/ellipse/plot/fit_ellipse_plot_util.py index 43bf81b5d..8ef3caf2b 100644 --- a/autogalaxy/ellipse/plot/fit_ellipse_plot_util.py +++ b/autogalaxy/ellipse/plot/fit_ellipse_plot_util.py @@ -1,236 +1,99 @@ -import copy, matplotlib -import numpy as np -import matplotlib.pyplot as plt - from astropy import units -from matplotlib.colors import LinearSegmentedColormap - - -def subplot_fit(fit_list): - # NOTE: - logscale = False - - # NOTE: - figure, axes = plt.subplots(nrows=1, ncols=2, figsize=(16, 8)) - - image = fit.data.native - image_temp = image - - # NOTE: - # image = copy.copy(sample.image) - # if sample.mask is not None: - # image[sample.mask.astype(bool)] = np.nan - # - # image_temp = copy.copy(sample.image) - # image_temp[~sample.mask.astype(bool)] = np.nan - # else: - # image_temp = None - - # NOTE: - def custom_colormap(): - # Define the colors - # colors = ['white', 'black', 'red'] - # positions = [0.0, 0.5, 1.0] - colors = ["grey", "white", "red", "darkred", "black"] - positions = [0.0, 0.25, 0.5, 0.75, 1.0] - - # Create the colormap - cmap = LinearSegmentedColormap.from_list( - "custom_colormap", list(zip(positions, colors)) - ) - - return cmap - - vmin = None - vmax = None - - if vmin is None: - vmin = 0.025 - if vmax is None: - vmax = 5.0 - - norm = matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax) - - im = axes[0].imshow( - np.log10(image) if logscale else image, - cmap="jet", - aspect="auto", - norm=norm, - ) - if image_temp is not None: - axes[0].imshow( - np.log10(image_temp) if logscale else image_temp, - cmap="jet", - aspect="auto", - alpha=0.5, - ) +import itertools +import matplotlib.pyplot as plt +import numpy as np - ellipse_list = [fit.ellipse for fit in fit_list] - list_of_angles = [] - list_of_y_fit = [] - list_of_y_errors_fit = [] - y_means = [] - y_stds = [] +def plot_ellipse_residuals( + array, + fit_list, + colors, + output, + for_subplot : bool = False +): - radii + color = itertools.cycle(colors) - for i, (a, parameters) in enumerate(zip(array, list_of_parameters)): - if i == 0: - m = sum(1 for value in parameters.values() if value is not None) + if for_subplot: + ax = plt.subplot(1, 2, 2) + else: - y_fit, y_errors_fit, (x, y), angles = sample.extract( - a=a, parameters=parameters, condition=extract_condition + fig = plt.figure( + figsize=(8, 8) ) - list_of_angles.append(angles) - list_of_y_fit.append(y_fit) - list_of_y_errors_fit.append(y_errors_fit) - if y_errors_fit is None: - # y_errors_fit = 0.05 * y_fit - raise NotImplementedError() + ax = fig.add_subplot(111) - # NOTE: - y_mean = np.nanmean(y_fit) - y_means.append(y_mean) + for i, fit in enumerate(fit_list): - # NOTE: - y_std = np.nanstd(y_fit) - y_stds.append(y_std) + angles = fit.ellipse.angles_from_x0_from(pixel_scale=array.pixel_scales[0]) - # NOTE: - axes[0].plot(x, y, marker="o", markersize=2.5, color="w") + color_plot = next(color) - # NOTE: - axes[1].errorbar( + plt.errorbar( angles * units.rad.to(units.deg), - y_fit, - yerr=y_errors_fit, + fit.data_interp, + yerr=fit.noise_map_interp, linestyle="None", marker="o", markersize=2.5, - color="black", + color=color_plot, + ) + plt.axhline( + np.nanmean(fit.data_interp), + linestyle=":", + color=color_plot, ) - axes[1].axhline(y_mean, linestyle=":", color="black") - - levels = np.sort(np.log10(y_means)) if logscale else np.sort(y_means) - axes[0].contour( - np.log10(image) if logscale else image, - # levels=y_means[::-1], - levels=levels, - colors="black", + + ax.set_xticks([0, 90, 180, 270, 360]) + + ax.set_xlabel( + r"$\phi$ (deg)", + fontsize=20, + ) + ax.set_ylabel( + r"$\rm I(\phi)$ [E/s]", + fontsize=20, + labelpad=-5., + ) + ax.tick_params(axis='y', labelsize=12.5) + ax.yaxis.tick_right() + ax.yaxis.set_label_position("right") + + ax.minorticks_on() + ax.tick_params( + axis='both', + which="major", + length=6, + width=2, + right=True, + top=True, + direction="in", + colors="black", labelsize=15 + ) + ax.tick_params( + axis='both', + which="minor", + length=3, + width=1, + right=True, + top=True, + direction="in", + colors="black", labelsize=15 ) - colors = [im.cmap(im.norm(level)) for level in levels][::-1] - for i, (angles, y_fit, y_errors_fit) in enumerate( - zip(list_of_angles, list_of_y_fit, list_of_y_errors_fit) - ): - axes[1].errorbar( - angles * units.rad.to(units.deg), - y_fit, - yerr=y_errors_fit, - linestyle="None", - marker="o", - markersize=2.5, - color=colors[i], - ) - # axes[0].plot( - # [247], - # [250], - # marker="o" - # ) - xticks = np.linspace(0, image.shape[1], 11) - yticks = np.linspace(0, image.shape[0], 11) - axes[0].set_xticks(xticks) - axes[0].set_yticks(xticks) - axes[1].set_xticks([0, 90, 180, 270, 360]) - axes[1].set_xlabel(r"$\phi$ (deg)", fontsize=15) - axes[1].set_ylabel(r"$\rm I(\phi)$ [E/s]", fontsize=15) - axes[1].tick_params(axis="y", labelsize=12.5) - axes[1].yaxis.tick_right() - axes[1].yaxis.set_label_position("right") - for i, ax in enumerate(axes): - ax.minorticks_on() - ax.tick_params( - axis="both", - which="major", - length=6, - right=True, - top=True, - direction="in", - colors="w" if i == 0 else "black", - ) - ax.tick_params( - axis="both", - which="minor", - length=3, - right=True, - top=True, - direction="in", - colors="w" if i == 0 else "black", - ) + ax.set_yscale("log") + + ax.set_title("Ellipse 1D Residuals", fontsize=20) + + if not for_subplot: + + bbox_original = output.bbox_inches + output.bbox_inches = None + + output.to_figure(structure=None, auto_filename="ellipse_residuals") + + output.bbox_inches = bbox_original + + plt.close() - axes[1].set_yscale("log") - - # text = axes[0].text( - # 0.05, - # 0.95, - # "model 1", - # horizontalalignment='left', - # verticalalignment='center', - # transform=axes[0].transAxes, - # fontsize=25, - # weight="bold", - # color="w" - # ) - # text.set_path_effects([patheffects.withStroke(linewidth=3, foreground='black')]) - - # text = axes[0].text( - # 0.7, - # 0.95, - # "NGC 2274", - # horizontalalignment='left', - # verticalalignment='center', - # transform=axes[0].transAxes, - # fontsize=25, - # weight="bold", - # color="w" - # ) - # text.set_path_effects([patheffects.withStroke(linewidth=3, foreground='black')]) - - if a_max is None: - axes[0].set_xlim(0, image.shape[1]) - axes[0].set_ylim(0, image.shape[0]) - else: - x0 = list_of_parameters[0]["x0"] - y0 = list_of_parameters[0]["y0"] - axes[0].set_xlim(x0 - a_max, x0 + a_max) - axes[0].set_ylim(y0 - a_max, y0 + a_max) - - # NOTE: - figure.subplots_adjust(left=0.05, right=0.95, bottom=0.075, top=0.95, wspace=0.0) - - y = [] - for i, (a, chi_squares_i) in enumerate(zip(array, chi_squares)): - n = len(chi_squares_i) - 1 - - # NOTE: - y_i = np.nansum(chi_squares_i[:-1]) / (n - m) - y.append(y_i) - - figure_stats, axes_stats = plt.subplots() - axes_stats.plot(array, y, linestyle="None", marker="o", color="black") - axes_stats.set_xscale("log") - axes_stats.set_yscale("log") - # directory = "./MASSIVE/metadata" - # filename = "{}/xy_model_default.numpy".format(directory) - # with open(filename, 'wb') as f: - # np.save(f, [x, y]) - # plt.show() - # exit() - - # NOTE: - # chi_squares_flattened = list(itertools.chain(*chi_squares)) - # plt.hist(chi_squares_flattened, bins=100, alpha=0.75) - - # plt.show();exit() diff --git a/autogalaxy/ellipse/plot/fit_ellipse_plotters.py b/autogalaxy/ellipse/plot/fit_ellipse_plotters.py index b6d99eeb7..9d0030c75 100644 --- a/autogalaxy/ellipse/plot/fit_ellipse_plotters.py +++ b/autogalaxy/ellipse/plot/fit_ellipse_plotters.py @@ -1,13 +1,11 @@ -import matplotlib.pyplot as plt import numpy as np -from typing import List - -from autoconf import conf +import math +from typing import List, Optional import autoarray as aa from autoarray import plot as aplt -from autoarray.plot.auto_labels import AutoLabels +from autogalaxy.ellipse.plot import fit_ellipse_plot_util from autogalaxy.ellipse.fit_ellipse import FitEllipse from autogalaxy.plot.abstract_plotters import Plotter @@ -18,6 +16,7 @@ from autogalaxy.plot.visuals.two_d import Visuals2D from autogalaxy.plot.include.two_d import Include2D +from autogalaxy.util import error_util class FitEllipsePlotter(Plotter): def __init__( @@ -51,6 +50,8 @@ def figures_2d( self, data: bool = False, disable_data_contours: bool = False, + ellipse_residuals: bool = False, + for_subplot : bool = False, suffix: str = "", ): """ @@ -71,6 +72,7 @@ def figures_2d( filename_tag = "" if data: + self.mat_plot_2d.contour = aplt.Contour( manual_levels=np.sort( [float(np.mean(fit.data_interp)) for fit in self.fit_list] @@ -93,7 +95,7 @@ def figures_2d( ellipse_list.append(aa.Grid2DIrregular.from_yx_1d(y=y, x=x)) visuals_2d = self.get_visuals_2d() + Visuals2D( - positions=ellipse_list, lines=ellipse_list + lines=ellipse_list ) self.mat_plot_2d.plot_array( @@ -107,3 +109,139 @@ def figures_2d( if disable_data_contours: self.mat_plot_2d.contour = contour_original + + if ellipse_residuals: + + try: + colors = self.mat_plot_2d.grid_plot.config_dict["c"] + except KeyError: + colors = "k" + + fit_ellipse_plot_util.plot_ellipse_residuals( + array=self.fit_list[0].dataset.data.native, + fit_list=self.fit_list, + colors=colors, + output=self.mat_plot_2d.output, + for_subplot=for_subplot, + ) + + def subplot_fit_ellipse(self): + + """ + Standard subplot of the attributes of the plotter's `FitEllipse` object. + """ + + self.open_subplot_figure(number_subplots=2) + + self.figures_2d(data=True) + self.figures_2d(ellipse_residuals=True, for_subplot=True) + + self.mat_plot_2d.output.subplot_to_figure(auto_filename="subplot_fit_ellipse") + self.close_subplot_figure() + +class FitEllipsePDFPlotter(Plotter): + def __init__( + self, + fit_pdf_list: List[FitEllipse], + mat_plot_1d: MatPlot1D = MatPlot1D(), + visuals_1d: Visuals1D = Visuals1D(), + include_1d: Include1D = Include1D(), + mat_plot_2d: MatPlot2D = MatPlot2D(), + visuals_2d: Visuals2D = Visuals2D(), + include_2d: Include2D = Include2D(), + sigma: Optional[float] = 3.0, + ): + super().__init__( + mat_plot_1d=mat_plot_1d, + visuals_1d=visuals_1d, + include_1d=include_1d, + mat_plot_2d=mat_plot_2d, + visuals_2d=visuals_2d, + include_2d=include_2d, + ) + + self.fit_pdf_list = fit_pdf_list + self.sigma = sigma + self.low_limit = (1 - math.erf(sigma / math.sqrt(2))) / 2 + + def get_visuals_1d(self) -> Visuals1D: + return self.visuals_1d + + def get_visuals_2d(self): + return self.get_2d.via_mask_from(mask=self.fit_pdf_list[0][0].dataset.mask) + + def subplot_ellipse_errors(self): + """ + Plots the individual attributes of the plotter's `FitEllipse` object in 1D. + + The API is such that every plottable attribute of the `FitEllipse` object is an input parameter of type bool of + the function, which if switched to `True` means that it is plotted. + + Parameters + ---------- + data + Whether to make a 1D plot (via `imshow`) of the image data. + disable_data_contours + If `True`, the data is plotted without the black data contours over the top (but the white contours + showing the ellipses are still plotted). + """ + + contour_original = self.mat_plot_2d.contour + self.mat_plot_2d.contour = False + + ellipse_centre_list = [] + fit_ellipse_list = [[] for _ in range(len(self.fit_pdf_list[0]))] + + for fit_list in self.fit_pdf_list: + + ellipse_centre_list.append(fit_list[0].ellipse.centre) + + for i, fit in enumerate(fit_list): + + points = fit.points_from_major_axis_from() + + x = points[:, 1] + y = points[:, 0] * -1.0 # flip for plot + + fit_ellipse_list[i].append( + aa.Grid2DIrregular.from_yx_1d(y=y, x=x) + ) + + print(i, len(x)) + + print() + + self.open_subplot_figure(number_subplots=len(fit_ellipse_list)) + + for i in range(len(fit_ellipse_list)): + + median_ellipse, [lower_ellipse, upper_ellipse] = error_util.ellipse_median_and_error_region_in_polar( + fit_ellipse_list[i], low_limit=self.low_limit, center=ellipse_centre_list[i] + ) + + # Unpack points + y_lower, x_lower = lower_ellipse[:, 0], lower_ellipse[:, 1] + y_upper, x_upper = upper_ellipse[:, 0], upper_ellipse[:, 1] + + # Close the contours + x_fill = np.concatenate([x_lower, x_upper[::-1]]) + y_fill = np.concatenate([y_lower, y_upper[::-1]]) + + visuals_2d = self.get_visuals_2d() + Visuals2D( + lines=median_ellipse, + fill_region=[y_fill, x_fill] + ) + + self.mat_plot_2d.plot_array( + array=self.fit_pdf_list[0][0].data, + visuals_2d=visuals_2d, + auto_labels=aplt.AutoLabels( + title=f"Ellipse Fit", + filename=f"subhplot_ellipse_errors", + ), + ) + + self.mat_plot_2d.output.subplot_to_figure(auto_filename="subplot_ellipse_errors") + self.close_subplot_figure() + + self.mat_plot_2d.contour = contour_original \ No newline at end of file diff --git a/autogalaxy/plot/__init__.py b/autogalaxy/plot/__init__.py index 81623d62f..24bc12594 100644 --- a/autogalaxy/plot/__init__.py +++ b/autogalaxy/plot/__init__.py @@ -93,3 +93,4 @@ from autogalaxy.galaxy.plot.galaxies_plotters import GalaxiesPlotter from autogalaxy.galaxy.plot.adapt_plotters import AdaptPlotter from autogalaxy.ellipse.plot.fit_ellipse_plotters import FitEllipsePlotter +from autogalaxy.ellipse.plot.fit_ellipse_plotters import FitEllipsePDFPlotter \ No newline at end of file diff --git a/autogalaxy/plot/mat_plot/two_d.py b/autogalaxy/plot/mat_plot/two_d.py index eebc44048..82f6de8b9 100644 --- a/autogalaxy/plot/mat_plot/two_d.py +++ b/autogalaxy/plot/mat_plot/two_d.py @@ -25,6 +25,7 @@ def __init__( output: Optional[aplt.Output] = None, array_overlay: Optional[aplt.ArrayOverlay] = None, contour: Optional[aplt.Contour] = None, + fill: Optional[aplt.Fill] = None, grid_scatter: Optional[aplt.GridScatter] = None, grid_plot: Optional[aplt.GridPlot] = None, vector_yx_quiver: Optional[aplt.VectorYXQuiver] = None, @@ -104,6 +105,8 @@ def __init__( Sets any annotations on the figure and customizes its appearance using `plt.annotate`. legend Sets whether the plot inclues a legend and customizes its appearance and labels using `plt.legend`. + fill + Sets the fill of the figure using `plt.fill` and customizes its appearance, such as the color and alpha. output Sets if the figure is displayed on the user's screen or output to `.png` using `plt.show` and `plt.savefig` array_overlay @@ -191,6 +194,7 @@ def __init__( xlabel=xlabel, text=text, annotate=annotate, + fill=fill, output=output, origin_scatter=origin_scatter, mask_scatter=mask_scatter, diff --git a/autogalaxy/plot/visuals/two_d.py b/autogalaxy/plot/visuals/two_d.py index 49f422979..3d9e425cf 100644 --- a/autogalaxy/plot/visuals/two_d.py +++ b/autogalaxy/plot/visuals/two_d.py @@ -17,6 +17,7 @@ def __init__( mesh_grid: aa.Grid2D = None, vectors: aa.VectorYX2DIrregular = None, patches: Union[ptch.Patch] = None, + fill_region: Optional[List] = None, array_overlay: aa.Array2D = None, light_profile_centres: aa.Grid2DIrregular = None, mass_profile_centres: aa.Grid2DIrregular = None, @@ -47,6 +48,7 @@ def __init__( mesh_grid=mesh_grid, vectors=vectors, patches=patches, + fill_region=fill_region, array_overlay=array_overlay, origin=origin, border=border, diff --git a/autogalaxy/util/error_util.py b/autogalaxy/util/error_util.py index 38650aedd..7c1810a2e 100644 --- a/autogalaxy/util/error_util.py +++ b/autogalaxy/util/error_util.py @@ -66,3 +66,86 @@ def quantile_profile_1d(profile_1d_list, q, weights=None): radial_quantile[radial_index] = quantile(x=radial_list, q=q, weights=weights)[0] return radial_quantile + + +def quantile_ellipse(ellipse_list, q): + """ + Compute the per-point quantile ellipse for a list of 2D ellipses. + """ + stacked = np.stack(ellipse_list, axis=0) # shape: (n_samples, n_points, 2) + return np.quantile(stacked, q=q, axis=0) # shape: (n_points, 2) + +def ellipse_median_and_error_region_via_quantile(ellipse_list, low_limit): + """ + Compute the median and confidence bounds for an ellipse shape. + """ + median_ellipse = quantile_ellipse(ellipse_list, q=0.5) + lower_ellipse = quantile_ellipse(ellipse_list, q=low_limit) + upper_ellipse = quantile_ellipse(ellipse_list, q=1.0 - low_limit) + + return median_ellipse, [lower_ellipse, upper_ellipse] + + +def cartesian_to_polar(yx, center=(0.0, 0.0)): + y, x = yx[..., 0], yx[..., 1] + y0, x0 = center + r = np.sqrt((x - x0) ** 2 + (y - y0) ** 2) + theta = np.arctan2(y - y0, x - x0) + # Wrap to [0, 2π] + theta = np.mod(theta, 2 * np.pi) + return r, theta + + +def polar_to_cartesian(r, theta, center=(0.0, 0.0)): + y0, x0 = center + x = r * np.cos(theta) + x0 + y = r * np.sin(theta) + y0 + return np.stack([y, x], axis=-1) + + +def ellipse_median_and_error_region_in_polar( + ellipse_list, low_limit, center=(0.0, 0.0) +): + """ + Computes median and error ellipses in (y,x) space using polar quantiles. + + Parameters + ---------- + ellipse_list : list of arrays, each of shape (n_points, 2) + A list of (y, x) coordinate ellipses. + sigma : float + The sigma confidence level (e.g. 1.0 for 68%). + center : tuple of float + The center (y0, x0) for the polar projection. + + Returns + ------- + (median, [lower, upper]) : tuple of np.ndarray + Each array has shape (n_points, 2), giving the (y, x) coordinates + of the median and error ellipses. + """ + + # Convert all ellipses to polar (r, θ), shape: (n_ellipses, n_points) + r_list = [] + theta_ref = None + + for ellipse in ellipse_list: + + r, theta = cartesian_to_polar(ellipse, center=center) + r_list.append(r) + + if theta_ref is None: + theta_ref = theta # Assume all ellipses use same angle ordering + + r_array = np.stack(r_list, axis=0) # shape: (n_ellipses, n_points) + + median_r = np.quantile(r_array, q=0.5, axis=0) + lower_r = np.quantile(r_array, q=low_limit, axis=0) + upper_r = np.quantile(r_array, q=1 - low_limit, axis=0) + + # Convert back to (y, x) + median_yx = polar_to_cartesian(median_r, theta_ref, center=center) + lower_yx = polar_to_cartesian(lower_r, theta_ref, center=center) + upper_yx = polar_to_cartesian(upper_r, theta_ref, center=center) + + return median_yx, [lower_yx, upper_yx] \ No newline at end of file From 395033fede6dbcabef5ced4ba71dc10ba08004e2 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 18 Jul 2025 14:06:36 +0100 Subject: [PATCH 2/2] ellipse fixes --- autogalaxy/config/visualize/plots.yaml | 1 + autogalaxy/ellipse/model/plotter_interface.py | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/autogalaxy/config/visualize/plots.yaml b/autogalaxy/config/visualize/plots.yaml index b26d6a05f..3d82db942 100644 --- a/autogalaxy/config/visualize/plots.yaml +++ b/autogalaxy/config/visualize/plots.yaml @@ -49,6 +49,7 @@ fit_interferometer: # Settings for plots of fits to inter fits_dirty_images: false # output dirty_images.fits showing the dirty image, noise-map, model-data, resiual-map, normalized residual map and chi-squared map? fit_ellipse: # Settings for plots of ellipse fitting fits (e.g. FitEllipse) + subplot_fit_ellipse : true # Plot subplot of all fit quantities for ellipse fits (e.g. the model data, residual-map, etc.)? data : true # Plot the data of the ellipse fit? data_no_ellipse: true # Plot the data without the black data ellipses, which obscure noisy data? ellipse_residuals: true # Plot the residuals of the ellipse fit? diff --git a/autogalaxy/ellipse/model/plotter_interface.py b/autogalaxy/ellipse/model/plotter_interface.py index 61412bd03..a17e18b80 100644 --- a/autogalaxy/ellipse/model/plotter_interface.py +++ b/autogalaxy/ellipse/model/plotter_interface.py @@ -105,6 +105,10 @@ def should_plot(name): ) + if should_plot("subplot_fit_ellipse"): + + fit_plotter.subplot_fit_ellipse() + fit_plotter.mat_plot_2d.use_log10 = True fit_plotter.figures_2d(data=should_plot("data"))