Improve Rainbow fitting

[karpov-sv]

This PR adds new sub-class for Rainbow, RainbowSymmetricFit that implements symmetric form of Bazin function for bolometric term, and slightly modified sigmoid function for temperature term, allowing for both rising and falling temperatures:

bol = A / (exp(-dt/rise_time) + exp(dt/fall_time)) 
temp = T_min + (T_max - T_min) / (1 + exp(dt / k_sig)

These are better suited for constrained minimization, as they do not require complex constraints combining several parameters as in current implementation (not actually implemented there, leading to nonsensical results sometimes).

New class also has options for fitting fixed temperature models, as well as rising only bolometric function (effectively just a sigmoid). Both may be toggled independently.
In principle, this implementation supersedes both basin.py and rising.py, and may replace them completely.

PR also modifies the base class by exposing some internals from Minuit - new field valid reflects whether the minimization succeeded or not, and new field minuit allows directly accessing Minuit instance, e.g. to see its diagnostic output.
I took a liberty of also increasing the verbosity level of Minuit itself (print_level=1) so that it prints some warnings during the minimization when it is problematic. Also, number of iterations is increased (ncall=10000, iterate=10), and default minimization strategy is changed to slightly more stable, but probably slower one (strategy=2).
It should be, of course, made optional, but I found no easy way to pass them as additional keyword parameters for the _eval method (too many levels of abstraction above it).

Related to the last point - as there is no easy way to add extra keywords to _eval method, I had to implement additional entry point for fitting with also getting back the parameter errors (instead of just making it an option). The method is fit_and_get_errors, with identical signature to __call__, and it returns the list containing arrays of parameters and their errors.

In general, the minimization on real data is not very stable, and strongly depends on initial parameters (as it is local minimizer only). I adjusted the heuristics for initial parameters to be slightly better, but it would be desirable to be able to also directly pass them during the feature evaluation. It comes down, again, to passing extra parameters to __call__ (as well as to proper normalization of the values there, but that's easily doable).

Tests are not (yet) implemented for new class, as it will probably require some changes anyway.

[codspeed-hq]

CodSpeed Performance Report

Merging #314 will not alter performance

_{Comparing karpov-sv:rainbow-symmetric (c56b7c0) with master (bc1039d)}

Summary

✅ 97 untouched benchmarks

[hombit]

Thank you, @karpov-sv! I'll look into code a bit later, here I'll try to address some of your general points.

... new field valid...

Related to the last point - as there is no easy way to add extra keywords to _eval method

I'm sorry that you've found code structure to be too abstract. The key idea in this library that user constructs feature extractor in advance and apply it to one, or many, light curves. This means that all the feature parameters must be passed in advance, they are currently handled as dataclasses.dataclass fields (represented by keyword arguments of the constructor). It also means that feature extractor instance is supposed to be immutable and should not change during the evaluation. Could I ask you to reimplement your work in this way, so .valid is not a field anymore, and all the feature extractor parameters are attributes and passed during the construction?

... but it would be desirable to be able to also directly pass them during the feature evaluation. It comes down, again, to passing extra parameters to call (as well as to proper normalization of the values there, but that's easily doable).

I understand this point, but I against it. The library is supposed to provide the even interface over feature extractors, so they could be combined and abstracted from evaluation step. In this particular case I see no way how user would approach with boundaries individual for each light curve, not using the light curve data itself. If we can do better job estimating boundaries from light curve, it would be nice eval does it for them. However, if it still doesn't fit your needs, I'm absolutely fine with adding more specific methods to the feature class, just in the same way you did with fit_and_get_errors.

I took a liberty of also increasing the verbosity level of Minuit itself (print_level=1)

All ideas you propose about Minuit sound very valid for me. print_level is the only thing I'd ask you to make parametric and defaulted to the previous value. I understand why you want it to be there, and probably understand why iminuit just prints instead of using proper logging, but I still would prefer to make the code silent by default.

[hombit]

Scaler has undo_scale which doesn't shift

light-curve-python/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py

Lines 47 to 48 in bc1039d

	def undo_scale(self, x):
	return x * self.scale

[karpov-sv]

Yes, but there is no way to use it without knowing the underlying data - _unscale_parameters method is defined in subclasses, and there is no way to know which parameter was transformed using which scaler from the base class apart from using it. So there is no clean way to apply just undo_scale method if the subclass uses undo_shift_scale (or even undo_shift_scale_band). Well, of course apart from defining new virtual method like _unscale_parameters_errors for the subclasses to implement - but that's probably too much?..

[hombit]

Ah, I see now what did you mean, it is about parameter errors, not about flux errors. Sorry for the confusion

[karpov-sv]

Yes, I understand your design choices and I agree with them in general. So I will modify the PR accordingly - as soon as I understand how to do it better.

One more point of consideration here, related to immutability - I am no expert in Minuit, but it seems to me that print_level is a global (static per-process) setting, so setting it on one instance will influence all other instances running. So - as you are planning to use Minuit wider across the package - what is the better way to go here? Fortunately, the rest of Minuit parameters (like strategy) are per-instance, so may be changed without side effects. print_level is not extremely important by itself, so if you wish to fully avoid side effects it may be excluded completely.

What is important is to have some sort of error reporting - that's why I need .valid field at least. But probably you have better approach already implemented for handling the errors in the fitters? As it is quite generic question, so it may e.g. raise some exception?.. I see that _eval_and_fill does handle some exceptions (and fills the results with fill_value then) - is it the way to go? Should I just raise ValueError (or better RuntimeError as it fits better?.. and modify the method to handle it as well?..) when the minimization fails?

Then, if error handling will be implemented this way - I may remove all mutable fields from the PR, and move the rest of features I need to the arguments of specific methods. (E.g. I might try to add some support for initial parameters, as for some tasks you have some educated guess about them, like rough peak time, or rough rise/fall timescales - but I do not yet know whether it is possible to do it in any generic way due to all this scaling/unscaling of parameters inside).

[hombit]

One more point of consideration here, related to immutability - I am no expert in Minuit, but it seems to me that print_level is a global (static per-process) setting, so setting it on one instance will influence all other instances running. So - as you are planning to use Minuit wider across the package - what is the better way to go here?

Let's just let user set it

Should I just raise ValueError (or better RuntimeError as it fits better?.. and modify the method to handle it as well?..) when the minimization fails?

Raising exception is fine. If user set fill_value=None the exception will propagate, otherwise fill_value: float will be used as an "na" value for this feature. I usually raise ValueError, but I don't have a strong opinion here, maybe we should just make a new exception type and use it instead.

due to all this scaling/unscaling of parameters inside

We do scaling/unscaling to make all parameters to be the same order of magnitude. We consider that user could use crazy units for inputs, let's say years for time and W/cm^2/Hz for fluxes, so it would be better to normalize everything before using the optimizer. This normalization is specific for this feature, so we can reimplement it in any other way, or even drop it if we believe that optimizer could handle these different orders of magnitudes properly.

[karpov-sv]

Ok, so here is new version, with all side effects removed (apart from single optional argument for fit_and_get_errors custom method that, if provided, sets Minuit print_level globally).
On fitting error the code raises RuntimeException which is better suited than ValueException. Thus, BaseMultiBandFeature is modified to also catch it.

Minuit parameters (ncall=10000, iter=10 and strategy=2) are hardcoded as they should be mostly harmless.

I did not yet come with any generic enough way to pass initial parameters to the fitter (as they should be normalized alongside with the data, which cannot be done in base class as the underlying model is unknown, and there is no method for that in subclasses either), so I will give up on that for now. Thus, probably the code is ready for your review now.

[hombit]

Sorry for the delay and thank you!

Do you have time to add some tests for the code? It is fine to say "no", I could make it myself a bit later

[karpov-sv]

Ok, I added some tests in a similar manner to existing Rainbow ones. Also incorporated the fixes you requested.

[karpov-sv]

The old test for standard RainbowFit stopped to converge with minuit.strategy=2 so I had to adjust the testing parameters a bit. Now it should pass.

[hombit]

It is now published as a part of 0.8.2 release