Update mm_rework to match mm_rework_lingfeng for comparison

[wei-lingfeng]

This pull request is just to show the updates that make mm_rework and mm_rework_lingfeng different. Since the plan is to go with mm_rework_lingfeng eventually, it does not need to be actually merged. I will open another pull request to merge mm_rework_lingfeng into dev and eventually into main.

1. startables.StarTable.combine_lists: Uncertainty Handling: Used sqrt of weighted variance as uncertainty. This value would be zero for stars with only 1 measurement. So the original code had to handle 1 epoch as a special case. But similarly, stars with 2 epochs with similar uncertainties are also going to have a close-to-zero uncertainty, as a line can go through 2 points, and so on with 3 or more epochs. This won't help us rescale the uncertainties to a larger value. Instead, it results in a smaller value with a lot more "special cases" that need to be handled. I don't think it's a good idea to make such peculiar assumptions. Upscaling uncertainties should be done outside of flystar. Change: Changed to error on the mean (error propagation), which does not have this issue. See next point for more reasons.
2. motion_model.Fixed.run_fit: Inconsistent uncertainty calculation: Fixed model uses weighted variance square root as uncertainty, while other models use scipy.curve_fit’s covariance matrix. This is inconsistent: scipy’s covariance matrix is essentially the same as error propagation, not weighted variance. (See code below as a test). To me, this is one more reason we should choose error propagation instead of weighted variance. Change: Changed everything to the error propagation equation (error on the mean).
3. motion_model.Parallax Inaccurate required points Currently it has n_pts_req=4. However, 3 pairs of x, y are enough to solve for a 5-parameter equation. If the user wants to use Parallax only for >=4 epochs, they can specify the motion_model_input column to do so. This should not be something the package imposes. Change n_pts_req=3.
4. motion_model.get_weights: Precision Problem: The back-and-forth conversion between x_wt = 1/xe**2, then sigma = 1/x_wt**0.5 would result in floating-point errors, and in multiple iterations of align, this error is cumulative. Change: Changed to motion_model.calc_sigma, which only returns sigma_x = xe or np.sqrt(xe), that can be input into scipy directly. Weight is only calculated when needed (most of the time not)
5. motion_model.Parallax.run_fit: BUG: sigma not taken square root in curve_fit. BUG: absolute_sigma is not used. Change: Now replaced everything with calc_sigma and added absolute_sigma
6. startables.StarTable.fit_velocities (now fit_motion_models): fit_star_idxs behavior: When fit_star_idxs is provided in align, i.e., only fit a subset of stars, do we want to wipe out all previous results? mm_rework will overwrite any existing values with nan. Change: I added the option to keep all existing values and only refit stars with the given index (by default)
7. parallax.parallax_in_direction: Minor Change: times=Times(mjd + 2400000.5, format=’jd’, scale=’tdb’) to times=Times(mjd, format=’mjd’, scale=’tdb’)
8. align.MosaicToRef/MosaicSelfRef: Default absolute_sigma: absolute_sigma=False previously, now absolute_sigma=True. We need to decide which one we want by default.
9. Known Issue: test_align would fail at test_MosaicSelfRef_vel. This failure is already present at this point. But mm_rework_lingfeng has no issues with this.

Example showing scipy’s curve_fit is equivalent to propagated uncertainty (error on the mean):

import numpy as np
from scipy.optimize import curve_fit

# 1. Define dummy data and explicit measurement errors
y_data = np.array([10.2, 9.8, 10.5])
sigma_y = np.array([0.1, 0.4, 0.2])   # Different errors for each point
x_dummy = np.zeros_like(y_data)       # x is required but ignored by a constant function

# --- Method A: Analytical Error Propagation ---
weights = 1.0 / (sigma_y ** 2)
c_prop = np.sum(y_data * weights) / np.sum(weights)
sigma_c_prop = np.sqrt(1.0 / np.sum(weights))

# --- Method B: SciPy Curve Fit ---
def constant_model(x, c):
    return c

# absolute_sigma=True is MANDATORY to match pure error propagation
popt, pcov = curve_fit(constant_model, x_dummy, y_data, sigma=sigma_y, absolute_sigma=True)

c_fit = popt[0]
sigma_c_fit = np.sqrt(pcov[0, 0]) # Square root of the variance element

# --- Display Results ---
print(f"Propagation:  Value = {c_prop:.5f}, Error = {sigma_c_prop:.5f}")
print(f"Curve Fit:    Value = {c_fit:.5f}, Error = {sigma_c_fit:.5f}")

Output:

Propagation:  Value = 10.23810, Error = 0.08729
Curve Fit:    Value = 10.23810, Error = 0.08729