From 477f73dd01dcde0a102df2e1d56dd9c269ad9275 Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Thu, 15 Feb 2024 03:09:04 +0100
Subject: [PATCH 01/11] Symmetric Bazin, with toggleable temperature evolution
 and rising part only

---
 .../features/rainbow/__init__.py              |   1 +
 .../features/rainbow/symmetric.py             | 200 ++++++++++++++++++
 2 files changed, 201 insertions(+)
 create mode 100644 light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py

diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/__init__.py b/light-curve/light_curve/light_curve_py/features/rainbow/__init__.py
index 9c4568d3..87bea2d2 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/__init__.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/__init__.py
@@ -1,2 +1,3 @@
 from .bazin import *
 from .rising import *
+from .symmetric import *
diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
new file mode 100644
index 00000000..d35996a0
--- /dev/null
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
@@ -0,0 +1,200 @@
+from dataclasses import dataclass
+from typing import Dict, List, Tuple
+
+import numpy as np
+
+from light_curve.light_curve_py.features.rainbow._base import BaseRainbowFit
+from light_curve.light_curve_py.features.rainbow._scaler import MultiBandScaler, Scaler
+from light_curve.light_curve_py.dataclass_field import dataclass_field
+
+__all__ = ["RainbowSymmetricFit"]
+
+
+@dataclass()
+class RainbowSymmetricFit(BaseRainbowFit):
+    """Multiband blackbody fit to the light curve using symmetric form of Bazin function
+
+    Note, that `m` and corresponded `sigma` are assumed to be flux densities.
+
+    Based on Russeil et al. 2023, arXiv:2310.02916.
+
+    Parameters
+    ----------
+    band_wave_cm : dict
+        Dictionary of band names and their effective wavelengths in cm.
+
+    with_baseline : bool, optional
+        Whether to include an offset in the fit, individual for each band.
+        If it is true, one more fit paramter per passband is added -
+        the additive constant with the same units as input flux.
+
+    with_temperature_evolution : bool, optional
+        Whether to include temperature evolution in the fit
+
+    with_rising_only : bool, optional
+        If false, only rising part of the Bazin function is used, making
+        it effectively a sigmoid function
+
+    Methods
+    -------
+    __call__(t, m, sigma, band, *, sorted=False, check=True, fill_value=None)
+        Evaluate the feature. Positional arguments are numpy arrays of the same length,
+        `band` must consist of the same strings as keys in `band_wave_cm`. If `sorted` is True,
+        `t` must be sorted in ascending order. If `check` is True, the input is checked for
+        NaNs and Infs. If `fill_value` is not None, it is used to fill the output array if
+        the feature cannot be evaluated.
+
+    model(t, band, *params)
+        Evaluate Rainbow model on the given arrays of times and bands. `*params` are
+        fit parameters, basically the output of `__call__` method but without the last
+        parameter (reduced Chi^2 of the fit). See parameter names in the `.name` attribute.
+    """
+
+    with_temperature_evolution: bool = dataclass_field(default=True, kw_only=True)
+    """Whether to include a temperature evolution in the fit."""
+    with_rise_only: bool = dataclass_field(default=False, kw_only=True)
+    """Whether to use sigmoid (rising part only) bolometric function in the fit."""
+
+    @staticmethod
+    def _common_parameter_names() -> List[str]:
+        return ["reference_time"]
+
+    def _bolometric_parameter_names(self) -> List[str]:
+        if self.with_rise_only:
+            return ["amplitude", "rise_time"]
+        else:
+            return ["amplitude", "rise_time", "fall_time"]
+
+    def _temperature_parameter_names(self) -> List[str]:
+        if self.with_temperature_evolution:
+            return ["Tmin", "Tmax", "k_sig"]
+        else:
+            return ["Tmin"]
+
+    def bol_func(self, t, params):
+        if self.with_rise_only:
+            t0, amplitude, rise_time = params[self.p.all_bol_idx]
+        else:
+            t0, amplitude, rise_time, fall_time = params[self.p.all_bol_idx]
+        dt = t - t0
+        result = np.zeros_like(dt)
+
+        # To avoid numerical overflows, let's only compute the exponents not too far from t0
+        if self.with_rise_only:
+            idx = (dt > -100*rise_time)
+            result[idx] = amplitude / (np.exp(-dt[idx] / rise_time) + 1)
+        else:
+            idx = (dt > -100*rise_time) & (dt < 100*fall_time)
+            result[idx] = amplitude / (np.exp(-dt[idx] / rise_time) + np.exp(dt[idx] / fall_time))
+
+        return result
+
+    def temp_func(self, t, params):
+        if self.with_temperature_evolution:
+            t0, T_min, T_max, k_sig = params[self.p.all_temp_idx]
+            dt = t - t0
+
+            result = np.zeros_like(dt)
+
+            # To avoid numerical overflows, let's only compute the exponent not too far from t0
+            idx1 = dt <= -100*k_sig
+            idx2 = (dt > -100*k_sig) & (dt < 100*k_sig)
+            idx3 = dt >= 100*k_sig
+
+            result[idx1] = T_min
+            result[idx2] = T_min + (T_max - T_min) / (1.0 + np.exp(dt[idx2] / k_sig))
+            result[idx3] = T_max
+        else:
+            t0, T_min = params[self.p.all_temp_idx]
+
+            return T_min
+
+    def _normalize_bolometric_flux(self, params) -> None:
+        # Internally we use amplitude of F_bol / <nu> instead of F_bol.
+        # It makes amplitude to be in the same units and the same order as
+        # the baselines and input fluxes.
+        params[self.p.amplitude] /= self.average_nu
+
+    def _denormalize_bolometric_flux(self, params) -> None:
+        params[self.p.amplitude] *= self.average_nu
+
+    def _unscale_parameters(self, params, t_scaler: Scaler, m_scaler: MultiBandScaler) -> None:
+        self._denormalize_bolometric_flux(params)
+
+        t0 = params[self.p.reference_time]
+        t0 = t_scaler.undo_shift_scale(t0)
+        params[self.p.reference_time] = t0
+
+        if self.with_rise_only:
+            amplitude, rise_time = params[self.p.bol_idx]
+            amplitude = m_scaler.undo_scale(amplitude)
+            rise_time = t_scaler.undo_scale(rise_time)
+            params[self.p.bol_idx] = amplitude, rise_time
+        else:
+            amplitude, rise_time, fall_time = params[self.p.bol_idx]
+            amplitude = m_scaler.undo_scale(amplitude)
+            rise_time = t_scaler.undo_scale(rise_time)
+            fall_time = t_scaler.undo_scale(fall_time)
+            params[self.p.bol_idx] = amplitude, rise_time, fall_time
+
+        if self.with_temperature_evolution:
+            T_min, T_max, k_sig = params[self.p.temp_idx]
+            k_sig = t_scaler.undo_scale(k_sig)
+            params[self.p.temp_idx] = T_min, T_max, k_sig
+
+    def _initial_guesses(self, t, m, band) -> Dict[str, float]:
+        t_rise = 0.1
+        t_fall = 0.1
+
+        # The amplitude here does not actually correspond to the amplitude of m values!
+        if self.with_baseline:
+            A = np.ptp(m)
+        else:
+            A = np.max(m)
+
+        params = {}
+
+        # Why do we have to strictly follow the order of parameters here?..
+        params["reference_time"] = t[np.argmax(m)]
+        params["amplitude"] = A
+        params["rise_time"] = t_rise
+
+        if not self.with_rise_only:
+            params['fall_time'] = t_fall
+
+        params['Tmin'] = 7000.0
+
+        if self.with_temperature_evolution:
+            params['Tmax'] = 0.0
+            params['k_sig'] = 1.0
+
+        return params
+
+    def _limits(self, t, m, band) -> Dict[str, Tuple[float, float]]:
+        t_amplitude = np.ptp(t)
+        if self.with_baseline:
+            m_amplitude = np.ptp(m)
+        else:
+            m_amplitude = np.max(m)
+
+        limits = {}
+
+        # Why do we have to strictly follow the order of parameters here?..
+        limits["reference_time"] = (np.min(t) - 10 * t_amplitude, np.max(t) + 10 * t_amplitude)
+        limits["amplitude"] = (0.0, 10 * m_amplitude)
+        limits["rise_time"] = (1e-4, 10 * t_amplitude)
+
+        if not self.with_rise_only:
+            limits['fall_time'] = (1e-4, 10 * t_amplitude)
+
+        limits['Tmin'] = (1e2, 1e6) # K
+
+        if self.with_temperature_evolution:
+            limits['Tmax'] = (1e2, 1e6) # K
+            limits['k_sig'] = (1e-4, 10.0 * t_amplitude)
+
+        return limits
+
+    def _baseline_initial_guesses(self, t, m, band) -> Dict[str, float]:
+        """Initial guesses for the baseline parameters."""
+        return {self.p.baseline_parameter_name(b): np.median(m[band == b]) for b in self.bands.names}

From 8a52168742da02d478af90501e6ac6d70eb340f8 Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Thu, 15 Feb 2024 03:53:21 +0100
Subject: [PATCH 02/11] Utility method for getting bolometric peak position
 added, and couple of typos fixed

---
 .../features/rainbow/symmetric.py             | 24 +++++++++++++++++--
 1 file changed, 22 insertions(+), 2 deletions(-)

diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
index d35996a0..b5db6e06 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
@@ -48,6 +48,9 @@ class RainbowSymmetricFit(BaseRainbowFit):
         Evaluate Rainbow model on the given arrays of times and bands. `*params` are
         fit parameters, basically the output of `__call__` method but without the last
         parameter (reduced Chi^2 of the fit). See parameter names in the `.name` attribute.
+
+    peak_time(*params)
+        Return bolometric peak time for given set of parameters
     """
 
     with_temperature_evolution: bool = dataclass_field(default=True, kw_only=True)
@@ -104,6 +107,8 @@ def temp_func(self, t, params):
             result[idx1] = T_min
             result[idx2] = T_min + (T_max - T_min) / (1.0 + np.exp(dt[idx2] / k_sig))
             result[idx3] = T_max
+
+            return result
         else:
             t0, T_min = params[self.p.all_temp_idx]
 
@@ -143,7 +148,10 @@ def _unscale_parameters(self, params, t_scaler: Scaler, m_scaler: MultiBandScale
             params[self.p.temp_idx] = T_min, T_max, k_sig
 
     def _initial_guesses(self, t, m, band) -> Dict[str, float]:
-        t_rise = 0.1
+        if self.with_rise_only:
+            t_rise = 1.0
+        else:
+            t_rise = 0.1
         t_fall = 0.1
 
         # The amplitude here does not actually correspond to the amplitude of m values!
@@ -165,7 +173,7 @@ def _initial_guesses(self, t, m, band) -> Dict[str, float]:
         params['Tmin'] = 7000.0
 
         if self.with_temperature_evolution:
-            params['Tmax'] = 0.0
+            params['Tmax'] = 10000.0
             params['k_sig'] = 1.0
 
         return params
@@ -198,3 +206,15 @@ def _limits(self, t, m, band) -> Dict[str, Tuple[float, float]]:
     def _baseline_initial_guesses(self, t, m, band) -> Dict[str, float]:
         """Initial guesses for the baseline parameters."""
         return {self.p.baseline_parameter_name(b): np.median(m[band == b]) for b in self.bands.names}
+
+    def peak_time(self, params) -> float:
+        """Returns true bolometric peak position for given parameters"""
+        if self.with_rise_only:
+            t0, amplitude, rise_time = params[self.p.all_bol_idx]
+
+            # It is not, strictly speaking, defined for rising only
+            return t0
+        else:
+            t0, amplitude, rise_time, fall_time = params[self.p.all_bol_idx]
+
+            return t0 + np.log(fall_time / rise_time) * rise_time * fall_time / (rise_time + fall_time)

From ae8d9d6ca9d70226a8b569e8c63b214f34974b84 Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Thu, 15 Feb 2024 04:54:47 +0100
Subject: [PATCH 03/11] Expose Minuit minimization result, and make it more
 verbose. Also, add method for getting parameters errors

---
 .../light_curve_py/features/rainbow/_base.py  | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/_base.py b/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
index 5580c2e3..a026f964 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
@@ -93,6 +93,8 @@ def __post_init__(self) -> None:
         else:
             self._lsq_model = self._lsq_model_no_baseline
 
+        self.valid = False
+
     def _initialize_minuit(self) -> None:
         self._check_iminuit()
 
@@ -269,11 +271,14 @@ def _eval(self, *, t, m, sigma, band):
             y=m,
             yerror=sigma,
         )
-        minuit = self.Minuit(least_squares, **initial_guesses)
-        minuit.migrad()
+        self.minuit = self.Minuit(least_squares, **initial_guesses)
+        # TODO: expose these parameters through function arguments
+        self.minuit.print_level = 1
+        self.minuit.migrad(ncall=10000, iterate=10)
 
-        reduced_chi2 = minuit.fval / (len(t) - self.size)
-        params = np.array(minuit.values)
+        self.valid = self.minuit.valid
+        reduced_chi2 = self.minuit.fval / (len(t) - self.size)
+        params = np.array(self.minuit.values)
 
         self._unscale_parameters(params, t_scaler, m_scaler)
         if self.with_baseline:
@@ -284,3 +289,9 @@ def _eval(self, *, t, m, sigma, band):
     # This is abstractmethod, but we could use default implementation while _eval is defined
     def _eval_and_fill(self, *, t, m, sigma, band, fill_value):
         return super()._eval_and_fill(t=t, m=m, sigma=sigma, band=band, fill_value=fill_value)
+
+    def fit_and_get_errors(self, t, m, sigma, band, *, sorted=None, check=True, fill_value=None):
+        params = self.__call__(t, m, sigma=sigma, band=band, sorted=sorted, check=check, fill_value=fill_value)
+        errors = np.array(self.minuit.errors)
+
+        return params, errors

From d2bfad85ef6cb981728b6903cac9949414863c5f Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Thu, 15 Feb 2024 06:28:13 +0100
Subject: [PATCH 04/11] Correctly unscale parameter errors

---
 .../light_curve_py/features/rainbow/_base.py     | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/_base.py b/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
index a026f964..cb5665e1 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
@@ -274,16 +274,29 @@ def _eval(self, *, t, m, sigma, band):
         self.minuit = self.Minuit(least_squares, **initial_guesses)
         # TODO: expose these parameters through function arguments
         self.minuit.print_level = 1
+        self.minuit.strategy = 2
         self.minuit.migrad(ncall=10000, iterate=10)
 
         self.valid = self.minuit.valid
         reduced_chi2 = self.minuit.fval / (len(t) - self.size)
         params = np.array(self.minuit.values)
+        self.errors = np.array(self.minuit.errors)
 
         self._unscale_parameters(params, t_scaler, m_scaler)
         if self.with_baseline:
             self._unscale_baseline_parameters(params, m_scaler)
 
+        # Unscale errors
+        # FIXME: is there any better but generic way to unscale all relevant errors without shifting?..
+        t_scaler.shift *= 0
+        m_scaler.shift *= 0
+        for _ in m_scaler.per_band_shift:
+            m_scaler.per_band_shift[_] = 0
+
+        self._unscale_parameters(self.errors, t_scaler, m_scaler)
+        if self.with_baseline:
+            self._unscale_baseline_parameters(self.errors, m_scaler)
+
         return np.r_[params, reduced_chi2]
 
     # This is abstractmethod, but we could use default implementation while _eval is defined
@@ -292,6 +305,5 @@ def _eval_and_fill(self, *, t, m, sigma, band, fill_value):
 
     def fit_and_get_errors(self, t, m, sigma, band, *, sorted=None, check=True, fill_value=None):
         params = self.__call__(t, m, sigma=sigma, band=band, sorted=sorted, check=check, fill_value=fill_value)
-        errors = np.array(self.minuit.errors)
 
-        return params, errors
+        return params, self.errors

From 6257e6df1b0aebe71bd035fcaf93df8358e33092 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 15 Feb 2024 05:33:15 +0000
Subject: [PATCH 05/11] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 .../features/rainbow/symmetric.py             | 28 +++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
index b5db6e06..7fb49457 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
@@ -3,9 +3,9 @@
 
 import numpy as np
 
+from light_curve.light_curve_py.dataclass_field import dataclass_field
 from light_curve.light_curve_py.features.rainbow._base import BaseRainbowFit
 from light_curve.light_curve_py.features.rainbow._scaler import MultiBandScaler, Scaler
-from light_curve.light_curve_py.dataclass_field import dataclass_field
 
 __all__ = ["RainbowSymmetricFit"]
 
@@ -84,10 +84,10 @@ def bol_func(self, t, params):
 
         # To avoid numerical overflows, let's only compute the exponents not too far from t0
         if self.with_rise_only:
-            idx = (dt > -100*rise_time)
+            idx = dt > -100 * rise_time
             result[idx] = amplitude / (np.exp(-dt[idx] / rise_time) + 1)
         else:
-            idx = (dt > -100*rise_time) & (dt < 100*fall_time)
+            idx = (dt > -100 * rise_time) & (dt < 100 * fall_time)
             result[idx] = amplitude / (np.exp(-dt[idx] / rise_time) + np.exp(dt[idx] / fall_time))
 
         return result
@@ -100,9 +100,9 @@ def temp_func(self, t, params):
             result = np.zeros_like(dt)
 
             # To avoid numerical overflows, let's only compute the exponent not too far from t0
-            idx1 = dt <= -100*k_sig
-            idx2 = (dt > -100*k_sig) & (dt < 100*k_sig)
-            idx3 = dt >= 100*k_sig
+            idx1 = dt <= -100 * k_sig
+            idx2 = (dt > -100 * k_sig) & (dt < 100 * k_sig)
+            idx3 = dt >= 100 * k_sig
 
             result[idx1] = T_min
             result[idx2] = T_min + (T_max - T_min) / (1.0 + np.exp(dt[idx2] / k_sig))
@@ -168,13 +168,13 @@ def _initial_guesses(self, t, m, band) -> Dict[str, float]:
         params["rise_time"] = t_rise
 
         if not self.with_rise_only:
-            params['fall_time'] = t_fall
+            params["fall_time"] = t_fall
 
-        params['Tmin'] = 7000.0
+        params["Tmin"] = 7000.0
 
         if self.with_temperature_evolution:
-            params['Tmax'] = 10000.0
-            params['k_sig'] = 1.0
+            params["Tmax"] = 10000.0
+            params["k_sig"] = 1.0
 
         return params
 
@@ -193,13 +193,13 @@ def _limits(self, t, m, band) -> Dict[str, Tuple[float, float]]:
         limits["rise_time"] = (1e-4, 10 * t_amplitude)
 
         if not self.with_rise_only:
-            limits['fall_time'] = (1e-4, 10 * t_amplitude)
+            limits["fall_time"] = (1e-4, 10 * t_amplitude)
 
-        limits['Tmin'] = (1e2, 1e6) # K
+        limits["Tmin"] = (1e2, 1e6)  # K
 
         if self.with_temperature_evolution:
-            limits['Tmax'] = (1e2, 1e6) # K
-            limits['k_sig'] = (1e-4, 10.0 * t_amplitude)
+            limits["Tmax"] = (1e2, 1e6)  # K
+            limits["k_sig"] = (1e-4, 10.0 * t_amplitude)
 
         return limits
 

From 85a2da74a599e3e254102fde9edb02b5bb987e41 Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Fri, 16 Feb 2024 12:31:01 +0100
Subject: [PATCH 06/11] Make the fitter pure, removing side-effects introduced
 earlier

---
 .../light_curve_py/features/_base.py          |  4 +-
 .../light_curve_py/features/rainbow/_base.py  | 53 ++++++++++---------
 .../features/rainbow/_scaler.py               | 11 ++++
 3 files changed, 41 insertions(+), 27 deletions(-)

diff --git a/light-curve/light_curve/light_curve_py/features/_base.py b/light-curve/light_curve/light_curve_py/features/_base.py
index cd3be57e..c74f6224 100644
--- a/light-curve/light_curve/light_curve_py/features/_base.py
+++ b/light-curve/light_curve/light_curve_py/features/_base.py
@@ -37,7 +37,7 @@ def _eval_and_fill(self, *, t, m, sigma, band, fill_value):
             if np.any(~np.isfinite(a)):
                 raise ValueError
             return a
-        except (ValueError, ZeroDivisionError) as e:
+        except (ValueError, ZeroDivisionError, RuntimeError) as e:
             if fill_value is not None:
                 return np.full(self.size, fill_value)
             raise e
@@ -142,7 +142,7 @@ def _eval_and_fill_single_band(self, *, t, m, sigma, fill_value):
             if np.any(~np.isfinite(a)):
                 raise ValueError
             return a
-        except (ValueError, ZeroDivisionError) as e:
+        except (ValueError, ZeroDivisionError, RuntimeError) as e:
             if fill_value is not None:
                 return np.full(self.size_single_band, fill_value)
             raise e
diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/_base.py b/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
index cb5665e1..6cec1e34 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/_base.py
@@ -93,8 +93,6 @@ def __post_init__(self) -> None:
         else:
             self._lsq_model = self._lsq_model_no_baseline
 
-        self.valid = False
-
     def _initialize_minuit(self) -> None:
         self._check_iminuit()
 
@@ -246,6 +244,14 @@ def _baseline_limits(self, t, m, band) -> Dict[str, Tuple[float, float]]:
         return limits
 
     def _eval(self, *, t, m, sigma, band):
+        params, errors = self._eval_and_get_errors(t=t, m=m, sigma=sigma, band=band)
+        return params
+
+    # This is abstractmethod, but we could use default implementation while _eval is defined
+    def _eval_and_fill(self, *, t, m, sigma, band, fill_value):
+        return super()._eval_and_fill(t=t, m=m, sigma=sigma, band=band, fill_value=fill_value)
+
+    def _eval_and_get_errors(self, *, t, m, sigma, band, print_level=None):
         # Initialize data scalers
         t_scaler = Scaler.from_time(t)
         m_scaler = MultiBandScaler.from_flux(m, band, with_baseline=self.with_baseline)
@@ -271,39 +277,36 @@ def _eval(self, *, t, m, sigma, band):
             y=m,
             yerror=sigma,
         )
-        self.minuit = self.Minuit(least_squares, **initial_guesses)
+        minuit = self.Minuit(least_squares, **initial_guesses)
         # TODO: expose these parameters through function arguments
-        self.minuit.print_level = 1
-        self.minuit.strategy = 2
-        self.minuit.migrad(ncall=10000, iterate=10)
+        if print_level is not None:
+            minuit.print_level = print_level
+        minuit.strategy = 2
+        minuit.migrad(ncall=10000, iterate=10)
 
-        self.valid = self.minuit.valid
-        reduced_chi2 = self.minuit.fval / (len(t) - self.size)
-        params = np.array(self.minuit.values)
-        self.errors = np.array(self.minuit.errors)
+        if not minuit.valid:
+            raise RuntimeError("Fitting failed")
+
+        reduced_chi2 = minuit.fval / (len(t) - self.size)
+        params = np.array(minuit.values)
+        errors = np.array(minuit.errors)
 
         self._unscale_parameters(params, t_scaler, m_scaler)
         if self.with_baseline:
             self._unscale_baseline_parameters(params, m_scaler)
 
         # Unscale errors
-        # FIXME: is there any better but generic way to unscale all relevant errors without shifting?..
-        t_scaler.shift *= 0
-        m_scaler.shift *= 0
-        for _ in m_scaler.per_band_shift:
-            m_scaler.per_band_shift[_] = 0
+        # We need to modify original scalers to only apply the scale, not shifts, to the errors
+        t_scaler.reset_shift()
+        m_scaler.reset_shift()
 
-        self._unscale_parameters(self.errors, t_scaler, m_scaler)
+        self._unscale_parameters(errors, t_scaler, m_scaler)
         if self.with_baseline:
-            self._unscale_baseline_parameters(self.errors, m_scaler)
-
-        return np.r_[params, reduced_chi2]
+            self._unscale_baseline_parameters(errors, m_scaler)
 
-    # This is abstractmethod, but we could use default implementation while _eval is defined
-    def _eval_and_fill(self, *, t, m, sigma, band, fill_value):
-        return super()._eval_and_fill(t=t, m=m, sigma=sigma, band=band, fill_value=fill_value)
+        return np.r_[params, reduced_chi2], errors
 
-    def fit_and_get_errors(self, t, m, sigma, band, *, sorted=None, check=True, fill_value=None):
-        params = self.__call__(t, m, sigma=sigma, band=band, sorted=sorted, check=check, fill_value=fill_value)
+    def fit_and_get_errors(self, t, m, sigma, band, *, sorted=None, check=True, **kwargs):
+        t, m, sigma, band = self._normalize_input(t=t, m=m, sigma=sigma, band=band, sorted=sorted, check=check)
 
-        return params, self.errors
+        return self._eval_and_get_errors(t=t, m=m, sigma=sigma, band=band, **kwargs)
diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py b/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py
index f5c64f5b..22e97d02 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py
@@ -47,6 +47,10 @@ def do_scale(self, x):
     def undo_scale(self, x):
         return x * self.scale
 
+    def reset_shift(self):
+        """Resets scaler shift to zero, keeping only the scale"""
+        self.shift *= 0
+
 
 @dataclass()
 class MultiBandScaler(Scaler):
@@ -83,3 +87,10 @@ def from_flux(cls, flux, band, *, with_baseline: bool) -> "MultiBandScaler":
 
     def undo_shift_scale_band(self, x, band):
         return x * self.per_band_scale[band] + self.per_band_shift[band]
+
+    def reset_shift(self):
+        """Resets scaler shift to zero, keeping only the scale"""
+        for _ in self.per_band_shift:
+            self.per_band_shift[_] = 0
+
+        super().reset_shift()

From f423c78f019068d178aa6b80e7aabf1e568b2c9e Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Tue, 20 Feb 2024 16:59:24 +0100
Subject: [PATCH 07/11] Changes from PR review

---
 .../light_curve_py/features/rainbow/_scaler.py       |  4 ++--
 .../light_curve_py/features/rainbow/symmetric.py     | 12 ++++++------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py b/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py
index 22e97d02..91f80154 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/_scaler.py
@@ -90,7 +90,7 @@ def undo_shift_scale_band(self, x, band):
 
     def reset_shift(self):
         """Resets scaler shift to zero, keeping only the scale"""
-        for _ in self.per_band_shift:
-            self.per_band_shift[_] = 0
+        for band in self.per_band_shift:
+            self.per_band_shift[band] = 0
 
         super().reset_shift()
diff --git a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
index 7fb49457..40dfb451 100644
--- a/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
+++ b/light-curve/light_curve/light_curve_py/features/rainbow/symmetric.py
@@ -109,10 +109,10 @@ def temp_func(self, t, params):
             result[idx3] = T_max
 
             return result
-        else:
-            t0, T_min = params[self.p.all_temp_idx]
 
-            return T_min
+        t0, T_min = params[self.p.all_temp_idx]
+
+        return T_min
 
     def _normalize_bolometric_flux(self, params) -> None:
         # Internally we use amplitude of F_bol / <nu> instead of F_bol.
@@ -214,7 +214,7 @@ def peak_time(self, params) -> float:
 
             # It is not, strictly speaking, defined for rising only
             return t0
-        else:
-            t0, amplitude, rise_time, fall_time = params[self.p.all_bol_idx]
 
-            return t0 + np.log(fall_time / rise_time) * rise_time * fall_time / (rise_time + fall_time)
+        t0, amplitude, rise_time, fall_time = params[self.p.all_bol_idx]
+
+        return t0 + np.log(fall_time / rise_time) * rise_time * fall_time / (rise_time + fall_time)

From 0ef8d4f694e5bc165be976717d4cc9652e17b519 Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Tue, 20 Feb 2024 17:58:18 +0100
Subject: [PATCH 08/11] Tests added

---
 .../features/test_rainbow_symmetric.py        | 185 ++++++++++++++++++
 1 file changed, 185 insertions(+)
 create mode 100644 light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py

diff --git a/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py b/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
new file mode 100644
index 00000000..19b196f8
--- /dev/null
+++ b/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
@@ -0,0 +1,185 @@
+import sys
+
+import numpy as np
+import pytest
+
+from light_curve.light_curve_py import RainbowSymmetricFit
+
+
+def bb_nu(wave_aa, T):
+    nu = 3e10 / (wave_aa * 1e-8)
+    return 2 * 6.626e-27 * nu**3 / 3e10**2 / np.expm1(6.626e-27 * nu / (1.38e-16 * T))
+
+
+@pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
+def test_noisy_no_baseline():
+    rng = np.random.default_rng(0)
+
+    band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
+
+    reference_time = 60000.0
+    amplitude = 1.0
+    rise_time = 5.0
+    fall_time = 30.0
+    Tmin = 5e3
+    Tmax = 15e3
+    k_sig = 4.0
+
+    t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, 1000))
+    band = rng.choice(list(band_wave_aa), size=len(t))
+    waves = np.array([band_wave_aa[b] for b in band])
+
+    temp = Tmin + (Tmax - Tmin) / (1.0 + np.exp((t - reference_time) / k_sig))
+    lum = amplitude / (np.exp(-(t - reference_time) / rise_time) + np.exp((t - reference_time) / fall_time))
+
+    flux = np.pi * bb_nu(waves, temp) / (5.67e-5 * temp**4) * lum
+    # S/N = 5 for minimum flux, scale for Poisson noise
+    flux_err = np.sqrt(flux * np.min(flux) / 5.0)
+    flux += rng.normal(0.0, flux_err)
+
+    feature = RainbowSymmetricFit.from_angstrom(band_wave_aa, with_baseline=False)
+
+    expected = [reference_time, amplitude, rise_time, fall_time, Tmin, Tmax, k_sig, 1.0]
+    actual = feature(t, flux, sigma=flux_err, band=band)
+
+    np.testing.assert_allclose(actual, expected, rtol=0.1)
+
+
+@pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
+def test_noisy_with_baseline():
+    rng = np.random.default_rng(0)
+
+    band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
+    average_nu = np.mean([3e10 / (w * 1e-8) for w in band_wave_aa.values()])
+
+    reference_time = 60000.0
+    amplitude = 1.0
+    rise_time = 10.0
+    fall_time = 30.0
+    Tmin = 5e3
+    Tmax = 15e3
+    k_sig = 4.0
+    baselines = {b: rng.exponential(scale=3 * amplitude / average_nu) for b in band_wave_aa}
+
+    t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, 1000))
+    band = rng.choice(list(band_wave_aa), size=len(t))
+    waves = np.array([band_wave_aa[b] for b in band])
+
+    temp = Tmin + (Tmax - Tmin) / (1.0 + np.exp((t - reference_time) / k_sig))
+    lum = amplitude / (np.exp(-(t - reference_time) / rise_time) + np.exp((t - reference_time) / fall_time))
+
+    flux = np.pi * bb_nu(waves, temp) / (5.67e-5 * temp**4) * lum + np.array([baselines[b] for b in band])
+    # S/N = 5 for minimum flux, scale for Poisson noise
+    # We make noise a bit smaller because the optimization is not perfect
+    flux_err = 0.1 * np.sqrt(flux * np.min(flux) / 5.0)
+    flux += rng.normal(0.0, flux_err)
+
+    feature = RainbowSymmetricFit.from_angstrom(band_wave_aa, with_baseline=True)
+
+    expected = [reference_time, amplitude, rise_time, fall_time, Tmin, Tmax, k_sig, *baselines.values(), 1.0]
+    actual = feature(t, flux, sigma=flux_err, band=band)
+
+    np.testing.assert_allclose(actual, expected, rtol=0.1)
+
+@pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
+def test_noisy_without_temperature_evolution():
+    rng = np.random.default_rng(0)
+
+    band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
+    average_nu = np.mean([3e10 / (w * 1e-8) for w in band_wave_aa.values()])
+
+    reference_time = 60000.0
+    amplitude = 1.0
+    rise_time = 10.0
+    fall_time = 30.0
+    Tmin = 10e3
+    baselines = {b: rng.exponential(scale=3 * amplitude / average_nu) for b in band_wave_aa}
+
+    t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, 1000))
+    band = rng.choice(list(band_wave_aa), size=len(t))
+    waves = np.array([band_wave_aa[b] for b in band])
+
+    temp = Tmin
+    lum = amplitude / (np.exp(-(t - reference_time) / rise_time) + np.exp((t - reference_time) / fall_time))
+
+    flux = np.pi * bb_nu(waves, temp) / (5.67e-5 * temp**4) * lum + np.array([baselines[b] for b in band])
+    # S/N = 5 for minimum flux, scale for Poisson noise
+    # We make noise a bit smaller because the optimization is not perfect
+    flux_err = 0.1 * np.sqrt(flux * np.min(flux) / 5.0)
+    flux += rng.normal(0.0, flux_err)
+
+    feature = RainbowSymmetricFit.from_angstrom(band_wave_aa, with_baseline=True, with_temperature_evolution=False)
+
+    expected = [reference_time, amplitude, rise_time, fall_time, Tmin, *baselines.values(), 1.0]
+    actual = feature(t, flux, sigma=flux_err, band=band)
+
+    np.testing.assert_allclose(actual, expected, rtol=0.1)
+
+@pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
+def test_noisy_with_rise_only():
+    rng = np.random.default_rng(0)
+
+    band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
+    average_nu = np.mean([3e10 / (w * 1e-8) for w in band_wave_aa.values()])
+
+    reference_time = 60000.0
+    amplitude = 1.0
+    rise_time = 10.0
+    Tmin = 5e3
+    Tmax = 15e3
+    k_sig = 4.0
+    baselines = {b: rng.exponential(scale=3 * amplitude / average_nu) for b in band_wave_aa}
+
+    t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * rise_time, 1000))
+    band = rng.choice(list(band_wave_aa), size=len(t))
+    waves = np.array([band_wave_aa[b] for b in band])
+
+    temp = Tmin + (Tmax - Tmin) / (1.0 + np.exp((t - reference_time) / k_sig))
+    lum = amplitude / (np.exp(-(t - reference_time) / rise_time) + 1)
+
+    flux = np.pi * bb_nu(waves, temp) / (5.67e-5 * temp**4) * lum + np.array([baselines[b] for b in band])
+    # S/N = 5 for minimum flux, scale for Poisson noise
+    # We make noise a bit smaller because the optimization is not perfect
+    flux_err = 0.1 * np.sqrt(flux * np.min(flux) / 5.0)
+    flux += rng.normal(0.0, flux_err)
+
+    feature = RainbowSymmetricFit.from_angstrom(band_wave_aa, with_baseline=True, with_rise_only=True)
+
+    expected = [reference_time, amplitude, rise_time, Tmin, Tmax, k_sig, *baselines.values(), 1.0]
+    actual = feature(t, flux, sigma=flux_err, band=band)
+
+    np.testing.assert_allclose(actual, expected, rtol=0.1)
+
+@pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
+def test_noisy_with_rise_only_notemp():
+    rng = np.random.default_rng(0)
+
+    band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
+    average_nu = np.mean([3e10 / (w * 1e-8) for w in band_wave_aa.values()])
+
+    reference_time = 60000.0
+    amplitude = 1.0
+    rise_time = 10.0
+    Tmin = 10e3
+    k_sig = 4.0
+    baselines = {b: rng.exponential(scale=3 * amplitude / average_nu) for b in band_wave_aa}
+
+    t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * rise_time, 1000))
+    band = rng.choice(list(band_wave_aa), size=len(t))
+    waves = np.array([band_wave_aa[b] for b in band])
+
+    temp = Tmin
+    lum = amplitude / (np.exp(-(t - reference_time) / rise_time) + 1)
+
+    flux = np.pi * bb_nu(waves, temp) / (5.67e-5 * temp**4) * lum + np.array([baselines[b] for b in band])
+    # S/N = 5 for minimum flux, scale for Poisson noise
+    # We make noise a bit smaller because the optimization is not perfect
+    flux_err = 0.1 * np.sqrt(flux * np.min(flux) / 5.0)
+    flux += rng.normal(0.0, flux_err)
+
+    feature = RainbowSymmetricFit.from_angstrom(band_wave_aa, with_baseline=True, with_rise_only=True, with_temperature_evolution=False)
+
+    expected = [reference_time, amplitude, rise_time, Tmin, *baselines.values(), 1.0]
+    actual = feature(t, flux, sigma=flux_err, band=band)
+
+    np.testing.assert_allclose(actual, expected, rtol=0.1)

From ccd7c2d324d3f22ed0c6672f3c441e57308c5a30 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 20 Feb 2024 16:59:01 +0000
Subject: [PATCH 09/11] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 .../light_curve_py/features/test_rainbow_symmetric.py      | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py b/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
index 19b196f8..d6d5b949 100644
--- a/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
+++ b/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
@@ -81,6 +81,7 @@ def test_noisy_with_baseline():
 
     np.testing.assert_allclose(actual, expected, rtol=0.1)
 
+
 @pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
 def test_noisy_without_temperature_evolution():
     rng = np.random.default_rng(0)
@@ -115,6 +116,7 @@ def test_noisy_without_temperature_evolution():
 
     np.testing.assert_allclose(actual, expected, rtol=0.1)
 
+
 @pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
 def test_noisy_with_rise_only():
     rng = np.random.default_rng(0)
@@ -150,6 +152,7 @@ def test_noisy_with_rise_only():
 
     np.testing.assert_allclose(actual, expected, rtol=0.1)
 
+
 @pytest.mark.skipif(sys.version_info < (3, 8), reason="iminuit requires Python >= 3.8")
 def test_noisy_with_rise_only_notemp():
     rng = np.random.default_rng(0)
@@ -177,7 +180,9 @@ def test_noisy_with_rise_only_notemp():
     flux_err = 0.1 * np.sqrt(flux * np.min(flux) / 5.0)
     flux += rng.normal(0.0, flux_err)
 
-    feature = RainbowSymmetricFit.from_angstrom(band_wave_aa, with_baseline=True, with_rise_only=True, with_temperature_evolution=False)
+    feature = RainbowSymmetricFit.from_angstrom(
+        band_wave_aa, with_baseline=True, with_rise_only=True, with_temperature_evolution=False
+    )
 
     expected = [reference_time, amplitude, rise_time, Tmin, *baselines.values(), 1.0]
     actual = feature(t, flux, sigma=flux_err, band=band)

From ce1bfb9e11efba046208231dcad24e2526358f3f Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Tue, 20 Feb 2024 18:07:15 +0100
Subject: [PATCH 10/11] Remove unused variable in the test

---
 .../tests/light_curve_py/features/test_rainbow_symmetric.py      | 1 -
 1 file changed, 1 deletion(-)

diff --git a/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py b/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
index d6d5b949..86c35ae2 100644
--- a/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
+++ b/light-curve/tests/light_curve_py/features/test_rainbow_symmetric.py
@@ -164,7 +164,6 @@ def test_noisy_with_rise_only_notemp():
     amplitude = 1.0
     rise_time = 10.0
     Tmin = 10e3
-    k_sig = 4.0
     baselines = {b: rng.exponential(scale=3 * amplitude / average_nu) for b in band_wave_aa}
 
     t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * rise_time, 1000))

From c56b7c0eb3a65d183828a99dbdc23a7728c946c6 Mon Sep 17 00:00:00 2001
From: karpov-sv <karpov.sv@gmail.com>
Date: Wed, 21 Feb 2024 13:32:02 +0100
Subject: [PATCH 11/11] "Fix" standard Rainbow test that stopped to converge
 with minuit.strategy=2

---
 light-curve/tests/light_curve_py/features/test_rainbow.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/light-curve/tests/light_curve_py/features/test_rainbow.py b/light-curve/tests/light_curve_py/features/test_rainbow.py
index eecd2673..e0887b68 100644
--- a/light-curve/tests/light_curve_py/features/test_rainbow.py
+++ b/light-curve/tests/light_curve_py/features/test_rainbow.py
@@ -64,12 +64,12 @@ def test_noisy_with_baseline():
     amplitude = 1.0
     rise_time = 10.0
     fall_time = 30.0
-    Tmin = 5e3
-    delta_T = 10e3
-    k_sig = 4.0
+    Tmin = 10e3
+    delta_T = 5e3
+    k_sig = 10.0
     baselines = {b: rng.exponential(scale=3 * amplitude / average_nu) for b in band_wave_aa}
 
-    t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, 1000))
+    t = np.sort(rng.uniform(reference_time - 5 * rise_time, reference_time + 5 * fall_time, 1000))
     band = rng.choice(list(band_wave_aa), size=len(t))
     waves = np.array([band_wave_aa[b] for b in band])