Hide scipy imports

Author: hombit

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

@@ -4,7 +4,6 @@
 from typing import Dict, List, Union
 
 import numpy as np
-from scipy.special import lambertw
 
 __all__ = [
     "bolometric_terms",
@@ -225,7 +224,6 @@ def value(t, t0, amplitude, rise_time):
 
     @staticmethod
     def initial_guesses(t, m, sigma, band):
-
         A = np.ptp(m)
         med_dt = median_dt(t, band)
 
@@ -328,6 +326,11 @@ def limits(t, m, sigma, band):
 
     @staticmethod
     def peak_time(t0, p):
+        try:
+            from scipy.special import lambertw
+        except ImportError:
+            raise ImportError("scipy is required for DoublexpBolometricTerm.peak_time, please install it")
+
         return t0 + np.real(-lambertw(p * np.exp(1)) + 1)
 
 

light-curve/light_curve/light_curve_py/minuit_ml.py

@@ -3,35 +3,10 @@
 from typing import Callable, Dict, Tuple
 
 import numpy as np
-from scipy.special import erf
-
-
-def logpdf(x, mu, sigma):
-    # We do not need second term as it does not depend on parameters
-    return -(((x - mu) / sigma) ** 2) / 2  # - np.log(np.sqrt(2*np.pi) * sigma)
-
-
-def barrier(x):
-    res = np.where(x > 0, 1 / x, np.inf)  # FIXME: naive barrier function
-
-    return res
-
-
-def logcdf(x):
-    # TODO: faster (maybe not so accurate, as we do not need it) implementation
-    # return norm.logcdf(x)
-
-    result = np.zeros(len(x))
-
-    idx = x < -5
-    result[idx] = -x[idx] ** 2 / 2 - 1 / x[idx] ** 2 - 0.9189385336 - np.log(-x[idx])
-    result[~idx] = np.log(0.5) + np.log1p(erf(x[~idx] / np.sqrt(2)))
-
-    return result
-
 
 try:
     from iminuit import Minuit
+    from scipy.special import erf
 except ImportError:
     MaximumLikelihood = None
 else:
@@ -63,21 +38,47 @@ def __call__(self, *par):
             ym = self.model(self.x, *par)
 
             if self.upper_mask is None:
-                result = -np.sum(logpdf(self.y, ym, self.yerror))
+                result = -np.sum(self.logpdf(self.y, ym, self.yerror))
             else:
                 # Measurements
-                result = -np.sum(logpdf(self.y[~self.upper_mask], ym[~self.upper_mask], self.yerror[~self.upper_mask]))
+                result = -np.sum(
+                    self.logpdf(self.y[~self.upper_mask], ym[~self.upper_mask], self.yerror[~self.upper_mask])
+                )
                 # Upper limits, Tobit model
                 # https://stats.stackexchange.com/questions/49443/how-to-model-this-odd-shaped-distribution-almost-a-reverse-j
                 result += -np.sum(
-                    logcdf((self.y[self.upper_mask] - ym[self.upper_mask]) / self.yerror[self.upper_mask])
+                    self.logcdf((self.y[self.upper_mask] - ym[self.upper_mask]) / self.yerror[self.upper_mask])
                 )
 
             # Barriers around parameter ranges
             # Scale is selected so that for the most of the range it is much smaller
             # than 0.5 which corresponds to 1-sigma errors
-            result += 0.0001 * np.sum(barrier((par - self.limits0) / self.limits_scale))
-            result += 0.0001 * np.sum(barrier((self.limits1 - par) / self.limits_scale))
+            result += 0.0001 * np.sum(self.barrier((par - self.limits0) / self.limits_scale))
+            result += 0.0001 * np.sum(self.barrier((self.limits1 - par) / self.limits_scale))
+
+            return result
+
+        @staticmethod
+        def logpdf(x, mu, sigma):
+            # We do not need the second term as it does not depend on parameters
+            return -(((x - mu) / sigma) ** 2) / 2  # - np.log(np.sqrt(2*np.pi) * sigma)
+
+        @staticmethod
+        def barrier(x):
+            res = np.where(x > 0, 1 / x, np.inf)  # FIXME: naive barrier function
+
+            return res
+
+        @staticmethod
+        def logcdf(x):
+            # TODO: faster (maybe not so accurate, as we do not need it) implementation
+            # return norm.logcdf(x)
+
+            result = np.zeros(len(x))
+
+            idx = x < -5
+            result[idx] = -x[idx] ** 2 / 2 - 1 / x[idx] ** 2 - 0.9189385336 - np.log(-x[idx])
+            result[~idx] = np.log(0.5) + np.log1p(erf(x[~idx] / np.sqrt(2)))
 
             return result