# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Tutorial: Using DFtoVW and exploring VW output
#
# ---
#
# The purpose of this tutorial is to show how to use the `DFtoVW` class to convert a pandas's DataFrame into a list of Vowpal Wabbit examples and to explore the outputs (model weights, VW output log) of the trained model. The VW output log is parsed using the class `VWLogParser` defined in this notebook.
#
# The task is to predict the concentration of [particulate matter](https://en.wikipedia.org/wiki/Particulates) (more specifically PM 2.5) in the atmosphere of 5 chinese cities. The original dataset contained 19 columns (targets, datetime and atmospheric features) and 167 358 observations.
# - For more details on the data, see the following UCI repository: https://archive.ics.uci.edu/ml/datasets/PM2.5+Data+of+Five+Chinese+Cities
# - For the associated academic papers, see [<NAME>, et al. "PM2. 5 data reliability, consistency, and air quality assessment in five Chinese cities."](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016JD024877) and [<NAME>, et al. "Assessing Beijing's PM2. 5 pollution: severity, weather impact, APEC and winter heating."](https://royalsocietypublishing.org/doi/10.1098/rspa.2015.0257)
#
# The data can be download from the following URL: https://archive.ics.uci.edu/ml/machine-learning-databases/00394/FiveCitiePMData.rar. A function `download_data` is available in this notebook to download and extract the data (but this step can be done manually if desired). The folder is containing the data is defined using the constant variable `DATA_FOLDER` and is `'PM_DATA'` by default.
#
# ### Notebook outline
#
# **1. Data**
#
# **2. Train a first model**
#
# **3. Visualizing model outputs**
#
# **4. Train a more complex model: using log transformed target and namespace interactions**
#
#
# ### Requirements
#
# The notebook was developed for VW 8.11.0.
#
# It should work with older versions (>= 8.10) except for one cell in section 4.3 where the attribute `name` of `Feature` is accessed.
#
# ---
#
#
# ### Importing packages/modules

# +
import os
from os.path import join
import re

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import statsmodels

from vowpalwabbit.DFtoVW import DFtoVW
from vowpalwabbit.pyvw import vw


# -

# ### Defining functions/class
#

# +
# Data import/download functions
def download_data(dest_dir="PM_DATA"):
    import requests
    from io import BytesIO
    from rarfile import RarFile
    
    URL_PM_CITIES = "https://archive.ics.uci.edu/ml/machine-learning-databases/00394/FiveCitiePMData.rar"

    print(f"Downloading data at {URL_PM_CITIES}")
    r = requests.get(URL_PM_CITIES)
    bcontent = BytesIO(r.content)
    rf = RarFile(bcontent)
    
    print(f"Extracting content in folder {repr(dest_dir)}")
    rf.extractall(dest_dir)

def import_pm_city_data(folder_path, verbose=True):
    df_lst = []
    for fname in os.listdir(folder_path):
        fpath = join(folder_path, fname)
        if verbose:
            print(f"Importing file: {fpath}")
        city_name = re.sub("pm$", repl="", string=re.search("^[a-z]+", string=fname.lower()).group())
        
        df_city = pd.read_csv(fpath)
        df_city_clean = (df_city
                         .assign(city=city_name)
                         .drop(columns=["No"] + [col for col in df_city.columns if ("PM" in col) and (col != "PM_US Post")])
                         .rename(columns={"PM_US Post": "PM"})
                         .dropna(subset=["PM"])
        )
        df_lst += [df_city_clean]
    
    df_city = (
        pd.concat(df_lst) # append dataframes
        .sample(frac=1, random_state=123) # shuffle
        .reset_index(drop=True)
    )
    
    return df_city

# VW output parsing function/class
class VWLogParser:
    """Parser for Vowpal Wabbit output log/diagnostic"""
    
    def __init__(self, file_path_or_list):
        """The file name or list of lines to parse"""
        if isinstance(file_path_or_list, (list, str)):
            self.file_path_or_list = file_path_or_list
        else:
            raise TypeError("Argument `fname` should be a str (for file path) or a list of log lines")
    
    def parse(self):
        """Parse the output from `vw` command, return dataframe/dictionnaries with the associated data."""
        # Init containers 
        self.table_lst = []
        self.params = {}
        self.metrics = {}

        self.inside_table = False
        self.after_table = False
        
        if isinstance(self.file_path_or_list, list):
            for row in self.file_path_or_list:
                self._parse_vw_row(row)
        else:
            with open(self.file_path_or_list, "r") as f:
                for row in f:
                    self._parse_vw_row(row)
            

        self.df = self._make_output_df(self.table_lst)

        return self.params, self.df, self.metrics

    def _cast_string(self, s):
        """Cast to float or int if possible"""
        try:
            out = float(s)
        except ValueError:
            out = s
        else:
            if out.is_integer():
                out = int(out)

        return out

    def _make_output_df(self, lst):
        """Make dataframe from the list"""
        # Make columns from first and second elements of the list
        columns = [f"{first_row}_{second_row}" for (first_row, second_row) in zip(*lst[:2])]

        df = pd.DataFrame(data=lst[2:], columns=columns)

        # Cast cols to appropriate types
        int_cols = ["example_counter", "current_features"]
        for col in int_cols:
            df[col] = df[col].astype(int)

        float_cols = df.columns.drop(int_cols)
        for col in float_cols:
            df[col] = df[col].astype(float)

        return df

    def _parse_vw_row(self, row):
        """Parse row and add parsed elements to instance attributes params, metrics and table_lst"""
        if "=" in row:
            param_name, value = [element.strip() for element in row.split("=", maxsplit=1)]
            if self.after_table:
                self.metrics[param_name] = self._cast_string(value)
            else:
                self.params[param_name] = self._cast_string(value)
        elif ":" in row:
            param_name, value = [element.strip() for element in row.split(":", maxsplit=1)]
            self.params[param_name] = self._cast_string(value)

        elif not self.after_table:
            if re.match("average\s+since", row):
                self.inside_table = True
            if row == "\n":
                self.inside_table = False
                self.after_table = True
            if self.inside_table:
                self.table_lst += [row.split()]


# Model weight inspection functions                
def get_feature_names(df):
    cat_names = get_cat_feature_names(df)
    num_names = df.select_dtypes(np.number).columns.tolist()
    
    return cat_names + num_names


def get_cat_feature_names(df):
    unique_values_cat = df.select_dtypes(object).apply(lambda s: s.unique()).to_dict()
    cat_names = [f"{key}={value}" 
                 for (key, unique_values) in unique_values_cat.items()
                 for value in unique_values
    ]
    
    return cat_names


def get_weight_from_name(model, feature_name, namespace_name=" "):
    space_hash = model.hash_space(namespace_name)
    feat_hash = model.hash_feature(feature_name, space_hash)
    return model.get_weight(feat_hash)


# -

# ### Constants

# +
# Data
DATA_FOLDER = "PM_DATA"

# Graphical
SUPTITLE_FONTSIZE = 20
SUPTITLE_FONTWEIGHT = "bold"
TITLE_FONTSIZE = 15
# -

# ---

# # 1. Data
#
# ## 1.1. Import

# +
if not (os.path.isdir(DATA_FOLDER) and len(os.listdir(DATA_FOLDER)) == 5):
    download_data(DATA_FOLDER)

df_city = import_pm_city_data(DATA_FOLDER)
# -

# The full label of some features (the non obvious one):
#
# - PM: PM2.5 concentration (ug/m^3)
# - DEWP: Dew Point (Celsius Degree)
# - TEMP: Temperature (Celsius Degree)
# - HUMI: Humidity (%)
# - PRES: Pressure (hPa)
# - cbwd: Combined wind direction
# - Iws: Cumulated wind speed (m/s)
# - precipitation: hourly precipitation (mm)
# - Iprec: Cumulated precipitation (mm)

# ## 1.2. Types

# The types of the columns:

df_city.dtypes

# Some columns of the dataframe, such as `year`, `month`, `day`, `hour`, have been imported as integer type but **should be treated as categorical** by the model. Hence, we convert these columns to categorical (`str`) type.

# +
to_cat_cols = ["year", "month", "day", "hour", "season"]

for col in to_cat_cols:
    df_city[col] = df_city[col].astype(str)
# -

# Also, we standardize the numerical variables so we can compare their weights.

for col in df_city.select_dtypes(np.number).columns.difference(["PM", "log_PM"]):
    df_city[col] = (df_city[col] - df_city[col].mean()) / df_city[col].std()

# ## 1.3. Highly correlated features

# The most correlated (> 0.99) variable(s):

(df_city.corr()
 .reset_index().melt(id_vars="index")
 .loc[lambda df: df["index"] < df["variable"]] # since corr matrix is symetrical
 .loc[lambda df: (0.99 < df.value) & (df.value != 1)]
)

# We drop the `Iprec` variable as it is almost perfectly correlated with the `precipitation` variable.

df_city.drop(columns=["Iprec"], inplace=True)

# # 2. Train a first model
#
# ## 2.1. Converting pandas's DataFrame to Vowpal Wabbit input format.
#
# We now use `DFtoVW` class to convert this dataframe to a suitable VW input format. One nice property of the class is that it will build the appropriate VW type (numerical or categorical) based on the pandas's column type.
#
# There are 2 usages for the class `DFtoVW`:
# - **Basic usage**, using the `DFtoVW.from_colnames` convenience method.
# - **Advanced usage**, that relies on the VW input format specification (see: https://github.com/VowpalWabbit/vowpal_wabbit/wiki/Input-format) and is build upon classes such as `Feature`, `Namespace` or `SimpleLabel`, `MulticlassLabel` etc..
#
# The current section illustrates the basic usage. Section 4 will present the advanced usage.
#
# To use the `DFtoVW.from_colnames` method: the user provides the target name as a string (parameter `y`), the features as a list of string (`x`) and the dataframe (`df`) that has to be converted to the method. The user can also specify the type of label (argument `label_type`, default being `simple_label` for the regression). This call will build a `DFtoVW` object. 
#
# To get the actual examples (as a list), the user call the `convert_df` method.

# +
y = "PM"
x = df_city.columns.difference(["PM", "log_PM"], sort=False).tolist()

print("label:", y)
print("features:", x)

converter = DFtoVW.from_colnames(df=df_city, y=y, x=x)
examples = converter.convert_df()
# -

# We can inspect the first few examples:

examples[:5]

# For categorical features, we see that VW format is `feature_name=feature_value` whereas for numerical feature the format is `feature_name:feature_value`. 
#
# Also note that:
# - for categorical variables, VW adds `:1` behind the scene. For instance `day=14` is equivalent to `day=14:1`
# - The `=` doesn't have any special meaning and another symbol could have been used. However it's quite standard to use `=`
#
# Finally, if the feature name provide by the user is not found of the dataframe, the class will raise a `ValueError`:

try:
    DFtoVW.from_colnames(df=df_city, y=y, x=["TEMP", "COLUMN_NOT_IN_DF"])
except Exception as e:
    print(type(e))
    print(e)

# ## 2.2. Define and train model

# We now define the VW model. Note that we enable logging and also set the progress parameter (`-P`) to 1 to log the information _for each iteration_.

# +
model = vw(P=1, enable_logging=True)

for ex in examples:
    model.learn(ex)

model.finish()
# -

# # 3. Visualizing model outputs
#
# ## 3.1. Retrieving model's parameters, losses/predictions and summary metrics from the log
#
# Since we enable logging in the model definition (subsection 2.2), we can inspect the model's log. We use the class `VWLogParser` defined above.
#
# The log is returned as a list of strings by the `vw.get_log` method. Below is the first 20 lines:

model.get_log()[:20]

# The class `VWLogParser` parse following log information:
# - the parameters (in the first lines)
# - the information available at each iteration (in our case since `P=1`)
# - the summary metrics at the end of the log
#
# The parsed information is stored in objects (`dict` and dataframe) that can be easily manipulated.

log_parser = VWLogParser(model.get_log())
params, df_iter, summary_metrics = log_parser.parse()

# +
print("Model's parameter")
display(params)
print("\n")

print("Information at each iteration")
display(df_iter)
print("\n")

print("Summary metrics")
display(summary_metrics)
# -

# ## 3.2. Visualizing losses across iterations + distribution of selected metrics

scatter_var = ["average_loss", "since_last"]
g = sns.relplot(data=df_iter[scatter_var + ["example_counter"]].melt(id_vars="example_counter"), 
                x="example_counter", y="value", col="variable",
                col_wrap=3, facet_kws={'sharey': False, 'sharex': True}, kind="scatter", s=5, height=5, aspect=1.25)
g.fig.suptitle('Scatter plot of losses', fontsize=SUPTITLE_FONTSIZE, fontweight=SUPTITLE_FONTWEIGHT)
g.set_titles("{col_name}", size=TITLE_FONTSIZE)
g.fig.tight_layout()

# In the following of the section, we consider the metrics for iterations after 50 000 and above. Before this threshold, the average loss is too high.

start_idx = 50_000

# +
distr_vars = ["current_label", "current_predict", "current_features", "example_weight"]

g = sns.displot(data=df_iter.loc[start_idx:, distr_vars].melt(), 
                x="value", col="variable", multiple="dodge", hue="variable", 
                common_bins=False, facet_kws=dict(sharex=False, sharey=False), col_wrap=4, height=5)
g.fig.suptitle('Distribution of selected metrics', fontsize=SUPTITLE_FONTSIZE, fontweight=SUPTITLE_FONTWEIGHT)
g.set_titles("{col_name}", size=TITLE_FONTSIZE)
g.fig.tight_layout()
# -

# We notice that the distribution of the predictions differs substantially from the one of the label. 
#
# ## 3.3. Visualizing the predictions of the model
#
# This section offers a visualization of the model's predictions and compare them with the labels (the truth). 

# +
error = df_iter.current_label - df_iter.current_predict

f, (ax1, ax2) = plt.subplots(figsize=(20, 7), ncols=2)
f.suptitle('Predictions and errors', fontsize=SUPTITLE_FONTSIZE, fontweight=SUPTITLE_FONTWEIGHT)

# Scatterplot pred vs truth
sns.scatterplot(data=df_iter.loc[start_idx:], x="current_predict", y="current_label", ax=ax1, s=4, alpha=.2)
ax1.set_title("Prediction vs truth", fontsize=TITLE_FONTSIZE)

# Adding x=y line
range_x = range(0, int(df_iter[["current_predict", "current_label"]].max().min()))
ax1.plot(range_x, range_x, linestyle=":", color="orange", linewidth=2.5)

# Histogram of errors
sns.histplot(error[start_idx:], ax=ax2)
ax2.set_title("Distribution of errors", fontsize=TITLE_FONTSIZE)

plt.show()
# -

#
# The model tends to undervalue the concentration of PM. Another way to see it is that the errors (difference between the label and the prediction) doesn't seem normally distributed (especially the right tail).

# ## 3.4. Visualizing learnt weights

# +
# Get VW feature names
feature_names = get_feature_names(df_city)

# Get weights from feature names
weights_df = pd.DataFrame([(name, get_weight_from_name(model, name), "=" in name) for name in feature_names], 
                          columns=["vw_feature_name", "weight", "is_cat"])

# Adding columns for easier visualization
weights_df["feature_name"] = weights_df.apply(lambda row: row.vw_feature_name.split("=")[0] if row.is_cat else row.vw_feature_name, axis=1)
weights_df["feature_value"] = weights_df.apply(lambda row: row.vw_feature_name.split("=")[1].zfill(2) if row.is_cat else row.vw_feature_name, axis=1)
weights_df.sort_values(["feature_name", "feature_value"], inplace=True)

# +
g = sns.catplot(data=weights_df[lambda df: df.is_cat], kind="bar", x="feature_value", y="weight", col="feature_name", hue="feature_name",
                col_wrap=3, sharex=False, sharey=False, aspect=1.5, dodge=False)
g.fig.suptitle('Feature weights (categorical features)', fontsize=SUPTITLE_FONTSIZE, fontweight=SUPTITLE_FONTWEIGHT)
g.set_titles("{col_name}", size=TITLE_FONTSIZE)

# Add horizontal bar at y=0
for ax in g.axes.flat:
    ax.axhline(0, color="gray", linestyle=":")
g.fig.tight_layout()
# -

# Based on the weights learnt by the model, the pollution will be higher:
#
# - for wind orientation "calm and variable" (CV)
# - for cities such as Beijing or Chengdu
# - for winter season/months
# - for certain day of month such as 21-23 (oddly)
#
# The pollution will tend to be lower: 
#
# - for year 2015
# - for hours around noon

# +
f, ax = plt.subplots(figsize=(12, 5))
ax = sns.barplot(data=weights_df[lambda df: df.is_cat == False], x="feature_name", y="weight")
ax.set_title("Feature weight (numerical features)", fontsize=TITLE_FONTSIZE)

# Set xlabels in bold, remove x-axis title
ax.set_xticklabels(ax.get_xticklabels(), fontweight=SUPTITLE_FONTWEIGHT)
ax.set(xlabel=None)

# Add horizontal line
ax.axhline(0, linestyle=":", color="gray")

plt.show()
# -

# Higher cumulated wind speed (Iws) and higher temperature (TEMP) are associated with lower pollution concentration. 
#
# Higher air pressure (PRES) is associated with higher pollution concentration.

# # 4. Train a more complex model: using log transformed target and namespace interactions
#
# In this section we show the **advanced usage** of `DFtoVW` class. To do so, we will need to import some specific classes/types from the DFtoVW module.

from vowpalwabbit.DFtoVW import SimpleLabel, Namespace, Feature, DFtoVW

# The following drawing explains **how to use these classes**:
#
# <div>
# <img src="DFtoVW_usage.png" align="center" width="1250"/>
# </div>

# ## 4.1. Using log transformed variable
#
# The target distribution is close to a log-normal distribution:

# +
f, (ax1, ax2) = plt.subplots(figsize=(12, 5), ncols=2)
f.suptitle('Distribution of the target (PM)', fontsize=SUPTITLE_FONTSIZE, fontweight=SUPTITLE_FONTWEIGHT)

df_city["PM"].hist(bins=40, ax=ax1)
ax1.set_title("No transformation", fontsize=TITLE_FONTSIZE)

np.log(df_city["PM"]).hist(bins=40, ax=ax2)
ax2.set_title("Log-transformed", fontsize=TITLE_FONTSIZE)
# -

# We decide to train the model on the log transformation of the target (`log_PM`). 
#
# For regression purpose, we use `SimpleLabel` class (more details on this type on the [Input Format page](https://github.com/VowpalWabbit/vowpal_wabbit/wiki/Input-format#simple) of the wiki) to represent the target.

df_city["log_PM"] = df_city["PM"].apply(np.log)
label = SimpleLabel("log_PM")

# ## 4.2. Defining namespaces
#
# Namespaces are defined using the `Namespace` class. They have a name and are formed with a (list of) `Feature`.

# +
# Datetime namespace
datetime_features = ["year", "month", "day", "hour", "season"]
ns_datetime = Namespace(features=[Feature(col) for col in datetime_features], name="datetime_ns")

# City namespace
ns_city = Namespace(features=Feature("city"), name="city_ns")

# Other namespace
other_features = df_city.columns.difference(["PM", "log_PM", "city"] + datetime_features, sort=False).tolist()
ns_other = Namespace(features=[Feature(col) for col in other_features], name="other")
# -

# ## 4.3. Converting to VW format and training model

converter_advanced = DFtoVW(df=df_city, namespaces=[ns_datetime, ns_other, ns_city], label=label)

# `Namespaces` and `Feature` objects can be accessed using `DFtoVW`'s instance attributes `namespaces` and `features`:

for namespace in converter_advanced.namespaces:
    print("namespace:", namespace.name)
    for feature in namespace.features:
        print("\tfeature:", feature.name)

# We can inspect the first few examples:

examples_advanced = converter_advanced.convert_df()
examples_advanced[:5]

# In this new model we will add interactions between the "city" namespace and the "other" namespace (which contains atmosphetic features).

# +
model_advanced = vw(arg_str="--redefine C:=city_ns --redefine O:=other_ns --interactions CO", P=1, enable_logging=True)

for ex in examples_advanced:
    model_advanced.learn(ex)

model_advanced.finish()
# -

# ## 4.4. Visualizing model output
#
# We transformed the labels and predictions using exponential function (since the target is log-transformed) so we can compare them with the quantities of section 3.

# +
params_advanced, df_iter_advanced, metrics_advanced = VWLogParser(model_advanced.get_log()).parse()

df_iter_advanced.current_label = np.exp(df_iter_advanced.current_label)
df_iter_advanced.current_predict = np.exp(df_iter_advanced.current_predict)

# +
error_advanced = df_iter_advanced.current_label - df_iter_advanced.current_predict

f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(figsize=(20, 12), ncols=2, nrows=2)
f.suptitle('Predictions, truth and errors', fontsize=SUPTITLE_FONTSIZE, fontweight=SUPTITLE_FONTWEIGHT)

sns.histplot(df_iter_advanced.current_label.iloc[start_idx:], ax=ax1, color="blue")
ax1.set_title("Distribution of current_label", fontsize=TITLE_FONTSIZE)

sns.histplot(df_iter_advanced.current_predict.iloc[start_idx:], ax=ax2, color="orange")
ax2.set_title("Distribution of current_predict", fontsize=TITLE_FONTSIZE)

sns.scatterplot(data=df_iter_advanced.iloc[start_idx:], 
                x="current_predict", y="current_label", ax=ax3, s=4, alpha=.2)
ax3.set_title("Prediction vs truth", fontsize=TITLE_FONTSIZE)

# Add x=y line
range_x = range(0, int(df_iter_advanced.loc[start_idx:, ["current_predict", "current_label"]].max().min()))
ax3.plot(range_x, range_x, linestyle=":", color="orange", linewidth=2.5)

sns.histplot(error_advanced, ax=ax4)
ax4.set_title("Distribution of errors", fontsize=TITLE_FONTSIZE)

plt.show()
# -

# For this model, the distribution of the predictions is more in line with the distribution of the labels.
#
# The errors of this model are also closer to a normal distribution, even though the model still undervalues some observations.
#
# ## 4.5. Comparing models performance
#
# In this subsection the model in section 2 is called "simple" and the current model (with interactions and log-transformed) is called "advanced".

# +
loss_type = 1 # 1 for L1-loss, 2 for L2-loss
n_iter = 10_000 # Number of iterations on which to compute the loss (for moving average and final value)

df_loss = (pd.concat([error, error_advanced], axis=1, keys=["simple", "advanced"])
           .apply(lambda x: np.power(np.abs(x), loss_type)))

f, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 5.5))
f.suptitle("Comparison of models' performances", fontweight=SUPTITLE_FONTWEIGHT, fontsize=SUPTITLE_FONTSIZE)

df_loss.rolling(n_iter).mean().plot(ax=ax1)
ax1.set_title(f"Moving average of loss (over {n_iter} iterations)", fontsize=15)

loss_last_it = df_loss.tail(n_iter).mean()
pct_diff_loss = 100 * (loss_last_it.advanced / loss_last_it.simple - 1)

loss_last_it.plot.bar(ax=ax2, color="gray")
ax2.set_title(f"Loss computed on the last {n_iter} iterations", fontsize=15)
ax2.text(.75, loss_last_it.advanced * 1.025, f"Δ% loss= {pct_diff_loss:.2f}", fontsize=12, fontweight=SUPTITLE_FONTWEIGHT)
ax2.set_xticklabels(ax2.get_xticklabels(), fontsize=13, rotation=0)

plt.show()
# -

# The model with interactions and log-transformed target offers a substantial improvement over the simple model: the **loss** (computed on the last 10 000 observations) **decrease by about 19%**.

StarcoderdataJupyterScriptsDedupFiltered