PyLIMA fitting procedures return swapped fluxes for certain telescope

[KKruszynska]

I have encountered an issue where TRF_fit and MCMC_fit (and possibly others) have swapped flux parameters for the first two telescopes.
Here are excerpts from a code that I use:

1. I load telescopes in this order: Gaia_G, OGLE_I, LCO_gp, LCO_ip

telescope_order = ['01_gaia_gsa_G.dat', 
                   '02_OGLE-2024-BLG-0034.dat',
                   '03_LCO_cleaned_gp.dat', 
                   '04_LCO_cleaned_ip.dat', 
                   # '05_ATLAS_binned_o.dat',
                   # '06_ATLAS_binned_c.dat',
                   # '07_blg607.13.v.25782.dat', 
                  ]

for event_telescope in telescope_order:

    try:
        lightcurve = np.loadtxt(file_paths + event_telescope, dtype=str)
        print(event_telescope)

        if 'gsa' in event_telescope:
            tel_name = event_telescope.split('.')[0]
            lightcurve = lightcurve[:, [0, 1, 2]].astype(float)
            t, m, e = lightcurve[:, 0], lightcurve[:, 1], lightcurve[:, 2]
            name = 'Gaia_' + tel_name.split('_')[-1].split('.')[0]
            lightcurve = np.vstack((t, m, e)).T
            telescope = telescopes.Telescope(name=name,
                                             lightcurve=lightcurve,
                                             lightcurve_names=['time', 'mag', 'err_mag'],
                                             lightcurve_units=['JD', 'mag', 'mag'],
                                             location='Space', spacecraft_name='Gaia')
        elif 'LCO' in event_telescope:
            tel_name = event_telescope.split('.')[0]
            lightcurve = lightcurve[:, [0, 1, 2]].astype(float)
            t, m, e = lightcurve[:, 0], lightcurve[:, 1], lightcurve[:, 2]
            name = 'LCO_' + tel_name.split('_')[-1].split('.')[0]
            lightcurve = np.vstack((t, m, e)).T
            telescope = telescopes.Telescope(name=name,
                                             lightcurve=lightcurve,
                                             lightcurve_names=['time', 'mag', 'err_mag'],
                                             lightcurve_units=['JD', 'mag', 'mag'],
                                             location='Earth')
    
        elif 'OGLE' in event_telescope:
            lightcurve = lightcurve[:, [0, 1, 2]].astype(float)
            t, m, e = lightcurve[:, 0], lightcurve[:, 1], lightcurve[:, 2]
            name = 'OGLE_I'
            lightcurve = np.vstack((t, m, e)).T
            telescope = telescopes.Telescope(name=name,
                                             lightcurve=lightcurve,
                                             lightcurve_names=['time', 'mag', 'err_mag'],
                                             lightcurve_units=['JD', 'mag', 'mag'],
                                             location='Earth')

        if (len(lightcurve) > 2):  # & ('gp' not in name):
            telescope.ld_a1 = 0.57
            telescope.define_limb_darkening_coefficients()
            my_event.telescopes.append(telescope)

        ite += 1
        print(name, np.min(lightcurve[:, 0]), np.min(lightcurve[:, 2]))

    except Exception as err:
        print(f'Unexpected {err=}, {type(err)=}')
        raise

2. I set up the event, and the model (USBL or PSBL), and I set the reference survey to OGLE_I (second telescope):

fancy = pyLIMA_fancy_parameters.StandardFancyParameters()
usbl = USBL_model.USBLmodel(my_event, fancy_parameters=fancy, parallax=['Full', 2460538.0],  blend_flux_parameter='ftotal')
my_event.name = 'OGLE_2024_BLG_0034_TRF_outer'
my_event.find_survey('OGLE_I')

3. I set up a fit with starting params, and I run the fit:

trf_fit = TRF_fit.TRFfit(usbl)
trf_fit.model_parameters_guess = [some starting parameters]
trf_fit.fit()

5. The resulting fit has good parameters and fluxes, but OGLE_I fluxes appear above Gaia_G fluxes in the results trf_fit.fit_results['best_model'].

This affects the plots, with OGLE_I and Gaia_G being separated from the plotted model, but LCO_gp and LCO_ip being aligned with it. See plot below:

When I manually swapped fluxes associated with Gaia_G and OGLE_I in the results tuple, I got the correct plot and a decent value of chi2. I am not sure if this is occurring for PSPL models, or if when I set find_survey to the first telescope.

[KKruszynska]

I checked. When I set find_survey to the first telescope, the issue doesn't appear.