diff --git a/action.yml b/action.yml index 7398eba3..efc6ec63 100644 --- a/action.yml +++ b/action.yml @@ -244,7 +244,7 @@ runs: # Get static plot list (from input) read -a plots_array_input <<< ${{ inputs.plots }} # Get dynamic plot list (from comment script) - read -a plots_array_dynamic <<< "$(python draft_comment.py plots)" + read -a plots_array_dynamic <<< "$(python scripts/draft_comment.py plots)" plots_array=("${plots_array_input[@]}" "${plots_array_dynamic[@]}") @@ -294,7 +294,7 @@ runs: # Draft comment # Note: The script uses many env variables. See the script for more details. - python draft_comment.py > $GITHUB_WORKSPACE/comment.txt + python scripts/draft_comment.py > $GITHUB_WORKSPACE/comment.txt shell: bash diff --git a/draft_comment.py b/scripts/draft_comment.py similarity index 77% rename from draft_comment.py rename to scripts/draft_comment.py index 6fc50e0e..6e52f9dd 100644 --- a/draft_comment.py +++ b/scripts/draft_comment.py @@ -4,6 +4,7 @@ Script can be called via command line or imported as a module. """ +import argparse import os import re from dataclasses import dataclass @@ -12,93 +13,8 @@ import numpy as np import pandas as pd +from metrics import min_max_normalized_mae, normalized_root_mean_square_error from numpy.typing import ArrayLike -import argparse - - -def min_max_normalized_mae(y_true: ArrayLike, y_pred: ArrayLike) -> float: - """ - Calculate the min-max normalized Mean Absolute Error (MAE). - - Parameters - ---------- - y_true : array-like - True values - - y_pred : array-like - Predicted values - - Returns - ------- - float: Min-max normalized MAE - - """ - y_true = np.array(y_true) - y_pred = np.array(y_pred) - - # Ignore -inf and inf values in y_true - y_true = y_true[np.isfinite(y_true)] - y_pred = y_pred[np.isfinite(y_pred)] - - # Calculate the absolute errors - abs_errors = np.abs(y_true - y_pred) - - # Check if all errors are the same - if np.all(abs_errors == abs_errors[0]): - return 0 # Return 0 if all errors are identical to avoid division by zero - - # Min-max normalization - min_error = np.min(abs_errors) - max_error = np.max(abs_errors) - - normalized_errors = (abs_errors - min_error) / (max_error - min_error) - - return np.mean(normalized_errors) - - -def mean_absolute_percentage_error( - y_true: ArrayLike, - y_pred: ArrayLike, - epsilon: float = 1e-5, - aggregate: bool = True, - ignore_inf=True, -) -> float: - """ - Calculate the Mean Absolute Percentage Error (MAPE). - - Parameters - ---------- - y_true : array-like - True values - y_pred : array-like - Predicted values - epsilon : float, optional (default=1e-10) - Small value to avoid division by zero - - Returns - ------- - float: MAPE - - """ - y_true = np.array(y_true) - y_pred = np.array(y_pred) - - # Ignore -inf and inf values - if ignore_inf: - y_true = y_true[np.isfinite(y_true)] - y_pred = y_pred[np.isfinite(y_pred)] - - # Avoid division by zero - y_true = y_true + epsilon - y_pred = y_pred + epsilon - - # Calculate the absolute percentage errors - mape = np.abs((y_true - y_pred) / y_true) - - if aggregate: - mape = np.mean(mape) - - return mape def create_numeric_mask(arr: ArrayLike) -> np.ndarray: @@ -164,8 +80,8 @@ def errors(self, branch_type: str) -> list: ) if len(logs) != 1: msg = ( - f"Expected exactly one log file in {branch_type}.snakemake/log " - "directory." + f"Expected exactly one log file in snakemake/log directory " + "({branch_type} branch)." ) raise ValueError(msg) @@ -178,11 +94,13 @@ def errors(self, branch_type: str) -> list: @property def sucessfull_run(self) -> bool: - """Check if run was successfull via errors in logs. + """ + Check if run was successfull via errors in logs. Returns ------- bool: True if run was successfull, False otherwise. + """ if self._sucessfull_run is None: self._sucessfull_run = not bool( @@ -191,6 +109,30 @@ def sucessfull_run(self) -> bool: return self._sucessfull_run +def get_deviation_df(df1: pd.DataFrame, df2: pd.DataFrame) -> pd.DataFrame: + nrmse_series = df1.apply( + lambda row: normalized_root_mean_square_error( + row.values, + df2.loc[row.name].values, + normalization="min-max", + ), + axis=1, + ) + pearson_series = df1.apply( + lambda row: np.corrcoef(row.values, df2.loc[row.name].values)[0, 1], + axis=1, + ).fillna(0) + + if df1.empty: + return pd.DataFrame(columns=["NRMSE", "Pearson"]) + else: + deviation_df = pd.DataFrame( + {"NRMSE": nrmse_series, "Pearson": pearson_series} + ).sort_values(by="NRMSE", ascending=False) + + return deviation_df + + class RunSuccessfull(CommentData): """Class to generate successfull run component.""" @@ -222,7 +164,7 @@ def __init__(self): if plot.split(".")[-1] in ["png", "jpg", "jpeg", "svg"] ] - self._variables_deviation_ds = None + self._variables_deviation_df = None self.plots_base_url = f"https://raw.githubusercontent.com/lkstrp/pypsa-validator/{self.plots_hash}/_validation-images/" @@ -238,12 +180,12 @@ def __init__(self): STATUS_NEW = ":warning: New" VARIABLES_FILE = "KN2045_Bal_v4/ariadne/exported_variables_full.xlsx" - VARIABLES_THRESHOLD = 5 + VARIABLES_THRESHOLD = 2 @property - def variables_deviation_ds(self): - if self._variables_deviation_ds is not None: - return self._variables_deviation_ds + def variables_deviation_df(self): + if self._variables_deviation_df is not None: + return self._variables_deviation_df vars1 = pd.read_excel(self.dir_main / self.VARIABLES_FILE) vars2 = pd.read_excel(self.dir_feature / self.VARIABLES_FILE) vars1 = vars1.set_index("Variable").loc[ @@ -255,23 +197,20 @@ def variables_deviation_ds(self): assert vars1.index.equals(vars2.index) - deviation = mean_absolute_percentage_error( - vars1, vars2, ignore_inf=False, aggregate=False - ) - deviation = pd.Series( - np.round(deviation, 2).mean(axis=1), index=vars1.index - ).sort_values(ascending=False) + deviation_df = get_deviation_df(vars1, vars2) # Filter for threshold - deviation = deviation.loc[deviation > self.VARIABLES_THRESHOLD] + deviation_df = deviation_df.loc[ + deviation_df["NRMSE"] > self.VARIABLES_THRESHOLD + ] - self._variables_deviation_ds = deviation - return self._variables_deviation_ds + self._variables_deviation_df = deviation_df + return self._variables_deviation_df @property def variables_plot_strings(self): plots = ( - self.variables_deviation_ds.index.to_series() + self.variables_deviation_df.index.to_series() .apply(lambda x: re.sub(r"[ |/]", "-", x)) .apply(lambda x: "ariadne_comparison/" + x + ".png") .to_list() @@ -286,15 +225,15 @@ def variables_comparison(self) -> str: ): return "" - df = self.variables_deviation_ds.apply(lambda x: f"{x:.2f}%") - df = pd.DataFrame(df, columns=["MAPE"]) + df = self.variables_deviation_df.map(lambda x: f"{x:.3f}") df.index.name = None return ( f"{df.to_html(escape=False)}\n" f"\n" - f"MAPE: Mean Absolute Percentage Error\n" - f"Threshold: MAPE > 5%\n" + f"NRMSE: Normalized (min-max) Root Mean Square Error\n" + f"Pearson: Pearson correlation coefficient\n" + f"Threshold: NRMSE > {self.VARIABLES_THRESHOLD}\n" f"Only variables reaching the threshold are shown. Find the equivalent " f"plot for all of them below.\n\n" ) @@ -327,7 +266,6 @@ def changed_variables_plots(self) -> str: @property def plots_table(self) -> str: """Plots comparison table.""" - rows: list = [] for plot in self.plots_list: url_a = self.plots_base_url + "main/" + plot @@ -346,14 +284,39 @@ def plots_table(self) -> str: ) return df.to_html(escape=False, index=False) + "\n" - def _read_to_dataframe(self, path: Path) -> pd.DataFrame: - """Read a file to a dataframe.""" - if path.suffix == ".csv": - return pd.read_csv(path) - elif path.suffix in [".xlsx", ".xls"]: - return pd.read_excel(path) + def _format_csvs_dir_files(self, df: pd.DataFrame) -> pd.DataFrame: + """Format csvs dir dataframes.""" + if "planning_horizon" in df.columns: + df = df.T.reset_index() + # Relevant header row + header_row_index = df.index[df.iloc[:, 0] == "planning_horizon"].tolist()[0] + # Set header + df.columns = df.iloc[header_row_index] + # Fill nan values in header + df.columns = [ + f"col{i+1}" if pd.isna(col) or col == "" or col is None else col + for i, col in enumerate(df.columns) + ] + # Remove all rows before header + df = df.iloc[header_row_index + 1 :] + + # Make non-numeric values the index + non_numeric = df.apply( + lambda col: pd.to_numeric(col, errors="coerce").isna().all() + ) + + if non_numeric.any(): + df = df.set_index(df.columns[non_numeric].to_list()) else: - return None + df = df.set_index("planning_horizon") + + # Make the rest numeric + df = df.apply(pd.to_numeric) + + df.index.name = None + df.columns.name = None + + return df @property def files_table(self) -> str: @@ -368,33 +331,38 @@ def files_table(self) -> str: index_str = "/".join(str(relative_path).split("/")[1:]) path_in_feature = self.dir_feature / relative_path - if not path_in_feature.exists(): - rows[file] = [index_str, "", self.STATUS_FILE_MISSING, "", ""] + if path_in_main.parent.name == "csvs" and path_in_main.suffix == ".csv": + df1 = pd.read_csv(path_in_main) + else: continue - df1 = self._read_to_dataframe(path_in_main) - if df1 is None: + if not path_in_feature.exists(): + rows[file] = [index_str, self.STATUS_FILE_MISSING, "", ""] continue - df2 = self._read_to_dataframe(path_in_feature) + else: + df2 = pd.read_csv(path_in_feature) + + df1 = self._format_csvs_dir_files(df1) + df2 = self._format_csvs_dir_files(df2) if df1.equals(df2): - rows[file] = [index_str, "", self.STATUS_EQUAL, "", ""] + rows[file] = [index_str, self.STATUS_EQUAL, "", ""] # Numeric type mismatch elif df1.apply(pd.to_numeric, errors="coerce").equals( df2.apply(pd.to_numeric, errors="coerce") ): - rows[file] = [index_str, "", self.STATUS_TYPE_MISMATCH, "", ""] + rows[file] = [index_str, self.STATUS_TYPE_MISMATCH, "", ""] # Nan mismatch elif not df1.isna().equals(df2.isna()): - rows[file] = [index_str, "", self.STATUS_NAN_MISMATCH, "", ""] + rows[file] = [index_str, self.STATUS_NAN_MISMATCH, "", ""] # Inf mismatch elif not df1.isin([np.inf, -np.inf]).equals( df2.isin([np.inf, -np.inf]) ): - rows[file] = [index_str, "", self.STATUS_INF_MISMATCH, "", ""] + rows[file] = [index_str, self.STATUS_INF_MISMATCH, "", ""] # Changed else: # Get numeric mask @@ -411,29 +379,30 @@ def files_table(self) -> str: # Check for changes in descriptive data df1_des = df1.copy() df2_des = df2.copy() - df1_des.loc[~numeric_mask] = np.nan - df2_des.loc[~numeric_mask] = np.nan + df1_des.loc[numeric_mask] = np.nan + df2_des.loc[numeric_mask] = np.nan # Check for changes in numeric data arr1_num = pd.to_numeric(df1.to_numpy()[numeric_mask]) arr2_num = pd.to_numeric(df2.to_numpy()[numeric_mask]) - nmae = min_max_normalized_mae(arr1_num, arr2_num) - mape = mean_absolute_percentage_error(arr1_num, arr2_num) + nrmse = normalized_root_mean_square_error( + arr1_num, arr2_num, normalization="min-max" + ) + mae_n = min_max_normalized_mae(arr1_num, arr2_num) if not df1_des.equals(df2_des): status = self.STATUS_CHANGED_NON_NUMERIC - elif nmae > 0.05 and mape > 0.05: + elif nrmse > 2 and mae_n > 0.05: status = self.STATUS_CHANGED_NUMERIC else: status = self.STATUS_ALMOST_EQUAL rows[file] = [ index_str, - f"{numeric_mask.mean():.1%}", status, - f"{nmae:.2f}", - f"{mape*100:.1f}%" if mape < 1 else f"{mape*100:.2e}%", + f"{nrmse:.3f}", + f"{mae_n:.2f}", ] # Loop through all files in feature dir to check for new files @@ -443,15 +412,19 @@ def files_table(self) -> str: relative_path = os.path.relpath(path_in_feature, self.dir_feature) index_str = "../" + "/".join(str(relative_path).split("/")[1:]) - df = self._read_to_dataframe(path_in_feature) - if df is None: + if ( + path_in_feature.parent.name == "csvs" + and path_in_feature.suffix == ".csv" + ): + df1 = pd.read_csv(path_in_main) + else: continue - if not path_in_feature.exists(): - rows[file] = [index_str, "", self.STATUS_NEW, "", ""] + if not path_in_main.exists(): + rows[file] = [index_str, self.STATUS_NEW, "", ""] # Combine and sort the results - df = pd.DataFrame(rows, index=["Path", "Numeric", "Status", "NMAE", "MAPE"]).T + df = pd.DataFrame(rows, index=["Path", "Status", "NRMSE", "MAE (norm)"]).T status_order = { self.STATUS_CHANGED_NUMERIC: 0, @@ -471,9 +444,9 @@ def files_table(self) -> str: return ( f"{df.to_html(escape=False)}\n" f"\n" - f"MAPE: Mean Absolute Percentage Error\n" - f"NMAE: Mean Absolute Error on Min-Max Normalized Data\n" - f"Status Threshold: NMAE > 0.05 and MAPE > 5%\n\n" + f"NRMSE: Normalized (min-max) Root Mean Square Error\n" + f"MAE (norm): Mean Absolute Error on normalized Data (min-max\n" + f"Status Threshold: MAE (norm) > 0.05 and NRMSE > 2\n" ) @property diff --git a/scripts/metrics.py b/scripts/metrics.py new file mode 100644 index 00000000..faacd26f --- /dev/null +++ b/scripts/metrics.py @@ -0,0 +1,159 @@ +""" +Helper module for calculating evaluation metrics. +""" + +import numpy as np +from numpy.typing import ArrayLike + + +def min_max_normalized_mae(y_true: ArrayLike, y_pred: ArrayLike) -> float: + """ + Calculate the min-max normalized Mean Absolute Error (MAE). + + Parameters + ---------- + y_true : array-like + True values + + y_pred : array-like + Predicted values + + Returns + ------- + float: Min-max normalized MAE + + """ + y_true = np.array(y_true) + y_pred = np.array(y_pred) + + # Ignore -inf and inf values in y_true + y_true = y_true[np.isfinite(y_true)] + y_pred = y_pred[np.isfinite(y_pred)] + + # Calculate the absolute errors + abs_errors = np.abs(y_true - y_pred) + + # Check if all errors are the same + if np.all(abs_errors == abs_errors[0]): + return 0 # Return 0 if all errors are identical to avoid division by zero + + # Min-max normalization + min_error = np.min(abs_errors) + max_error = np.max(abs_errors) + + normalized_errors = (abs_errors - min_error) / (max_error - min_error) + + return np.mean(normalized_errors) + + +def mean_absolute_percentage_error( + y_true: ArrayLike, + y_pred: ArrayLike, + epsilon: float = 1e-9, + aggregate: bool = True, + ignore_inf=True, +) -> float: + """ + Calculate the Mean Absolute Percentage Error (MAPE). + + Parameters + ---------- + y_true : array-like + True values + y_pred : array-like + Predicted values + epsilon : float, optional (default=1e-9) + Small value to avoid division by zero + + Returns + ------- + float: MAPE + + """ + y_true = np.array(y_true) + y_pred = np.array(y_pred) + + # Ignore -inf and inf values + if ignore_inf: + y_true = y_true[np.isfinite(y_true)] + y_pred = y_pred[np.isfinite(y_pred)] + + # Avoid division by zero + y_pred = y_pred + epsilon + + # Calculate the absolute percentage errors + mape = np.abs((y_true - y_pred) / y_true) + + if aggregate: + mape = np.mean(mape) + + return mape + + +def normalized_root_mean_square_error( + y_true: ArrayLike, + y_pred: ArrayLike, + ignore_inf: bool = True, + normalization: str = "min-max", + fill_na: float = 0, + epsilon: float = 1e-9, +) -> float: + """ + Calculate the Normalized Root Mean Square Error (NRMSE). + + Parameters + ---------- + y_true : array-like + True values + y_pred : array-like + Predicted values + ignore_inf : bool, optional (default=True) + If True, ignore infinite values in the calculation + normalization : str, optional (default='min-max') + Method of normalization. Options: 'mean', 'range', 'iqr', 'min-max' + epsilon : float, optional (default=1e-9) + Small value to add to normalization factor to avoid division by zero + + Returns + ------- + float: NRMSE + + """ + y_true_ = np.array(y_true) + y_pred_ = np.array(y_pred) + + if np.array_equal(y_true_, y_pred_): + return 0 + + # Fill NaN values + y_true_ = np.nan_to_num(y_true_, nan=fill_na) + y_pred_ = np.nan_to_num(y_pred_, nan=fill_na) + + # Ignore -inf and inf values + if ignore_inf: + mask = np.isfinite(y_true_) & np.isfinite(y_pred_) + y_true_ = y_true_[mask] + y_pred_ = y_pred_[mask] + + # Calculate the squared errors + squared_errors = (y_true_ - y_pred_) ** 2 + + # Calculate RMSE + rmse = np.sqrt(np.mean(squared_errors)) + + # Normalize RMSE + if normalization == "mean": + normalization_factor = np.mean(y_true_) + elif normalization == "range": + normalization_factor = np.ptp(y_true_) + elif normalization == "iqr": + q75, q25 = np.percentile(y_true_, [75, 25]) + normalization_factor = q75 - q25 + elif normalization == "min-max": + normalization_factor = np.max(y_true_) - np.min(y_true_) + epsilon + else: + raise ValueError("Invalid normalization method.") + + nrmse = rmse / normalization_factor + + return nrmse