lmfit-py
lmfit-py copied to clipboard
lmfit
[BUG] Caching of ci in ModelResult.conf_interval() to simple
Description
ModelResult caches the result of conf_inveral(), but it does not check if the kwargs are identical.
mod = LmfitModel()
res = mod.fit(data)
res.conf_interval()
# oh no, maxiter is to low
# Lets rerun with higher maxiter
res.conf_interval(maxiter=1000) # does nothing
@Tillsten Hm, I guess the intention is sort of obvious, but did maxiter get exceeded by the fit? Should conf_interval further exceed that?
Currently, ModelResul.conf_interval is just passing **kws to the conf_interval function. So either maxiter needs to be handled specially (and the whole, what is an iteration? thing - we can count function evaluations) or we need to know where every argument would go... I'm not sure how to best address this. Suggestions, Opinions?
I am not sure what you mean by "ModelResult caches the result of conf_inveral()" - I don't see anything special happening there or? I have never been really a fan of the apparent possibility to change parameters for CI calls. In my opinion this run is tied to the Minimizer and the MinimizerResult, so why can you not just to the fit again with the parameters you want (or just set it higher to start with from the beginning)? I doubt that there are many cases that this is such a time-consuming thing that it is really a limiting factor, but perhaps I am not understanding the issue correctly.
@reneeotten @Tillsten looking at the code for the confidence.conf_interval() function and confidence.ConfidenceInterval class more closely, they just use maxiter internally as the maximum number of iterations (per parameter) to bracket the confidence level. In one place the default is set to 50, in another place 200. ModelResult.conf_interval() should just be passing that value to confidence.conf_interval() function.
It looks like it does to me.
Adding some print statements to a reasonably well-behaved test case (where error bars were estimated), it actually needed like 5 iterations, but it was seeing the maxiter value sent to ModelResult.conf_interval().
A more complete example would be helpful.
Sorry for the intial bad bug description.
As Matt mentioned, the maxiter is the iteration used in the ci calculation, I just took one of the possible keyword arguments of the ci calculation.
The mentioned issues is that the second call does return the result from the first one. The conf_interval() methods internally caches the result in ci_out, and if available returns it. The only issue with that, that it does not take different kwargs into account. The easy workaround is to set ci_out to None.
To clarify this issue: Here's the code for ModelResult.conf_interval()
def conf_interval(self, **kwargs):
"""Calculate the confidence intervals for the variable parameters.
Confidence intervals are calculated using the
:func:`confidence.conf_interval` function and keyword arguments
(`**kwargs`) are passed to that function. The result is stored in
the :attr:`ci_out` attribute so that it can be accessed without
recalculating them.
"""
if self.ci_out is None:
self.ci_out = conf_interval(self, self, **kwargs)
return self.ci_out
The problem is that the result of the conf_interval() call is "memoized" by being stored in self.ci_out, but the signature used to generate self.ci_out is not memoized. So if conf_interval is called twice with different kwargs you will (unexpectedly, and in my opinion, incorrectly) get the same result twice.
There are a few possible solutions. The most sophisticated would be to memoize the calling signature in addition to the results. So that would be to set self.ci_out_kwargs = kwargs and check self.ci_out is None or self.ci_out_kwargs is not kwargs or something as the condition rather than just self.ci_out is None, though there may be issues with hashability of the kwargs.
A simpler brute force option would be a recalculate flag which either defaults to true or false. In this case if recalculate=False then the memoized data would be used if present. If recalculate=True the memoized data would never be used.
a third bruter solution would be to just always recalculate the confidence intervals.
My opinion: I don't think this is a good use case for memoization at all. For my use cases confidence interval calculation is never anywhere near prohibitive enough that I can't call it whenever I like. But even it was I would handle the memoization on the user end (rather than lmfit backend) by only calling conf_interval() once and storing the data and not making subsequent calls to conf_interval(). If code within lmfit is making multiple calls to conf_interval() necessitating the memoization then that might be a bigger issue. I haven't researched if this is the case.
See my suggestion which duplicates this bug for another example of surprising results due to this memoization: https://github.com/lmfit/lmfit-py/discussions/796.
Oh, in many of my usecases the CI calculation takse quite some time, still I am fine with removing the caching, since the result is still saved on the object. Alternativly, maybe just putting lru_cache decorator on the method works?
Btw, some of the problematic use cases are properly due not well tought out values in the CI calculation, since I did not really tested their influcence. In addition, the output of the ci function is rather unwidely. Adding the result to the parameters themself would also make sense in retrospect.
@jagerber48 @Tillsten The output of the calc_all_ci method is a dict with keywords of sigma and values that are the (somewhat complex) set of data for those confident intervals.
I would suggest that it is not a violation of the API or actually unreasonable to cache the values, but then always recalculate them and then update and return that cached dictionary. In that case, it's always an update. Then increasing maxiter would have an effect, and doing 1-sigma and then 2-sigma would work as expected, and the user would have to explicitly clear previous results.
Does that seem reasonable to you two?
@newville having the method save it's result into self.ci_out but then always recalculating when ModelResult.conf_interval() is called (and updating self.ci_out with the new data) seems reasonable to me.
The only thing I don't exactly understand is what you mean by "and the user would have to explicitly clear previous results", but I assume you just mean that if the user has saved the calculated confidence interval into my_conf_intervals outside the scope of lmfit they would need to overwrite their own variable my_conf_intervals with the new results in which case I agree. I just want to make sure we won't have to do ModelResult.ci_out = None to get results to update which is my current workaround. That is the user won't have to monkey with lmfit internals at all.
@jagerber48 By "and the user would have to explicitly clear previous results" I meant that if the user wanted to clear out any previous calculations, they would need to either explicitly set ModelResult.ci_out = None (and that seems kludgy) or maybe better to have a new keyword argument to ModelResult.conf_interval() - perhaps reset=True, with reset=False being the default (better suggestions for that name welcome) that would first clear any previous calculations cached in ci_out.
FWIW, I would suggest that we can declare the current behavior a bug that we want to fix, and not worry about backward compatibility.
@newville I guess I'm still not fully following your logic.
Here's my suggested form for the conf_interval() method:
def conf_interval(self, **kwargs):
"""Calculate the confidence intervals for the variable parameters.
Confidence intervals are calculated using the
:func:`confidence.conf_interval` function and keyword arguments
(`**kwargs`) are passed to that function. The result is stored in
the :attr:`ci_out` attribute so that it can be accessed without
recalculating them.
"""
self.ci_out = conf_interval(self, self, **kwargs)
return self.ci_out
I just removed the if self.ci_out is None:. The results are always cached in self.ci_out so the user can access them if they want (and they can also be accessed by ModelResult.dumps() which is the only I see self.ci_out used in the class) but they are never accessed by the conf_interval() method, which is fine if we're ok with always recalculating when conf_interval() is called.
I'm not sure what the behavior would be for your reset=False flag. Would it be this:
def conf_interval(self, reset=True, **kwargs):
ci_out = conf_interval(self, self, **kwargs)
if self.ci_out is None or reset:
self.ci_out = ci_out
return ci_out
If so the reset=False behavior seems strange to me. I'd expect the results to be updated any time the recalculation is performed. If the user wants to save old results they can do that on their end.
Just to post a third code block for reference:
def conf_interval(self, use_cached=True, **kwargs):
new_kwargs = self.conf_interval_kwargs != kwargs
if new_kwargs and use_cached:
print(f'Warning: cannot use cached values for conf_interval() because new kwargs are being used')
calculate = self.ci_out is None or new_kwargs or not use_cached
if calculate
self.ci_out = conf_interval(self, self, **kwargs)
self.conf_interval_kwargs = kwargs
return self.ci_out
This is I guess the most sophisticated option that allows lmfit to save effort sometimes but also won't give the unexpected behavior that we were seeing of using the cached value even when new kwargs are used. The unexpected behavior would not be requesting to use the cached value but actually recalculating due to new kwargs, since it's unexpected I add a warning.
Nonetheless, I still prefer the simple first option where conf_interval() never uses the cached value and any "effort saving" is the user's responsibility.
@jagerber48 I was sort of thinking of something like:
def conf_interval(self, reset=True, **kwargs):
ci_out = conf_interval(self, self, **kwargs)
if self.ci_out is None or reset:
self.ci_out = ci_out
else:
self.ci_out.update(ci_out)
return self.ci_out
Does that look reasonable? That tries to express "new calls of conf_interval() update the existing ci_out unless the user asks to reset it".
I believe that is separate from #797 (I have not looked at it closely yet: in meetings!), but I maybe they are related.
@newville Ah I see, that seems fine. For my use case I only use the sigmas kwarg so the update behavior just makes it so that if I changes the sigmas that I calculate it just updates with new ones without deleting old which is fine. I guess the only strange behavior would be if I calculate sigma=1 with one set of settings (like other kwargs) and then sigma=2 with other kwargs. Then the different sigmas would be calculated with different settings which may be confusing.
But that's getting deep into edge cases and I'm not really worried about it, especially if the user is explicitly calling for that behavior with reset=False. I'd be fine with the code snippet you posted.
I would be fine and would prefer just overwriting ci_out in case of another call to conf_inf.
@Tillsten @jagerber48 I am ok with "always overwrite" or "overwrite by default but allow updating" or "update by default but allow overwriting". I would probably put them in that order as "best", "fine", and "definitely acceptable".
I thought @jagerber48 was asking to be able to update (first do sigma=1, then do sigma=2, have a printout show both), but that may not be correct. Like, "save earlier results" could be expected to be in user code, but the printed table is so nice ;). I guess we could document that if the user saves earlier results and then does an "update" of that dictionary, they can manage the "update vs overwrite" themselves.
Because that saving of state was the previous behavior (but with the caching part sort of buggy), it seems fine to defer to that, but I think it is also OK to say that always overwriting is clearer and maybe "less surprising".
Maybe @reneeotten or someone else has a "design preference"?
@newville ah yeah maybe I wasn't clear with my example. With my first do sigma=1 then do sigma=2 I was just surprised that the printout was still showing only sigma=1 alone without any sigma=2. My expected behavior was to see the printout show only sigma=2.
Though I can see why someone might want or expect to see both sigma=1 and sigma=2. Nonetheless that's not my preference since it's just more complex, hard to control, and in my opinion, outside of the scope of lmfit. It should be up to the user to handle those sorts of effort savings.
I'm on board with the simple always overwrite solution. i.e.
def conf_interval(self, **kwargs):
self.ci_out = conf_interval(self, self, **kwargs)
return self.ci_out
As you say, the user can manage update vs overwrite themselves.
@jagerber48 Ah, thanks! So, I think we all agree that "always overwrite self.ci_out" is preferred. It's also approximately the easiest!
If you would like to make a PR for that, I would be on board. Maybe we should wait until #797 is merged?