AutoProf fails to fit isophotes of face-on galaxies

[hgao-astro]

The optimization of isophotal ellipticity fails when the initialized ellipticity is close to 0 or 1. Because the perturbation was done after converting the ellipticity to another parameter space of $[-\infty, +\infty]$ using _inv_x_to_eps and then scaling back using _x_to_eps. The functional form makes any finite perturbation at $\epsilon \rightarrow 0$ and $\epsilon \rightarrow 1$ (corresponds to $-\infty$ and $+\infty$ in the perturbed parameter space) negligible. Since there are few galaxies with global ellipticity close to 1, the problem mostly affects nearly face-on galaxies.

AutoProf/autoprof/pipeline_steps/Isophote_Fit.py

Lines 628 to 631 in 8db0fd2

	perturbations[-1][i]["ellip"] = _x_to_eps(
	_inv_x_to_eps(perturbations[-1][i]["ellip"])
	+ np.random.normal(loc=0, scale=perturb_scale)
	)

I assume that the purpose of such scaling is to do optimization in a continuous parameter space. However, since the perturbation was done in a Monte-Carlo manner, I presume that it is not critical to do such scaling. I try directly perturbing the ellipticity and when it hits its parameter boundary I simply reflect it against the boundary.

tmp_eps = perturbations[-1][i]["ellip"] + np.random.normal(
    loc=0, scale=perturb_scale
)
if tmp_eps < 0:
    tmp_eps = -tmp_eps % 1
    # tmp_eps = 0
elif tmp_eps > 1:
    tmp_eps = 1 - tmp_eps % 1
    # tmp_eps = 1
perturbations[-1][i]["ellip"] = tmp_eps

This temporary solution seems working. Let me know if you have other ideas. I would be happy to initiate a pull request when ready. Below I show a failed example and a successful fit using the new code above.

[hgao-astro]

btw, during initialization, the optimization of global ellipticity also involves _inv_x_to_eps and _x_to_eps.

AutoProf/autoprof/pipeline_steps/Isophote_Initialize.py

Lines 285 to 308 in 8db0fd2

	res = minimize(
	lambda e, d, r, p, c, n, msk: sum(
	list(
	_fitEllip_loss(_x_to_eps(e[0]), d, r * m, p, c, n, msk)
	for m in np.linspace(0.8, 1.2, 5)
	)
	),
	x0=_inv_x_to_eps(ellip),
	args=(
	dat,
	circ_ellipse_radii[-2],
	phase,
	results["center"],
	results["background noise"],
	mask,
	),
	method="Nelder-Mead",
	options={
	"initial_simplex": [
	[_inv_x_to_eps(ellip) - 1 / 15],
	[_inv_x_to_eps(ellip) + 1 / 15],
	]
	},
	)

I have no evidence this will cause an issue. Please weigh in on this. Thanks!

[ConnorStoneAstro]

Hi @hgao-astro,

Wow, this is great work! To me this seems clearly better for face-on galaxies. My only concern is for edge-on systems with ellip close to 1. In that case the shape becomes very sensitive to small changes in the ellipticity and so part of the remapping of _x_to_eps was to effectively decrease the step size near edge-on. Perhaps we could get the best of both worlds by adding an if statement which does the _x_to_eps scheme when ellip > 0.5 and doing your reflection scheme when ellip < 0.5. Both method are basically equivalent around ellip = 0.5 so it should be a pretty smooth transition. What do you think about that?

I agree for the initialization it is fine to just keep it since its more the optimization that encounters the issue.

[hgao-astro]

Thanks for the prompt reply. Let me give it a try and see if it works as expected. btw, I forgot to mention that I used a perturbation scale = 0.1 to achieve successful face-on fits, and have not yet evaluated the impact on computation time.

[hgao-astro]

Thanks for your patience. I have been occupied by job searching these days. I will try to get to it in Jan. 2024.

[ConnorStoneAstro]

@hgao-astro Good luck with your job search! I look forward to hearing from you in the new year!