Fit unable to determiine number of fit parameters

[RJaBi]

Describe the issue:

When fitting a correlator using Obs and fits.least_squares directly, I would like to be able to have a singularly defined fit function regardless of the number of excited state terms I choose to include. Hence I might have a function such as the one below (or the vectorised form in the 'code example')

def fcn(p, x):
    """                                                                   
    sum_{alpha} A_alpha (* exp(-E_alpha * t) + exp(-E_\alpha * (NT -t)))  
    """
    NTerms = int(len(p)/2)
    A = p[0:NTerms]
    E_P = p[NTerms:]
    return (sum(ai * (np.exp(-Ei * x) + np.exp(-Ei * (NT - x)))
                for ai, Ei in zip(A, E_P))) / len(A)

This does not work. In fits.py, pyerrors tries to determine the number of fit parameters by iterating over the number of fit parameters until the function works basically. My function defined above will not raise the required IndexError or TypeError and so pyerrors quits out when it tries to check against the number of priors.

My first thought to fixing this (besides just rewriting my function which I don't really want to do) would be an optional argument (a single value or a dict for combined fits) for the number of fit parameters. Alternatively, one could just use the number of priors and/or initial parameters instead.

Code example:

def fcn(p, x):
    """                                                                                                        
    sum_{alpha} A_alpha (* exp(-E_alpha * t) + exp(-E_\alpha * (Nt -t)))                                       
    """
                                                                          
    # Assumes first half of terms are A                                                                        
    # second half are E                                                                                        
    NTerms = int(len(p)/2)
    A = np.array(p[0:NTerms])[:, np.newaxis]   # shape (n, 1)                                                  
    E_P = np.array(p[NTerms:])[:, np.newaxis]    # shape (n, 1)
    # This if statement handles the case where x is a single value rather than an array                                        
    if isinstance(x, np.float64) or isinstance(x, np.int64) or isinstance(x, float)  or isinstance(x, int):
        x = np.array([x])[np.newaxis, :]             # shape (1, m)                                            
    else:
        x = np.array(x)[np.newaxis, :]             # shape (1, m)                                              
    exp_term = np.exp(-E_P * x) + np.exp(-E_P * (NT - x))  # shape (n, m)                                      
    weighted_sum = A * exp_term                            # shape (n, m)                                      
    return np.mean(weighted_sum, axis=0)       # shape(m)

Error message:

Runtime information:

system Linux
python 3.12.3
pyerrors 2.12.0
numpy 1.26.4
scipy 1.13.0
matplotlib 3.8.4
pandas 2.2.2

[s-kuberski]

Hi Ryan,

thanks for reporting this issue. In the past, I have also encountered cases where the automatic test for the number of fit parameters created issues which I had to resolve by adapting the fit function by throwing errors if it was not called with the desired number of parameters.

So I agree that we could change something to keep the original behavior for the standard user (who does not have to think about the number of parameters at all) but still allow to set it explicitly. The two approaches that you have suggested, passing the number of parameters or inferring them from the lenght of the initial guess, also came to my mind.

I think that adding an optional parameter that allows to fix the number of parameters is the cleaner option and it could be added with a minimal change. For the other option, the user would be forced to think about what to pass as initial guess (of course, they could choose the same standard parameters as it is done inside of the function), as passing a bad initial guess could lead to problems with convergence.

If we implemented the first case, the user could still do something like

initial_guess = [1., .2, -1., .2]
n_parms = len(initial_guess)
fitp = pe.fits.least_squares(x, y, fitf, initial_guess=initial_guess, n_parms= n_parms)

to get the same behavior as in the second.

This could be quickly added to both fit routines. Any further comments @fjosw?

[fjosw]

Sounds good to me. Adding the parameter sounds like the cleaner solution to me. Do you want to have a look @s-kuberski or should I do this when I find the time?

[s-kuberski]

I can do it :)

[s-kuberski]

Hi @RJaBi ,
one more question concerning

or a dict for combined fits

As I see it, there should always only be a single number of parameters of the fit, which is the global one. Because some fit parameters are most likely shared in a combined fit, passing a dict would be ambiguous. Also, I think that the proper assignments would have to be done anyways at the stage of defining the fit functions that enter the combined fit.

But maybe I am missing something. Do you have something in mind where passing a single integer number would not resolve your issue?

[RJaBi]

Ah I see what you mean. I haven't looked into combined fits, and was merely looking at how the code loops over the keys in the dicts when determining fit parameters (but I probably have this slightly wrong cause I didn't try any).

In any case, a single integer parameter suits me fine.

Many thanks
Ryan

[s-kuberski]

We've now added the functionality in #269

[fjosw]

I just triggered a release for version 2.15 which contains this feature. You should now be able to get this new version via pypi from @RJaBi

[RJaBi]

Many thanks!

This does seem to be working and I got it from the new version (on conda-forge actually).

Thanks again,
Ryan