<jupyter_start><jupyter_text>SetupLoad example data and prepare feature normalization.<jupyter_code>from __future__ import annotations

import logging
from typing import Dict, Tuple

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.datasets import load_boston
from sklearn.metrics import r2_score

from pr3.pursuit import PracticalProjectionPursuitRegressor
sns.set_style('whitegrid')
%matplotlib inline

RANDOM_SEED = 2021
TRG_RATIO = 0.75
np.random.seed(RANDOM_SEED)
boston = load_boston()
print(boston.DESCR)
xcols = boston.feature_names
ycol = "MEDV"

df = pd.DataFrame(
    data=boston.data,
    columns=xcols,
)
df[ycol] = boston.target
trg_idxs = np.random.binomial(1, p=TRG_RATIO, size=df.shape[0]).astype(bool)
trg_df = df.iloc[trg_idxs, :].copy()
tst_df = df.iloc[~trg_idxs, :].copy()
from dataclasses import dataclass

@dataclass
class FeatureNormalizer:
    logarithm: bool = False
    winsorize: bool = False
    zscore: bool = False
    
    _logarithm_cols: Dict[int, float] = None
    _winsorize_extremes: Dict[int, Tuple[float, float]] = None
    _zscore_stats: Dict[int, Tuple[float, float]] = None
    
    HEAVY_TAILED_SKEW: float = 2.0
    CONTINUOUS_UNIQUE_COUNT: int = 5
    LOG_SUMMAND_QUANTILE: float = 0.005
    EXTREME_QUANTILE: float = 0.005
    
    def fit(self, x: np.ndarray) -> FeatureNormalizer:
        x = x.copy()
        if self.logarithm:
            x = self._logarithm_fit(x)._logarithm_transform(x)
        if self.winsorize:
            x = self._winsorize_fit(x)._winsorize_transform(x)
        if self.zscore:
            x = self._zscore_fit(x)._zscore_transform(x)
            
        return self
        
    def transform(self, x: np.ndarray) -> np.ndarray:
        x = x.copy()
        if self.logarithm:
            x = self._logarithm_transform(x)
        if self.winsorize:
            x = self._winsorize_transform(x)
        if self.zscore:
            x = self._zscore_transform(x)
            
        return x
        
    def _logarithm_fit(self, x: np.ndarray) -> FeatureNormalizer:
        skews = ((x - x.mean(axis=0)) ** 3.0).mean(axis=0) / x.var(axis=0) ** 1.5
        self._logarithm_cols = {
            col: np.quantile(x[x[:, col] > 0, col], self.LOG_SUMMAND_QUANTILE)
            for col, skew in enumerate(skews)
            if skew > self.HEAVY_TAILED_SKEW
            and len(np.unique(x[:, col])) > self.CONTINUOUS_UNIQUE_COUNT
            and all(x[:, col] >= 0)
        }
        return self
        
    def _winsorize_fit(self, x: np.ndarray) -> FeatureNormalizer:
        lows = np.quantile(x, q=self.EXTREME_QUANTILE, axis=0)
        highs = np.quantile(x, q=1 - self.EXTREME_QUANTILE, axis=0)
        self._winsorize_extremes = dict(zip(range(x.shape[1]), zip(lows, highs)))
        return self
        
    def _zscore_fit(self, x: np.ndarray) -> FeatureNormalizer:
        mns = np.mean(x, axis=0)
        sds = np.std(x, axis=0)
        self._zscore_stats = dict(zip(range(x.shape[1]), zip(mns, sds)))
        return self
        
    def _logarithm_transform(self, x: np.ndarray) -> np.ndarray:
        if self._logarithm_cols is None:
            raise AttributeError("Log transform not yet fit on training data.")
        for col, quantile in self._logarithm_cols.items():
            x[:, col] = np.log(quantile + x[:, col])
        return x
    
    def _winsorize_transform(self, x: np.ndarray) -> np.ndarray:
        if self._winsorize_extremes is None:
            raise AttributeError("Winsorization transform not yet fit on training data.")
        for col, extremes in self._winsorize_extremes.items():
            x[:, col] = np.clip(x[:, col], extremes[0], extremes[1])
        return x
    
    def _zscore_transform(self, x: np.ndarray) -> np.ndarray:
        if self._zscore_stats is None:
            raise AttributeError("Z-score transform not yet fit on training data.")
        for col, stats in self._zscore_stats.items():
            x[:, col] = (x[:, col] - stats[0]) / stats[1]
        return x<jupyter_output><empty_output><jupyter_text>Model fitting We fit our projection pursuit regression below, where the key contributor to some loose form of "interpretability" is the sparsity constraint introduced by the least angle regression used for projection vector optimization. That is, by limiting the projection vector to have three nonzero coordinates (as specified by the argument `max_iter=3`), it becomes easier to to understand the meaning of each one-dimensional projection, and therefore also to understand the contribution of each ridge function.<jupyter_code>f = FeatureNormalizer(logarithm=True, winsorize=False, zscore=True)
trg_x = f.fit(trg_df[xcols].values).transform(trg_df[xcols].values)
trg_y = trg_df[ycol].values

ppr = PracticalProjectionPursuitRegressor(
    n_stages=5,
    learning_rate=1.0,
    ridge_function_class="polynomial",
    ridge_function_kwargs=dict(degree=3),
    projection_optimizer_class="least_angle",
    projection_optimizer_kwargs=dict(max_iter=3),
    random_state=RANDOM_SEED,
).fit(trg_x, trg_y)
ppr.plot_losses()
tst_df['yhat'] = ppr.predict(f.transform(tst_df[xcols].values))
print(f"Test R2: {r2_score(tst_df[ycol], tst_df['yhat']):0.3f}")<jupyter_output><empty_output><jupyter_text>Model visualization Below, we visualize the learned ridge functions (the nonlinear regression estimates in the one-dimensional projected space). Note that each stage fits against the residuals from previous stages, hence the learned functions do not appear to be good fits to projected data (except in the first stage). Furthermore, any apparent "gap" in the fit represents the component of variance explained by earlier stages of training.It can be very tempting to develop _post hoc_ "just so stories" upon viewing these plots; it may be safer to register any hypotheses about interpretation ahead of generating the plots below.<jupyter_code>ppr.plot(
    trg_x, 
    trg_y, 
    feature_names=xcols, 
    fig_height=2.5, 
    fig_width=5.0, 
    scatter_sample_ratio=0.5,
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10142552842864165, 'fraction_numerical': 0.00950352285761101, 'mean_word_length': 3.3198867657466384, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '103007', 'license': 'MIT', 'n_tokens_mistral': 2104, 'n_tokens_neox': 2055, 'n_words': 527, 'path': 'example.ipynb', 'repo_name': 'mikss/pr3', 'size': '6103'}