Merge pull request #15 from AlexHls/fix-model-syntax

Fix model syntax

Author: AlexHls

sccala-env.txt

@@ -5,7 +5,7 @@ george
 emcee
 astropy
 specutils
-pystan
+pystan==3.9.0
 corner
 cloudpickle
 itsdangerous

src/sccala/collect_scm_data.py

@@ -22,6 +22,9 @@ def main(args):
     vel_sys = args.vel_sys
     col_sys = args.col_sys
     ae_sys = args.ae_sys
+    rho = args.rho
+    rho_calib = args.rho_calib
+    error_mode = args.error_mode
 
     df = sl.load_data(
         sne_list,
@@ -37,6 +40,9 @@ def main(args):
         vel_sys=vel_sys,
         col_sys=col_sys,
         ae_sys=ae_sys,
+        rho=rho,
+        rho_calib=rho_calib,
+        error_mode=error_mode,
     )
 
     return df
@@ -115,6 +121,25 @@ def cli():
         help="Value of the systematic a/e uncertainty.",
         type=float,
     )
+    parser.add_argument(
+        "-r",
+        "--rho",
+        help="Correlation between the color and magnitude uncertainties. Default: 1.0",
+        default=1.0,
+        type=float,
+    )
+    parser.add_argument(
+        "--rho_calib",
+        help="Correlation between the color and magnitude uncertainties for calibrator SNe. Default: 0.0",
+        default=0.0,
+        type=float,
+    )
+    parser.add_argument(
+        "--error_mode",
+        help="Mode to calculate asymmetric errors. Default: mean",
+        default="mean",
+        choices=["mean", "max", "min"],
+    )
 
     args = parser.parse_args()
 

src/sccala/interplib/epoch_interp.py

@@ -86,7 +86,7 @@ def __init__(
 
         self.data = np.array(data)
         self.data_error = np.array(data_error)
-        self.tkde = np.array(tkde)
+        self.tkde = tkde
         self.mjd = np.array(mjd)
 
         self.snname = snname
@@ -97,7 +97,7 @@ def __init__(
         self.reg_max = reg_max
         self.extrapolate = extrapolate
 
-        self.toe = float(np.percentile(tkde, 50.0))
+        self.toe = float(np.percentile(tkde.resample(10000), 50.0))
 
         # Confert dates to restframe
         self.time = (mjd - self.toe) / (1 + red)
@@ -216,7 +216,7 @@ def diagnostic_plot(self, diagnostic, target, flux_interp=False):
                 color="red",
                 label="Median",
             )
-            ax1.set_xlabel("v(t={:.1f})d) (km/s)".format(self.dates[plotind]))
+            ax1.set_xlabel("v(t={:.1f}d) (km/s)".format(self.dates[plotind]))
             ax1.set_title(
                 r"v({:.1f}) = {:.2f} +{:.2f}/ -{:.2f} km/s | {:s}".format(
                     self.dates[plotind],
@@ -242,8 +242,16 @@ def diagnostic_plot(self, diagnostic, target, flux_interp=False):
                 color="red",
             )
             ax2.axhline(self.median[plotind] / conv, color="red")
-            lower = self.dates[plotind] + self.toe - np.percentile(self.tkde, 15.87)
-            upper = self.dates[plotind] + self.toe - np.percentile(self.tkde, 84.13)
+            lower = (
+                self.dates[plotind]
+                + self.toe
+                - np.percentile(self.tkde.resample(10000), 15.87)
+            )
+            upper = (
+                self.dates[plotind]
+                + self.toe
+                - np.percentile(self.tkde.resample(10000), 84.13)
+            )
             ax2.axvspan(
                 lower, upper, alpha=0.3, color="blue", label="1$\sigma$ (68.26%)"
             )