Unveiling trends in renewable energy with Linear Regression 🌍🔋

# Library Imports
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.feature_selection import SelectKBest, f_classif, RFE, VarianceThreshold
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression, Lasso
from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error, make_scorer
from sklearn.model_selection import cross_val_score
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.preprocessing import KBinsDiscretizer
import os
import numpy as np

training_path_training_set = os.path.join("data","Training_set_augmented.csv")
public_test_set = os.path.join("data","Public_Test_Set.csv")
target = 'Production (GWh)'

df = pd.read_csv(training_path_training_set)
test_set = pd.read_csv(public_test_set)
print(df.shape)

X = df.drop(columns=target)
y = df[target]

class LogTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, columns=None):
        self.columns = columns  # You can specify which columns to apply the log transformation to
    
    def fit(self, X, y=None):
        return self  # No fitting necessary
    
    def transform(self, X):
        # Check if input is a Pandas DataFrame or NumPy array
        if isinstance(X, np.ndarray):
            # If input is a NumPy array, we don't have columns, just indices
            X_copy = X.copy()

            # Handle NaN by replacing with column's median (for NumPy array)
            # If column-wise operations are needed, it's better to convert to DataFrame
            X_copy = np.where(np.isnan(X_copy), np.nanmedian(X_copy, axis=0), X_copy)

            # Handle infinity by replacing with the column's max value
            X_copy = np.where(np.isinf(X_copy), np.nanmax(X_copy, axis=0), X_copy)

            # Apply log transformation to all elements
            X_copy = np.log1p(X_copy)  # Log1p avoids log(0) by computing log(x + 1)

        elif isinstance(X, pd.DataFrame):
            # If input is a DataFrame, use pandas' functionality
            X_copy = X.copy()

            # Handle NaN by filling with the median of the column
            X_copy = X_copy.fillna(X_copy.median())

            # Handle infinity by replacing with the column's max value
            X_copy = X_copy.replace(np.inf, X_copy.max())

            # Apply log transformation to specified columns or all columns
            if self.columns is not None:
                for col in self.columns:
                    X_copy[col] = np.log1p(X_copy[col])  # Apply log1p to specified columns
            else:
                X_copy = np.log1p(X_copy)  # Apply log1p to all columns

        return X_copy

# Custom Binning Transformer to return DataFrame with column names
class BinningTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, bins=3):
        self.bins = bins
    
    def fit(self, X, y=None):
        return self  # No fitting necessary for binning
    
    def transform(self, X):
        # Ensure that the data is numeric
        X_copy = X.copy()  # Working on a copy of the data
        
        # Apply binning to each numeric column
        for i in range(X_copy.shape[1]):
            column_data = X_copy[:, i]  # Accessing column i as an array
            if np.issubdtype(column_data.dtype, np.number):  # Apply binning only to numeric columns
                discretizer = KBinsDiscretizer(n_bins=self.bins, encode='ordinal', strategy='uniform')
                binned = discretizer.fit_transform(column_data.reshape(-1, 1))  # Reshape for binning
                X_copy[:, i] = binned.flatten()  # Flatten to match original shape
        
        return X_copy  # Return the transformed NumPy array

# Custom Ratio Feature Transformer to return DataFrame with column names
class RatioFeatureTransformer(BaseEstimator, TransformerMixin):
    def fit(self, X, y=None):
        return self  # No fitting necessary
    
    def transform(self, X):
        X_copy = X.copy()
        # Create ratio features
        X_copy['GDP_to_Energy_Consumption'] = X_copy['GDP'] / X_copy['Energy Consumption']
        X_copy['Energy_Storage_to_Capacity'] = X_copy['Energy Storage Capacity'] / X_copy['Installed Capacity (MW)']
        
        # Ensure the result is a Pandas DataFrame with column names
        return X_copy

# Custom Time Feature Transformer to return DataFrame with column names
class TimeFeatureTransformer(BaseEstimator, TransformerMixin):
    def fit(self, X, y=None):
        return self  # No fitting necessary
    
    def transform(self, X):
        X_copy = X.copy()
        # Create time-based features
        X_copy['Year_Mod5'] = X_copy['Year'] % 5  # Capture cyclical trends based on year
        X_copy['Year_Diff'] = X_copy['Year'] - X_copy['Year'].min()  # Difference from the earliest year
        
        # Ensure the result is a Pandas DataFrame with column names
        return X_copy

numeric_features = X.select_dtypes(exclude='object').columns

log_transform_features = [
    'Biomass Availability',
    'CO2 Emissions',
    'Electricity Prices',
    'Energy Consumption',
    'Energy Storage Capacity',
    'Energy Subsidies',
    'GDP',
    'Government Policies',
    'Innovation Index',
    'Installed Capacity (MW)',
    'Investments (USD)',
    'Population',
    'Public-Private Partnerships in Energy',
    'R&D Expenditure',
    'Renewable Energy Patents',
    'Solar Irradiance',
    'Wind Speed'
]

passthrough_features = [
    'Control of Corruption',
    'Ease of Doing Business',
    'Energy Market Liberalization',
    'Grid Integration Capability',
    'Hydro Potential',
    'Innovation Index',
    'Political Stability',
    'Proportion of Energy from Renewables',
    'Regulatory Quality',
    'Renewable Energy Targets',
]


year_features = ['Year']
ratio_features = ['GDP', 'Energy Consumption', 'Energy Storage Capacity', 'Installed Capacity (MW)']
ignore_features = year_features + ratio_features

num_remaining_features = [col for col in numeric_features if col not in log_transform_features]
num_remaining_features = [col for col in num_remaining_features if col not in ratio_features]
num_remaining_features = [col for col in num_remaining_features if col not in passthrough_features]
num_remaining_features = [col for col in num_remaining_features if col not in year_features]


cat_features = X.select_dtypes(include='object').columns

### TRANSFORMERS
log_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('log', LogTransformer())
])


time_transformer = Pipeline(steps=[
    ('time', TimeFeatureTransformer()),
])

binning_ratio_transformer = Pipeline(steps=[
    ('ratio', RatioFeatureTransformer()),
    ('imputer', SimpleImputer(strategy='median')),  # Handle missing numerical values
    ('binning', BinningTransformer())
])

num_transformer = Pipeline([
    ('imputer', SimpleImputer(strategy='median')),  # Handle missing numerical values
    ('scaler', StandardScaler())  # Scale the remaining numerical features
])

cat_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='constant', fill_value='missing')),
    ('onehot', OneHotEncoder(handle_unknown='ignore'))
])

### PIPELINE
preprocessor = ColumnTransformer(
    transformers=[
        ('log', log_transformer, log_transform_features),
        ('passthrough', 'passthrough', passthrough_features),
        ('binning_ratio', binning_ratio_transformer, ratio_features),
        ('year', time_transformer, ['Year']),
        ('num', num_transformer, num_remaining_features),
        ('cat', cat_transformer, cat_features)
    ]
)

# Evaluating three different models: Linear Regression, Regularized Linear Regression, and Random Forest
for model in [LinearRegression(), Lasso(), RandomForestRegressor()]:
    pipeline = Pipeline(steps=[
        ('preprocessor', preprocessor),
        ('classifier', model)
    ])
    
    cv_scores = cross_val_score(
        pipeline,
        X,
        y,
        cv=5,
        scoring='neg_mean_squared_error',
        n_jobs=-1
    )
    rmse_scores = np.sqrt(-cv_scores)
    rmse = rmse_scores.mean()
    print(rmse)

# Elected to submit a model using Linear Regression

pipeline = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('classifier', LinearRegression())
])

pipeline.fit(X, y)

X_test = test_set.drop('ID', axis = 1)

y_pred = pipeline.predict(X_test)

prediction_df = pd.DataFrame({'ID': test_set['ID'], 'Predicted Production (GWh)': y_pred})
prediction_df.head()