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

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/afit-csce623-master/demos/blob/main/demo_ridge_lasso_elasticnet.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + id="51jGVasovC4n"
# This tutorial developed by <NAME>
# Apr 30, 2021

# Developed in support of education and research activities at 
# the Air Force Institute of Technology and the CSCE 623
# Machine Learning course

# For changes or questions, contact <EMAIL>

# Public Domain

# This is free and unencumbered software released into the public domain.

# Anyone is free to copy, modify, publish, use, compile, sell, or
# distribute this software, either in source code form or as a compiled
# binary, for any purpose, commercial or non-commercial, and by any
# means.

# In jurisdictions that recognize copyright laws, the author or authors
# of this software dedicate any and all copyright interest in the
# software to the public domain. We make this dedication for the benefit
# of the public at large and to the detriment of our heirs and
# successors. We intend this dedication to be an overt act of
# relinquishment in perpetuity of all present and future rights to this
# software under copyright law.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.

# For more information, please refer to <https://unlicense.org>

# + [markdown] id="RqRv7gnYXNHF"
# # Ridge, Lasso, and ElasticNet Demo

# + [markdown] id="JUcs5fH98-Lu"
# This notebook demonstrates the use of Ridge, Lasso, and ElasticNet on a linear regression dataset.

# + [markdown] id="Q8dnwK5nXSr2"
# ## Generate Data

# + [markdown] id="1ykyhqHJ-xYc"
# Here we generate regression data with 8 features, but only 3 features are needed to determine the target `y`.

# + id="O3XoZQtRM6hl"
import numpy as np
import pandas as pd
from sklearn.datasets import make_regression
import time

from IPython.display import Markdown as md
from IPython.display import display, Math, Latex

# initialize random state and regression parameters
random_state = np.random.RandomState(101)
n_features = 8
n_informative = 3
bias_data = random_state.uniform(0,10)
noise = 2

# Build a regression task
X, y, coef_data = make_regression(n_samples=6000, 
                                  n_features=n_features, n_informative=n_informative, 
                                  bias=bias_data, noise=noise,
                                  random_state=random_state, coef=True)


# + [markdown] id="ATRKQBA3XUkd"
# ## Helper Functions

# + [markdown] id="oQJREDcY-52E"
# The following helper functions will be used throughout this tutorial.

# + [markdown] id="cz4DvDr2Xhw8"
# ### Generate Model Function

# + [markdown] id="wffxOZIE-9SL"
# This function generates a model string of the form $y = \beta_0 + \beta_1 x_1 + \dots + \beta_n x_n$ where $n$ is the number of coefficients in `betas`.

# + id="uOZbGERCOnlE"
def generate_model_function(beta0, betas, feature_names):
    if np.isclose(beta0, 0):
        model_function = f'$y = '
    else:
        model_function = f'$y = {beta0:.3f} + '

    for beta, feature_name in zip(betas, feature_names):
        model_function += f'{beta:.3f}{feature_name} + '
    
    return model_function[:-3] + '$'


# + [markdown] id="_GIccjVpXlli"
# ### Clean Up Model

# + [markdown] id="taJNjfp7_bMA"
# Some models don't clean-up very small coefficients, leaving values looking like $ \dots + 0.000 x_3 + \dots $. This function will filter out coefficients that are near zero.

# + id="pBC7nnL4Qg6R"
def clean_up_model(coefs, names):
    non_zero_coef_filter = np.logical_not(np.isclose(coefs, 0))
    non_zero_coefs = coefs[non_zero_coef_filter]
    non_zero_names = names[non_zero_coef_filter]

    return non_zero_coefs, non_zero_names



# + [markdown] id="4j8f2or6Xp5T"
# ### Generate Clean Model

# + [markdown] id="BbcGUnDb_y7n"
# This function combines the above two functions to return a cleaned-up model string.

# + id="jNEoRvPhf1EL"
def generate_clean_model(beta0, coefs, names):
    coefs_, names_ = clean_up_model(coefs, names)
    return generate_model_function(beta0, coefs_, names_)


# + [markdown] id="Ato9lwsaXs-_"
# ### Powerset Function

# + [markdown] id="RarJVNi1_27H"
# This function generates the power set of a list, excluding the empty set.

# + id="dMnCU3E2rmQw"
# powerset function

from itertools import chain, combinations

def powerset(iterable, j=1, k=999):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)

    # start with 1, as we don't need the empty set
    p = list(chain.from_iterable(combinations(s, r) for r in range(1, len(s)+1)))
    return [item for item in p if len(item) >= j and len(item) <= k]

feature_powerset = powerset(range(n_features))


# + [markdown] id="ETAw-RhgbIqN"
# ### Time String

# + id="4qgmHLi_bL3C"
def convert_to_time(seconds):
    return time.strftime('%H:%M:%S', time.gmtime(seconds))


# + [markdown] id="00omhviges2K"
# ### Progress String

# + id="0SAwp5yoeu_I"
def print_progress_str(start_time, trial_time, total_trials, current_trial_idx, feature_set_name_str, best_score, current_score):

    time_elapsed_str = convert_to_time(time.perf_counter() - start_time)
    time_trial_sec = time.perf_counter() - trial_time

    status_str = f'{time_elapsed_str} {time_trial_sec:.2f}s {current_trial_idx} of {total_trials}: {feature_set_name_str} {current_score}'
    remaining_time = (time.perf_counter() - start_time) / current_trial_idx * (total_trials - current_trial_idx)
    
    if (current_score < best_score):
        print('\r', status_str)
    else:
        print('\r', status_str, f'Time Remaining: {time.strftime("%H:%M:%S", time.gmtime(remaining_time))}', end='')


# + [markdown] id="5tmXjkZyAAqn"
# ## Prepatory actons

# + [markdown] id="k6OBz3oiAIma"
# Here, we split the train and test datasets

# + id="WuphFhLhNZfd"
from sklearn.model_selection import train_test_split

# split and sequester test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state=random_state)

# + [markdown] id="HTjhNsAWAOMM"
# We set some variables that we will reuse throught the notebook.

# + id="JoEsu_RLNRgt"
kfold_size = 5

# create feature names -- these are generic feature names x_1, x_2, ..., x_25
feature_names = np.array([f'x_{{{idx+1}}}' for idx in range(n_features)])

# + [markdown] id="58IQ3r_9TExT"
# ## Ridge

# + [markdown] id="2mYuY_ayYxpo"
# ### Set alpha values

# + [markdown] id="aclnC-OAIum5"
# We will set alpha values. With Ridge regression, we will include 0.
#
#

# + id="QvJ5dcFJYzD3"
alpha_min = 0.0000001
alpha_max = 100000
alpha_num = int(abs(np.log10(alpha_min)) + abs(np.log10(alpha_max)) + 1)
alphas = np.geomspace(alpha_min, alpha_max, num=alpha_num)
alphas = np.insert(alphas, 0, 0)

# + [markdown] id="kCrdAxa0Yhry"
# ### Ridge - manual cross validation

# + [markdown] id="Qk2x6oTDI46y"
# The following code uses the Ridge regression model for each alpha value. We keep track of each RMSE, and select the model that yielded the lowest RMSE. We compare the model with the model that generated the original dataset.

# + colab={"base_uri": "https://localhost:8080/", "height": 98} id="ObNmaKPxWQBs" outputId="fe2a99a9-c000-4e1e-abe3-e09ce9462b3c"
from sklearn.linear_model import Ridge
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error

df_ridge = pd.DataFrame(columns=['alpha', 'rmse', 'r2'])

for alpha in alphas:
    ridge = Ridge(alpha=alpha, fit_intercept=True, random_state=random_state)
    scores = cross_val_score(ridge, X_train, y_train, scoring='neg_root_mean_squared_error', cv=kfold_size)
    ridge.fit(X_train, y_train)
    df_ridge.loc[len(df_ridge.index)] = [alpha, -np.mean(scores), ridge.score(X_test, y_test)]

best_ridge = {}
best_ridge['alpha'] = df_ridge.loc[np.argmin(df_ridge['rmse'])]['alpha']
best_ridge['model'] = Ridge(alpha=best_ridge['alpha'], fit_intercept=True, random_state=random_state)
best_ridge['model'].fit(X_train, y_train)
best_ridge['r2'] = best_ridge['model'].score(X_test, y_test)
y_pred = best_ridge['model'].predict(X_test)
best_ridge['rmse'] = mean_squared_error(y_test, y_pred, squared=False)

display(md(f'Best Ridge with $\\alpha$={best_ridge["alpha"]}: ' +
           f'{(generate_clean_model(best_ridge["model"].intercept_, best_ridge["model"].coef_, feature_names))}, ' +
           f'RMSE: {best_ridge["rmse"]:.4f}, $R^2$: {best_ridge["r2"]:.4f}'))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="vdQgNT4wYlDo"
# ### RidgeCV

# + [markdown] id="p0CclZFwJNKV"
# The following code uses the RidgeCV to automatically find the alpha value that yields the smallest RMSE. We compare the model with the model that generated the original dataset.

# + colab={"base_uri": "https://localhost:8080/", "height": 94} id="bWx9570fYm-E" outputId="06d5a895-9b24-4889-9c4a-e93b0661f294"
from sklearn.linear_model import RidgeCV

ridgeCV = RidgeCV(alphas=alphas,
                  scoring='neg_root_mean_squared_error',
                  cv=kfold_size)
ridgeCV.fit(X_train, y_train)
ridgeCV_r2 = ridgeCV.score(X_test, y_test)
y_pred = ridgeCV.predict(X_test)
ridgeCV_rmse = mean_squared_error(y_test, y_pred, squared=False)
display(md(f'RidgeCV with $\\alpha$={ridgeCV.alpha_}: ' +
           f'{(generate_clean_model(ridgeCV.intercept_, ridgeCV.coef_, feature_names))}, ' +
           f'RMSE: {ridgeCV_rmse:.4f}, $R^2$: {ridgeCV_r2:.4f}'))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="OGVfIW4ubks2"
# ### Ridge with subsets of features

# + [markdown] id="N75D307tJWaR"
# Neither Ridge nor RidgeCV is able to fully eliminate unimportant features. Here, we test the power set of features to find the set of features and alpha value that generates the lowest RMSE.

# + colab={"base_uri": "https://localhost:8080/", "height": 655} id="eiP13kfDbnF_" outputId="5538d710-aad6-4585-fb40-f041ae660313"
# use power set

df_ridge = pd.DataFrame(columns=['feature_set', 'alpha', 'rmse']) #, 'r2'])

trial_time = time.perf_counter()
start_time = time.perf_counter()
total_trials = len(feature_powerset) * len(alphas)
current_trial_idx = 0
best_rmse = 1000

for idx_feature, feature_set in enumerate(feature_powerset):

    feature_set_str = ','.join(str(e) for e in feature_set) # current feature set as string
    feature_set_name_str = ','.join(str(e) for e in feature_names[np.array(feature_set)]) # feature set names as string

    for idx_alpha, alpha in enumerate(alphas):
        current_trial_idx += 1

        ridge = Ridge(alpha=alpha, fit_intercept=True, random_state=random_state)
        scores = cross_val_score(ridge, X_train[:,np.array(feature_set)], y_train, scoring='neg_root_mean_squared_error', cv=kfold_size)

        # store details in dataframe for later accessing
        df_ridge.loc[len(df_ridge.index)] = [feature_set_str, alpha, -np.mean(scores)]

        # print to show progress
        print_progress_str(start_time, trial_time, total_trials, current_trial_idx, feature_set_name_str, best_rmse, -np.mean(scores))

        # update best_rmse
        best_rmse = -np.mean(scores) if -np.mean(scores) < best_rmse else best_rmse

        trial_time = time.perf_counter()

print('\r', f'Total time: {convert_to_time(time.perf_counter() - start_time)}')

best_ridge = {}
best_ridge['params'] = df_ridge.loc[np.argmin(df_ridge['rmse'])]
best_ridge['features'] = [int(e) for e in best_ridge['params']['feature_set'].split(',')]
best_ridge['feature_names'] = feature_names[np.array(best_ridge['features'])]
best_ridge['alpha'] = best_ridge['params']['alpha']

best_ridge['model'] = Ridge(alpha=best_ridge['alpha'])
best_ridge['model'].fit(X_train[:,np.array(best_ridge['features'])], y_train)
best_ridge['r2'] = best_ridge['model'].score(X_test[:,np.array(best_ridge['features'])], y_test)
y_pred = best_ridge['model'].predict(X_test[:,np.array(best_ridge['features'])])
best_ridge['rmse'] = mean_squared_error(y_test, y_pred, squared=False)

display(md(f'Best Ridge with $\\alpha$={best_ridge["alpha"]}: ' +
           f'{(generate_clean_model(best_ridge["model"].intercept_, best_ridge["model"].coef_, best_ridge["feature_names"]))}, ' + 
           f'RMSE: {best_ridge["rmse"]:.4f}, ' +
           f'$R^2$: {best_ridge["r2"]:.4f}'
           ))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="81tbIKecZErl"
# ## Lasso

# + [markdown] id="qihdwBP9ZHV_"
# ### Set alpha values

# + [markdown] id="JjTE1niELOI_"
# We will set alpha values. With Lasso regression, we will not include 0.

# + id="1fltmQvKZWBl"
alpha_min = 0.0000001
alpha_max = 100000
alpha_num = int(abs(np.log10(alpha_min)) + abs(np.log10(alpha_max)) + 1)
alphas = np.geomspace(alpha_min, alpha_max, num=alpha_num)

# Lasso is unreliable with alpha=0, use LinearRegression instead when alpha=0 is desired
#alphas = np.insert(alphas, 0, 0)

# + [markdown] id="LFgXQuwfZLlm"
# ### Lasso - manual cross validation

# + [markdown] id="C0gTkM7iLSl-"
# The following code uses the Lasso regression model for each alpha value. We keep track of each RMSE, and select the model that yielded the lowest RMSE. We compare the model with the model that generated the original dataset.

# + colab={"base_uri": "https://localhost:8080/", "height": 77} id="jka9W-bxZf-V" outputId="4e2384a7-5397-4275-958e-e911341238d1"
from sklearn.linear_model import Lasso
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error

df_lasso = pd.DataFrame(columns=['alpha', 'rmse', 'r2'])

for alpha in alphas:
    lasso = Lasso(alpha=alpha, fit_intercept=True, random_state=random_state)
    scores = cross_val_score(lasso, X_train, y_train, scoring='neg_root_mean_squared_error', cv=kfold_size)
    lasso.fit(X_train, y_train)
    df_lasso.loc[len(df_lasso.index)] = [alpha, -np.mean(scores), lasso.score(X_test, y_test)]

best_lasso = {}
best_lasso['alpha'] = df_lasso.loc[np.argmin(df_lasso['rmse'])]['alpha']
best_lasso['model'] = Lasso(alpha=best_lasso['alpha'], fit_intercept=True, random_state=random_state)
best_lasso['model'].fit(X_train, y_train)
best_lasso['r2'] = best_lasso['model'].score(X_test, y_test)
y_pred = best_lasso['model'].predict(X_test)
best_lasso['rmse'] = mean_squared_error(y_test, y_pred, squared=False)

display(md(f'Best Ridge with $\\alpha$={best_lasso["alpha"]}: ' +
           f'{(generate_clean_model(best_lasso["model"].intercept_, best_lasso["model"].coef_, feature_names))}, ' +
           f'RMSE: {best_lasso["rmse"]:.4f}, $R^2$: {best_lasso["r2"]:.4f}'))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="UwVSU7dEZOhF"
# ### LassoCV

# + [markdown] id="RMVm5z1SLYwY"
# The following code uses the LassoCV to automatically find the alpha value that yields the smallest RMSE. We compare the model with the model that generated the original dataset.

# + colab={"base_uri": "https://localhost:8080/", "height": 77} id="-UP8k_u3ZGWE" outputId="dce772a7-e4c7-4599-e3f1-2ffc7d3d4595"
from sklearn.linear_model import LassoCV

lassoCV = LassoCV(alphas=alphas,
                #   scoring='neg_root_mean_squared_error',
                  cv=kfold_size)
lassoCV.fit(X_train, y_train)
lassoCV_r2 = lassoCV.score(X_test, y_test)
y_pred = lassoCV.predict(X_test)
lassoCV_rmse = mean_squared_error(y_test, y_pred, squared=False)
display(md(f'LassoCV with $\\alpha$={lassoCV.alpha_}: ' +
           f'{(generate_clean_model(lassoCV.intercept_, lassoCV.coef_, feature_names))}, ' +
           f'RMSE: {lassoCV_rmse:.4f}, $R^2$: {lassoCV_r2:.4f}'))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="ZobiGFXhcIVi"
# ### Lasso with subsets of features

# + [markdown] id="qnkfxZvPLnGd"
# Though it performs better than Ridge, neither Lasso nor LassoCV is able to fully eliminate unimportant features. Here, we test the power set of features to find the set of features and alpha value that generates the lowest RMSE.

# + colab={"base_uri": "https://localhost:8080/", "height": 417} id="RG91g4QsDrdl" outputId="5426c18e-18b2-4b96-f119-6d2e34cc2bdf"
# use power set

df_lasso = pd.DataFrame(columns=['feature_set', 'alpha', 'rmse'])

trial_time = time.perf_counter()
start_time = time.perf_counter()
total_trials = len(feature_powerset) * len(alphas)
current_trial_idx = 0
best_rmse = 1000

for feature_set in feature_powerset:

    feature_set_str = ','.join(str(e) for e in feature_set) # current feature set as string
    feature_set_name_str = ','.join(str(e) for e in feature_names[np.array(feature_set)]) # feature set names as string

    for alpha in alphas:
        current_trial_idx += 1

        lasso = Lasso(alpha=alpha, fit_intercept=True, random_state=random_state)
        scores = cross_val_score(lasso, X_train[:,np.array(feature_set)], y_train, scoring='neg_root_mean_squared_error', cv=kfold_size)

        # store details in dataframe for later accessing
        df_lasso.loc[len(df_lasso.index)] = [feature_set_str, alpha, -np.mean(scores)]

        # print to show progress
        print_progress_str(start_time, trial_time, total_trials, current_trial_idx, feature_set_name_str, best_rmse, -np.mean(scores))

        # update best_rmse
        best_rmse = -np.mean(scores) if -np.mean(scores) < best_rmse else best_rmse

        trial_time = time.perf_counter()

print('\r', f'Total time: {convert_to_time(time.perf_counter() - start_time)}')

best_lasso = {}
best_lasso['params'] = df_lasso.loc[np.argmin(df_lasso['rmse'])]
best_lasso['features'] = [int(e) for e in best_lasso['params']['feature_set'].split(',')]
best_lasso['feature_names'] = feature_names[np.array(best_lasso['features'])]
best_lasso['alpha'] = best_lasso['params']['alpha']

best_lasso['model'] = Lasso(alpha=best_lasso['alpha'])
best_lasso['model'].fit(X_train[:,np.array(best_lasso['features'])], y_train)
best_lasso['r2'] = best_lasso['model'].score(X_test[:,np.array(best_lasso['features'])], y_test)
y_pred = best_lasso['model'].predict(X_test[:,np.array(best_lasso['features'])])
best_lasso['rmse'] = mean_squared_error(y_test, y_pred, squared=False)

display(md(f'Best Lasso with $\\alpha$={best_lasso["alpha"]}: ' +
           f'{(generate_clean_model(best_lasso["model"].intercept_, best_lasso["model"].coef_, best_lasso["feature_names"]))}, ' + 
           f'RMSE: {best_lasso["rmse"]:.4f}, ' +
           f'$R^2$: {best_lasso["r2"]:.4f}'
           ))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="U7zE5eoVeDzX"
# ## ElasticNet

# + [markdown] id="1b6kgVtaeaLT"
# ### Set alpha values

# + [markdown] id="qNeSzHUTLw4I"
# We will set alpha values. With ElasticNet regression, which includes Lasso when l1_ratio=1, we will not include 0.

# + id="3nUnghLXeb7e"
alpha_min = 0.0000001
alpha_max = 100000
alpha_num = int(abs(np.log10(alpha_min)) + abs(np.log10(alpha_max)) + 1)
alphas = np.geomspace(alpha_min, alpha_max, num=alpha_num)

# Lasso is unreliable with alpha=0, use LinearRegression instead when alpha=0 is desired
#alphas = np.insert(alphas, 0, 0)

# + [markdown] id="BU6kNhXwekAb"
# ### Set lr_ratio values

# + [markdown] id="vuDJQsPhME9R"
# Set l1_ratio values. This is that amount of l1 regularization the linear model will attempt to perform (l1_ratio = 1.0 corresponds with Lasso, and l1_ratio = 0.0 corresponds with Ridge).

# + id="XyT63444eGle"
l1_ratios = np.linspace(0.01, 1.0 ,num=100)

# + [markdown] id="4ztjUcEvjwTx"
# ### ElasticNetCV

# + [markdown] id="VlK4GqaGMDh2"
# The following code uses the ElasticNetCV to automatically find the alpha and l1_ratio values that yield the smallest RMSE. We compare the model with the model that generated the original dataset.

# + colab={"base_uri": "https://localhost:8080/", "height": 77} id="K7nDCe_Vf29r" outputId="4bcabe94-c32a-477b-f0c5-51a4f2ebabf8"
from sklearn.linear_model import ElasticNetCV

elastic_netCV = ElasticNetCV(alphas=alphas,
                      l1_ratio = l1_ratios,
                #   scoring='neg_root_mean_squared_error',
                  cv=kfold_size,
                  max_iter=1000)
elastic_netCV.fit(X_train, y_train)
elastic_netCV_r2 = elastic_netCV.score(X_test, y_test)
y_pred = elastic_netCV.predict(X_test)
elastic_netCV_rmse = mean_squared_error(y_test, y_pred, squared=False)
display(md(f'ElasticCV with $\\alpha$={elastic_netCV.alpha_}, $l1$ ratio={elastic_netCV.l1_ratio_}: ' +
           f'{(generate_clean_model(elastic_netCV.intercept_, elastic_netCV.coef_, feature_names))}, ' +
           f'RMSE: {elastic_netCV_rmse:.4f}, $R^2$: {elastic_netCV_r2:.4f}'))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

# + [markdown] id="yQ_mQ0KvjzUf"
# ### ElasticNet with subsets of features

# + [markdown] id="Qir406pFMsQ5"
# ElasticNetCV's performance in this case is nearly identical to that of Lasso.Here, we test the power set of features to find the set of features, alpha, and l1_ratio values that generate the lowest RMSE.

# + id="x4FMmp1Sj25p" colab={"base_uri": "https://localhost:8080/", "height": 434} outputId="b688e83e-0ba0-4eed-b9fc-a50f8f9d8cbb"
# use power set

from sklearn.linear_model import ElasticNet

df_elastic_net = pd.DataFrame(columns=['feature_set', 'alpha', 'rmse'])

trial_time = time.perf_counter()
start_time = time.perf_counter()
total_trials = len(feature_powerset) * len(alphas) * len(l1_ratios)
current_trial_idx = 0
best_rmse = 1000

for feature_set in feature_powerset:

    feature_set_str = ','.join(str(e) for e in feature_set) # current feature set as string
    feature_set_name_str = ','.join(str(e) for e in feature_names[np.array(feature_set)]) # feature set names as string

    for alpha in alphas:
        for l1_ratio in l1_ratios:
            current_trial_idx += 1

            # create model
            elastic_net = ElasticNet(alpha=alpha, fit_intercept=True, random_state=random_state)

            # calculate rmse scores
            scores = cross_val_score(elastic_net, X_train[:,np.array(feature_set)], y_train, scoring='neg_root_mean_squared_error', cv=kfold_size)

            # store details in dataframe for later accessing
            df_elastic_net.loc[len(df_elastic_net.index)] = [feature_set_str, alpha, -np.mean(scores)]

            # print to show progress
            print_progress_str(start_time, trial_time, total_trials, current_trial_idx, feature_set_name_str, best_rmse, -np.mean(scores))

            # update best_rmse
            best_rmse = -np.mean(scores) if -np.mean(scores) < best_rmse else best_rmse

            trial_time = time.perf_counter()

print('\r', f'Total time: {convert_to_time(time.perf_counter() - start_time)}')

best_elastic_net = {}
best_elastic_net['params'] = df_elastic_net.loc[np.argmin(df_elastic_net['rmse'])]
best_elastic_net['features'] = [int(e) for e in best_elastic_net['params']['feature_set'].split(',')]
best_elastic_net['feature_names'] = feature_names[np.array(best_elastic_net['features'])]
best_elastic_net['alpha'] = best_elastic_net['params']['alpha']

best_elastic_net['model'] = ElasticNet(alpha=best_elastic_net['alpha'])
best_elastic_net['model'].fit(X_train[:,np.array(best_elastic_net['features'])], y_train)
best_elastic_net['r2'] = best_elastic_net['model'].score(X_test[:,np.array(best_elastic_net['features'])], y_test)
y_pred = best_elastic_net['model'].predict(X_test[:,np.array(best_elastic_net['features'])])
best_elastic_net['rmse'] = mean_squared_error(y_test, y_pred, squared=False)

display(md(f'Best ElasticNet with $\\alpha$={best_elastic_net["alpha"]}, $l1$ ratio={elastic_netCV.l1_ratio_}: ' +
           f'{(generate_clean_model(best_elastic_net["model"].intercept_, best_elastic_net["model"].coef_, best_elastic_net["feature_names"]))}, ' + 
           f'RMSE: {best_elastic_net["rmse"]:.4f}, ' +
           f'$R^2$: {best_elastic_net["r2"]:.4f}'
           ))
display(md(f'Original Data: {(generate_clean_model(bias_data, coef_data, feature_names))}'))

StarcoderdataJupyterScriptsDedupFiltered