Update temperature.py (#414)

* Update temperature.py

Improved lower bounds and initial guess for t_color

* Update bolometric.py

Improved time scale parameters for linexp and doublexp

* Update test_rainbow.py

Fixed tests

* Update temperature.py

Homogenize t_color guess for DelayedSigmoid

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Konstantin Malanchev <[email protected]>

light-curve/light_curve/light_curve_py/features/rainbow/bolometric.py

@@ -227,18 +227,14 @@ def value(t, t0, amplitude, rise_time):
     def initial_guesses(t, m, sigma, band):
 
         A = np.ptp(m)
+        med_dt = median_dt(t, band)
 
         # Compute points after or before maximum
         peak_time = t[np.argmax(m)]
         after = t[-1] - peak_time
         before = peak_time - t[0]
 
-        # Peak position as weighted centroid of everything above zero
-        idx = m > np.median(m)
-        # Weighted centroid sigma
-        dt = np.sqrt(np.sum((t[idx] - peak_time) ** 2 * m[idx] / sigma[idx]) / np.sum(m[idx] / sigma[idx]))
-        # Empirical conversion of sigma to rise/rise times
-        rise_time = dt / 2
+        rise_time = 100 * med_dt
         rise_time = rise_time if before >= after else -rise_time
 
         initial = {}
@@ -297,39 +293,36 @@ def value(t, t0, amplitude, time1, time2, p):
     @staticmethod
     def initial_guesses(t, m, sigma, band):
         A = np.ptp(m)
+        med_dt = median_dt(t, band)
 
         # Naive peak position from the highest point
         t0 = t[np.argmax(m)]
-        # Peak position as weighted centroid of everything above zero
-        idx = m > np.median(m)
-        # t0 = np.sum(t[idx] * m[idx] / sigma[idx]) / np.sum(m[idx] / sigma[idx])
-        # Weighted centroid sigma
-        dt = np.sqrt(np.sum((t[idx] - t0) ** 2 * m[idx] / sigma[idx]) / np.sum(m[idx] / sigma[idx]))
 
         # Empirical conversion of sigma to rise/fall times
-        time1 = 10 * dt
-        time2 = 10 * dt
+        time1 = 50 * med_dt
+        time2 = 50 * med_dt
 
         initial = {}
         initial["reference_time"] = t0
         initial["amplitude"] = A
         initial["time1"] = time1
         initial["time2"] = time2
-        initial["p"] = 1
+        initial["p"] = 0.1
 
         return initial
 
     @staticmethod
     def limits(t, m, sigma, band):
         t_amplitude = np.ptp(t)
         m_amplitude = np.ptp(m)
+        med_dt = median_dt(t, band)
 
         limits = {}
         limits["reference_time"] = (np.min(t) - 10 * t_amplitude, np.max(t) + 10 * t_amplitude)
         limits["amplitude"] = (0.0, 10 * m_amplitude)
-        limits["time1"] = (1e-1, 2 * t_amplitude)
-        limits["time2"] = (1e-1, 2 * t_amplitude)
-        limits["p"] = (0, 100)
+        limits["time1"] = (med_dt, 2 * t_amplitude)
+        limits["time2"] = (med_dt, 2 * t_amplitude)
+        limits["p"] = (1e-4, 10)
 
         return limits
 
@@ -338,6 +331,15 @@ def peak_time(t0, p):
         return t0 + np.real(-lambertw(p * np.exp(1)) + 1)
 
 
+def median_dt(t, band):
+    # Compute the median distance between points in each band
+    dt = []
+    for b in np.unique(band):
+        dt += list(t[band == b][1:] - t[band == b][:-1])
+    med_dt = np.median(dt)
+    return med_dt
+
+
 bolometric_terms = {
     "sigmoid": SigmoidBolometricTerm,
     "bazin": BazinBolometricTerm,

light-curve/light_curve/light_curve_py/features/rainbow/temperature.py

@@ -111,21 +111,24 @@ def value(t, t0, temp_min, temp_max, t_color):
 
     @staticmethod
     def initial_guesses(t, m, sigma, band):
+        med_dt = median_dt(t, band)
+
         initial = {}
         initial["Tmin"] = 7000.0
         initial["Tmax"] = 10000.0
-        initial["t_color"] = 1.0
+        initial["t_color"] = 10 * med_dt
 
         return initial
 
     @staticmethod
     def limits(t, m, sigma, band):
         t_amplitude = np.ptp(t)
+        med_dt = median_dt(t, band)
 
         limits = {}
         limits["Tmin"] = (1e3, 2e6)  # K
         limits["Tmax"] = (1e3, 2e6)  # K
-        limits["t_color"] = (1e-4, 10 * t_amplitude)
+        limits["t_color"] = (2 * med_dt, 10 * t_amplitude)
 
         return limits
 
@@ -160,27 +163,39 @@ def value(t, t0, Tmin, Tmax, t_color, t_delay):
 
     @staticmethod
     def initial_guesses(t, m, sigma, band):
+        med_dt = median_dt(t, band)
+
         initial = {}
         initial["Tmin"] = 7000.0
         initial["Tmax"] = 10000.0
-        initial["t_color"] = 1.0
+        initial["t_color"] = 10 * med_dt
         initial["t_delay"] = 0.0
 
         return initial
 
     @staticmethod
     def limits(t, m, sigma, band):
         t_amplitude = np.ptp(t)
+        med_dt = median_dt(t, band)
 
         limits = {}
         limits["Tmin"] = (1e3, 2e6)  # K
         limits["Tmax"] = (1e3, 2e6)  # K
-        limits["t_color"] = (1e-4, 10 * t_amplitude)
+        limits["t_color"] = (2 * med_dt, 10 * t_amplitude)
         limits["t_delay"] = (-t_amplitude, t_amplitude)
 
         return limits
 
 
+def median_dt(t, band):
+    # Compute the median distance between points in each band
+    dt = []
+    for b in np.unique(band):
+        dt += list(t[band == b][1:] - t[band == b][:-1])
+    med_dt = np.median(dt)
+    return med_dt
+
+
 temperature_terms = {
     "constant": ConstantTemperatureTerm,
     "sigmoid": SigmoidTemperatureTerm,

light-curve/tests/light_curve_py/features/test_rainbow.py

@@ -72,7 +72,7 @@ def test_noisy_all_functions_combination():
     bolometric_names = ["bazin", "sigmoid", "linexp", "doublexp"]
     bolometric_params = [bazin_parameters, sigmoid_parameters, linexp_parameters, doublexp_parameters]
 
-    Tsigmoid_parameters = [5e3, 15e3, 4.0]  # Tmin  # Tmax  # t_color
+    Tsigmoid_parameters = [5e3, 15e3, 10]  # Tmin  # Tmax  # t_color
 
     constant_parameters = [1e4]  # T