Quantile Regression updates

Author: robert-robison

Models/Code/quantileRegression/qr_eval2.ipynb

@@ -112,17 +112,16 @@ def _check_x(
 
 
 def _feature_transform(
-    x: np.ndarray, feature_names: list
+    x: np.ndarray, feature_names: list, convert_theta: bool = True
 ) -> Tuple[np.ndarray, np.ndarray]:
-    # Convert theta
-    idx = feature_names.index("theta")
-    x[:, idx] = np.mod(x[:, idx], 2 * np.pi)
+    if convert_theta:
+        # Convert theta
+        idx = feature_names.index("theta")
+        x[:, idx] = np.mod(x[:, idx], 2 * np.pi)
 
     # x to constant
     keep_feats = ["theta", "theta_d", "x_d", "force_in", "alpha"]
     return x, np.array([i for i, f in enumerate(feature_names) if f in keep_feats])
-    # idx = feature_names.index("x")
-    # return x, np.delete(np.arange(x.shape[1]), idx)
 
 
 class QuantileRegressionSimulator:
@@ -155,6 +154,7 @@ def __init__(
         alpha_dist_params: dict | None = None,
         random_state: Union[int, None] = None,
         smoother: str = "butterdiff",
+        convert_theta: bool = True,
     ):
         self.nrow = x.shape[0]
         self.n_feat = 2
@@ -244,6 +244,7 @@ def __init__(
         self.freq = freq
         self.smoother = smoother
         self.model_params = model_params
+        self.convert_theta = convert_theta
 
     def train(self) -> List[lgb.Booster]:
         self.models = []
@@ -253,7 +254,9 @@ def train(self) -> List[lgb.Booster]:
             y = self.datasets[i].get_label()
 
             # feature transform
-            X, f_idx = _feature_transform(X, self.feature_names_ + ["alpha"])
+            X, f_idx = _feature_transform(
+                X, self.feature_names_ + ["alpha"], convert_theta=self.convert_theta
+            )
 
             # Train data
             dtrain = lgb.Dataset(X[:, f_idx], label=y, free_raw_data=False)
@@ -299,7 +302,9 @@ def simulate_paths(
             )
 
             # Feature transform
-            x, f_idx = _feature_transform(x, self.feature_names_ + ["alpha"])
+            x, f_idx = _feature_transform(
+                x, self.feature_names_ + ["alpha"], convert_theta=self.convert_theta
+            )
 
             # Make 2nd derivative predictions
             for i in range(self.n_feat):
@@ -317,3 +322,49 @@ def simulate_paths(
                 )
 
         return preds
+
+    def predict_single(self, X: pd.DataFrame, levels=list[int]) -> pd.DataFrame:
+        # Levels to quantiles
+        alpha = [(1 - lev / 100) / 2 for lev in levels]
+        q = alpha + [1 - a for a in alpha]
+        q_names = [f"{pref}_{lev}" for pref in ["lower", "upper"] for lev in levels]
+
+        # Convert X
+        x = _check_x(X, self.feature_names_)
+
+        # Make predictions
+        res_list = []
+        for q_, q_name in zip(q, q_names):
+            x_ = np.hstack((x, np.repeat(q_, x.shape[0]).reshape(-1, 1)))
+
+            # feature transform
+            x_, f_idx = _feature_transform(
+                x_, self.feature_names_ + ["alpha"], convert_theta=self.convert_theta
+            )
+
+            for i, model in enumerate(self.models):
+                var = ["theta", "x"][i]
+                tmp = pd.DataFrame()
+
+                # Make 2nd derivative predictions
+                p = model.predict(x_[:, f_idx])
+                tmp["pred"] = p
+                tmp["variable"] = f"{var}_d"
+                tmp["quantile"] = q_name
+                tmp["t"] = X.index
+                res_list.append(tmp.copy())
+
+                # Make/propogate first derivative predictions?
+                tmp["pred"] = self.dt * X[f"{var}_d"].to_numpy()
+                tmp["variable"] = var
+                tmp["quantile"] = q_name
+                tmp["t"] = X.index
+                res_list.append(tmp)
+        all_preds = pd.concat(res_list)
+
+        # Pivot wide
+        all_preds = pd.pivot_table(
+            all_preds, index=["variable", "t"], columns="quantile", values="pred"
+        ).reset_index()
+
+        return all_preds

Models/Code/quantileRegression/qr_simulation.py

@@ -13,14 +13,15 @@
 
 outpath = Path.home() / "Box/NASA_Figures/data"
 inpath = Path.cwd() / "../../../Data/cartpoleData"
+evalpath = Path.cwd() / "../../../Results/evaluation/predictions"
 
 model = QuantileRegressionSimulator
 
 params: dict[str, Any] = {
-    "name": "qr_exp2",  # Custom unique name used for saving predictions, parameters
+    "name": "qr_exp21",  # Custom unique name used for saving predictions, parameters
     "model_name": model.__name__,
     "model_params": {
-        "m_factor": 100,
+        "m_factor": 10,
         "freq": 1 / 4,
         "alpha_dist": "beta",
         "alpha_dist_params": {
@@ -29,26 +30,34 @@
         },
         "dt": 0.01,
         "model_params": {
-            "num_iterations": 500,
-            "learning_rate": 1e-2,
+            "num_iterations": 1000,
+            "learning_rate": 1e-3,
         },
-        "smooth_derv_est": False,
-        "smoothing_samples": 100,
-        "smoothing_perc": 0.95,
+        "smooth_derv_est": True,
+        "smoothing_samples": None,
+        "smoothing_perc": 1.0,
         "smoother": "meandiff",
+        "convert_theta": True,
     },
     # Which datasets
-    "datasets": ["det"],  # det, low_noise, high_noise only options
+    "datasets": [
+        "det",
+        "low_noise",
+        "high_noise",
+    ],  # det, low_noise, high_noise only options
     # Validation parameters
     # Remainder of train always validated (unless train_seconds == 400)
     # Others must be specified
     "valid_train": True,
-    "valid_valid": False,
-    "valid_test": False,
-    "train_seconds": 400,  # Use first __ seconds of data to train, rest for val
+    "valid_valid": True,
+    "valid_test": True,
+    "train_seconds": 100,  # Use first __ seconds of data to train
+    "val_train_start": 400,  # Start val after __ seconds, same as train_seconds if None
+    "val_train_seconds": 0,  # Use __ seconds of remaining data to val, None = all
     "n_sims": 100,
     "levels": [50, 80, 95],
     "var_names": ["theta", "x", "theta_d", "x_d"],
+    "eval_modes": ["single", "multi"],
     "random_state": 6,
 }
 
@@ -65,20 +74,35 @@
         if resp.lower()[0] == "n":
             sys.exit()
 
-    # Propogate param
+    # Propogate param, set params
     params["model_params"]["random_state"] = params["random_state"]
 
-    # Read in data
+    # Loop through datasets
+    all_sim_data = []
     for dname in params["datasets"]:
         # Get all relevent datasets
         valid_sets, valid_starts = {}, {}
         data = pd.read_csv(inpath / f"{dname}_train.csv", index_col="t")
         train = data.loc[: params["train_seconds"]].copy()
 
         # Get validation sets
-        if len(train) != len(data):
-            valid_sets["val_train"] = data.loc[params["train_seconds"] :].copy()
-            skiprows = round(params["train_seconds"] * params["model_params"]["dt"])
+        if params["val_train_seconds"] is None:
+            params["val_train_seconds"] = 500 - params["train_seconds"]
+
+        if params["val_train_seconds"] > 0:
+            sp = (
+                params["train_seconds"]
+                if params["val_train_start"] is None
+                else params["val_train_start"]
+            )
+
+            # Save validation data
+            valid_sets["val_train"] = data.loc[
+                sp : (sp + params["val_train_seconds"])
+            ].copy()
+
+            # Get correct starting point
+            skiprows = round(sp / params["model_params"]["dt"])
             valid_starts["val_train"] = (
                 pd.read_csv(
                     inpath / "det_train.csv",
@@ -99,10 +123,8 @@
                 .to_numpy()
             )
         if params["valid_valid"]:
-            valid_sets["valid"] = pd.read_csv(
-                inpath / f"{dname}_val.csv", index_col="t"
-            )
-            valid_starts["valid"] = (
+            valid_sets["val"] = pd.read_csv(inpath / f"{dname}_val.csv", index_col="t")
+            valid_starts["val"] = (
                 pd.read_csv(
                     inpath / "det_val.csv",
                     nrows=1,
@@ -124,7 +146,7 @@
             )
 
         # Add valid starts to parameters
-        params["valid_starts"] = valid_starts
+        params["valid_starts"] = {k: list(v) for k, v in valid_starts.items()}
 
         # Train model
         print("Training Model")
@@ -141,47 +163,70 @@
         # Simulate over validation segments
         sim_data_list = []
         for name, val_data in valid_sets.items():
-            # Simulate trajectories
-            sims = curr_model.simulate_paths(
-                valid_starts[name],
-                force=val_data.force_in.to_numpy(),
-                n=params["n_sims"],
-                steps=val_data.shape[0],
-            )  # nsims x nsteps x 4
-
-            for i, var in enumerate(params["var_names"]):
-                # Caculate quantiles
-                sim_df = pd.DataFrame(
-                    np.quantile(sims[..., i], axis=0, q=q).T,
-                    columns=q_names,
-                    index=val_data.index,
-                )
-                sim_df["mean"] = sims[..., i].mean(axis=0)
-                sim_df["actual"] = val_data[var]
-                sim_df["name"] = name
-                sim_df["variable"] = var
-                sim_df["t"] = sim_df.index
-                sim_df = sim_df.reset_index(drop=True)
-                sim_data_list.append(sim_df)
+            for eval_mode in params["eval_modes"]:
+                if eval_mode == "multi":
+                    # Simulate trajectories
+                    sims = curr_model.simulate_paths(
+                        valid_starts[name],
+                        force=val_data.force_in.to_numpy(),
+                        n=params["n_sims"],
+                        steps=val_data.shape[0],
+                    )  # nsims x nsteps x 4
+
+                    for i, var in enumerate(params["var_names"]):
+                        # Caculate quantiles
+                        sim_df = pd.DataFrame(
+                            np.quantile(sims[..., i], axis=0, q=q).T,
+                            columns=q_names,
+                            index=val_data.index,
+                        )
+                        sim_df["mean"] = sims[..., i].mean(axis=0)
+                        sim_df["actual"] = val_data[var]
+                        sim_df["name"] = name
+                        sim_df["variable"] = var
+                        sim_df["t"] = sim_df.index
+                        sim_df["eval_mode"] = eval_mode
+                        sim_df = sim_df.reset_index(drop=True)
+                        sim_data_list.append(sim_df)
+                elif eval_mode == "single":
+                    sims = curr_model.predict_single(val_data, levels=params["levels"])
+                    for var in params["var_names"]:
+                        sim_df = sims[sims.variable == var].copy()
+                        sim_df.index = val_data.index
+                        sim_df["actual"] = val_data[var]
+                        sim_df["name"] = name
+                        sim_df["t"] = sim_df.index
+                        sim_df["eval_mode"] = eval_mode
+                        sim_df = sim_df.reset_index(drop=True)
+                        sim_data_list.append(sim_df)
 
         sim_data = pd.concat(sim_data_list)
         sim_data["noise"] = dname
 
         print("Saving predictions, parameters, and model")
 
-        # Save predictions
-        sim_data.to_csv(
-            outpath / f"validation/predictions/{params['name']}.csv", index=False
-        )
-
-        # Save parameters
-        serializable_params = make_serializable(params | {"noise": dname})
-        with open(outpath / f"validation/parameters/{params['name']}.json", "w") as f:
-            f.write(json.dumps(serializable_params, indent=4))
+        # Store predictions
+        all_sim_data.append(sim_data)
 
         # Save model
         del curr_model.datasets
         with open(
-            outpath / f"validation/model_objects/{params['name']}.pkl", "wb"
+            outpath / f"validation/model_objects/{params['name']}_{dname}.pkl", "wb"
         ) as outp:
             pickle.dump(curr_model, outp, pickle.HIGHEST_PROTOCOL)
+
+    # Save predictions
+    all_data = pd.concat(all_sim_data)
+    all_data.to_csv(
+        outpath / f"validation/predictions/{params['name']}.csv", index=False
+    )
+
+    # Save to experiment eval directory as well
+    all_data[
+        (all_data.t - all_data.groupby("name")["t"].transform("min")) <= 10
+    ].to_csv(evalpath / f"{params['name']}.csv", index=False)
+
+    # Save parameters
+    serializable_params = make_serializable(params)
+    with open(outpath / f"validation/parameters/{params['name']}.json", "w") as f:
+        f.write(json.dumps(serializable_params, indent=4))