HG12 fitting in photometry module can produce non-physical non-monotonic phase function fits

[hhsieh00]

For some data sets, apparently particular sparse data sets sampling objects with significant lightcurve variations, fitting to a HG12 model can produce a non-monotonic phase function fit (e.g., where the maximum brightness of the object is at a non-zero phase angle), which is unphysical (maximum illumination of an object as seen by the observer should result in the maximum brightness of the object). Currently sbpy will deliver a warning if the best-fit phase function is non-monotonic, but will not enforce monotonicity for the best-fit solution delivered as output. A basic fix would be to add functionality in the photometry module to enforce monotonicity in all phase function solutions, but this appears to be a fundamental issue with the HG1G2 and HG12 systems themselves, where non-monotonic solutions are possible, at least for data sets like the ones described above, even though the relevant papers confirm that such solutions are unphysical. H. Hsieh is planning to contact Karri Muinonen and Antti Penttila (developers of the HG1G2 and HG12 systems) about this issue, and will keep J.-Y. Li in the loop.

[hhsieh00]

Actually, upon further investigation, it seems that the Penttila et al. (2016) HG1G2 and HG12 models have some enforcement of monotonicity that does not seem to be fully implemented in the sbpy photometry module, so looks like this is a potential area for improvement in sbpy. Going back to the Penttilä et al. (2016) paper, I found that they define constraints that 0 <= G1,G2 and 1-G1-G2<=1 (equation 5), and also mentions adding an additional "penalty factor" if the parameter values don't meet those constraints in Section 3.3. I also noted that their online tool uses the COBYLA algorithm (which is a "derivative-free simplex optimization" algorithm) which is not available as an option in astropy.modeling.fitting, but is available as an option in scipy.optimization.minimize, but looks like it would require some additional work to make it behave like one of the Fitter options in astropy.modeling.fitting. For what it's worth, I tried SimplexLSQFitter as the fitting routine coupled with sbpy's HG12_Pen16 model and produced similar (i.e., non-physical) results as LevMarLSQFitter. JY says he is not sure if those constraints are implemented in the current sbpy code so will take a look in the next few weeks and see if there is a straightforward way to add them if not, particularly since they have been defined already.

[jianyangli]

I just took a quick look at the relevant papers and the code. The constraints G1, G2 >= 0 and 1 - G1 - G2 <= 1 are not implemented in the current sbpy code. The fit to the models follows astropy.modeling's approach using a Fitter class. The model as implemented is not linearized as in Eq. 18 of Muinonen et al. (2010). So the fit is always performed in a non-linear fashion.

Based on the simple relationship between (G1, G2) and G12 as in Muinonen et al. (2010) and Penttilä et al. (2016), a similar constraint of 0 <= G12 <= 1 can be applied to G12, which would fulfill the constraints for G1 and G2. I tried that for the Penttilä 2016 HG12 model with a LevMarLSQFitter. But the fit would just return the best-fit parameter at the boundary (G12 = 1) for the data that would otherwise result in a non-monotonic phase function (with G12 = 4.21). The same happened no matter whether I fitted in magnitude or in reflectance scale. However, the best-fit parameter G12 = 1 did have a smaller $\chi^2$ than the best-fit parameters returned by the online tool. I think we need some input from Karri and Antti regarding how their tools work and which approach should be adopted.

[hhsieh00]

Thanks for the follow-up, Jian-Yang! Do you want to contact Karri and Antti or should I?

[jianyangli]

Let me contact them and point them to this thread.

[hhsieh00]

Just to add some more information for Karri and Antti, I have attached the data file that produces a non-monotonic best-fit solution when fit using the current sbpy HG12_Pen16 model, but results in a more reasonable-looking monotonic best-fit solution when fit using the online tool at http://h152.it.helsinki.fi/HG1G2/ , at least when applying the non-linear constrained fit option. Currently, sbpy's HG12_Pen16 implementation produces a fit that looks like the H,G1,G2 fit produced by the online tool with this data file using the linear, unconstrained fit option (turns down near zero phase angle).

data_p2019a3.txt

[anttipenttila]

Hi. The photometric systems (older HG or newer HG1G2-based ones) do not really 'bend' to all cases that easily. And yes, without constraining the parameters, one can receive best fits that are non-monotonic and therefore nonphysical. In my own works, and in the online fit tool, I usually use constrained fitting routines to ensure that the resulting fit is 'physical'. Particularly, I use 'COBYLA' algorithm in my Javascript or Fortran implementations. The limits for G1 and G2, if one wants monotonic functions, are actually a bit complicated, but if they both are between 0 and 1, and their sum also, then everything is fine. The original G12-parameter from Muinonen et al. is also a bit more cumbersome with actual limits, but again, 0<=G12<=1 will work, although small negative values can still give monotonic behavior. The newer G12* is easier, because it is designed for 0<=G12<=1 -condition.

I don't really know Astropy fitter functions that well. Personally, I have used scipy.optimize.minimize with 'COBYLA' method and with boundary conditions. Something like this, for example, can fit HG12*-system using scipy.optimize:

# Photometric system
def pf(xdeg,par1,par2):
    return sbpy.photometry.HG12_Pen16.evaluate(xdeg*np.pi/180,par1,par2)

# Sum-of-squares
def sse_fun(x,data):
    return sum([(pf(d[0],x[0],x[1])-d[1])**2 for d in data])

# Fit function
def fit_fun(data, x0=[19,0.12]):
    cv = ({'type': 'ineq', 'fun': lambda x: x[1]},
          {'type': 'ineq', 'fun': lambda x: 1-x[1]})
    res = scipy.optimize.minimize(sse_fun,x0,args=data,constraints=cv,method='COBYLA')
    return res

[anttipenttila]

BTW, from a quick look into astropy, I see that there is the SLSQPLSQFitter, which uses internally scipy.optimize SLSQP method. That should be suitable, and I have tested that my example above works fine if you change the 'COBYLA' method in minimization to 'SLSQP' method.

[mkelley]

Great! I was hoping something already implemented would work.

[anttipenttila]

I made a try to do non-constrained and constrained fits with astropy Fitters. I almost succeeded! Constraints in HG12 and HG12 work well. The HG1G2-system would need a bit more complicated constraints. The 'ineqcons' would be the correct one, but they are so poorly documented in astropy that I don't know how to use them. Attached is my example Jupyter notebook (as PDF) where I fit one dataset with astropy functions using constrained/non-constrained HG1G2, HG12, HG12 systems.

Fitting HG1G2 using astropy.pdf

[jianyangli]

Thank you for looking into this @anttipenttila ! Sorry that I didn't have time to respond earlier.

I looked through your notebook and the pdf, and agreed with your aprpoach with both scipy.optimize and astropy.modeling. Actually, your astropy approach is the same as my previous trial. There is a slightly straightforward way to add constraints to the models. The code block below is close to what I used before for testing. I couldn't figure out how to do ineqcons for HG1G2 correctly (with SLSQPLSQFitter), so I just used 0<=G1, G2<=1 below with LevMarLSQFitter.

import numpy as np
from astropy.io import ascii
from astropy.modeling.fitting import LevMarLSQFitter
from sbpy.photometry import HG1G2, HG12, HG12_Pen16

# read in Henry's data
data = ascii.read('data_p2019a3.txt')

# initialize a fitter
fitter = LevMarLSQFitter()

# fit HG12_Pen16 model
print('Fit HG12* in Penttilä 2016')
# initialize a model
m0 = HG12_Pen16(bounds={'G12': [0, 1]})
# fit model
m = fitter(m0, np.deg2rad(data['col1']), data['col2'])
# print out results
print('RMS = {:.8f}'.format((fitter.fit_info['fvec']**2).sum()))
print('Fitted H = {:.4f}, G12 = {:.6f}'.format(m.H.value, m.G12.value))
print('    meaning G1 = {:.6f}, G2 = {:.6f}\n'.format(m._G12_to_G1(m.G12), m._G12_to_G2(m.G12)))

# fit the original HG12 model
print('Fit original HG12')
m0 = HG12(bounds={'G12': [-0.0818919, 0.909714]})
m = fitter(m0, np.deg2rad(data['col1']), data['col2'])
print('RMS = {:.8f}'.format((fitter.fit_info['fvec']**2).sum()))
print('Fitted H = {:.4f}, G12 = {:.6f}'.format(m.H.value, m.G12.value))
print('    meaning G1 = {:.6f}, G2 = {:.6f}\n'.format(m._G12_to_G1(m.G12), m._G12_to_G2(m.G12)))

# fit HG1G2 model
print('Fit HG1G2')
m0 = HG1G2(bounds={'G1': [0, 1], 'G2': [0, 1]})
m = fitter(m0, np.deg2rad(data['col1']), data['col2'])
print('RMS = {:.8f}'.format((fitter.fit_info['fvec']**2).sum()))
print('Fitted H = {:.4f}, G1 = {:.6f}, G2 = {:.6f}'.format(m.H.value, m.G1.value, m.G2.value))

I think the results are close enough to what you had in the notebook. However, the problem is that the returned parameters are not the same as what your online tool would return, no matter what method is used, astropy or scipy.optimization. Still using Henry's data, below is a screenshot of what the online tool returns. Could you please let me know how your online tool implement the constrained fit? I hope we can reconcile the difference between the python approach and your online tool.

[anttipenttila]

Thanks Jian-Yang for your notes. The online fit is using the Javascript COBYLA library for constrained nonlinear minimization. I did some debugging, and it seems that in this (a bit difficult) case, the parameters for the fitting routing needed some tuning. I have now increased both the maximum number of iteration rounds, and also the initial size of the simplex in the algorithm. These changes have been uploaded to the online tool site, and you should see now that they are almost the same as with Python fitting.

[jianyangli]

Great! Thank you @anttipenttila .

I think I figured out how to make ineqcons work for astropy fitter. Basically when define the constraints in lambda functions, we should add *args in the input of the lambda func to accomodate other parameters that might be passed to the constraint function.

m0 = HG1G2(bounds={'G1': [0, 1], 'G2': [0, 1]}, ineqcons=[lambda x, *args: 1 - x[1] - x[2]])
fitter = SLSQPLSQFitter()
m = fitter(m0, phase, mag)

Note that MarLevLSQFitter doesn't support inequality constraints. So we'll have to use SLSQPLSQFitter here. I may make this fitter the default in sbpy.photometry.

Glad that we are now close to a solution. I'll try to implement it in sbpy.photometry when I get more time, and let you review it. Thanks.

[jianyangli]

@anttipenttila , @hhsieh00 , I submitted a PR #366 that added the default constraints to parameters for all four IAU HG-series models. This change should be transparent to your existing code. So you don't need to change anything, but just try to rerun what you have, and the resulting model should be constrained fit automatically. Note that you should use astropy.modeling.fitting.SLSQPLSQFitter, at least for the HG1G2 model, because LevMarLSQFitter doesn't support inequality constraints.

@anttipenttila , do you mind reviewing this PR for us? Thank you.

[anttipenttila]

Not a frequent user of advanced git features, so I don't really know how I should review. I think I was successful in pulling the correct version (gh pr checkout 366), but not sure what should I do then...

[jianyangli]

Hi @anttipenttila , at the functionality level, after you pull and checkout the branch on your local computer, you can link it to your existing python, and test it by running whatever test codes you have, such as the ipython notebooks that you shared earlier in PDF, and see how the models work now.

Also, if you are interested in doing some review at the code level, that'll be great. The simplest way to do that is to click the PR link above in this thread to open the PR, and then click "Files changed" tab right below the title. Then you can see the changes that this PR has compared to the "main" branch. You can then go through the changes and add your comments wherever you like by clicking the + sign after the line numbers. After you review, you can leave review comment by clicking the "Review changes" button on the far right of the line after the tabs buttons. You can either check "comments", "approve" or "request changes" checkbox. You can also add general comments in the "Conversation" tab like what you do here for issues.

Thank you for your input to this issue and the PR.