visualize_model_apply.py

# -*- coding: utf-8 -*-
"""
Created on Wed Jan 29 19:08:59 2025

@author: luisfernando
"""

import numpy as np
import pandas as pd
from joblib import load
import matplotlib.pyplot as plt
from scipy.stats import qmc

import sys
import os
import time
import yaml

'''
Start the process below. The order of the functions is not chronological.
'''
START_PROCESS = time.time()

# Function 1: Plot the graphs of the estimated demand
def plot_results_with_diff(time_intervals_main, predictions_main, 
                           time_intervals_backup, predictions_backup, Diff,
                           archetype_label, use_diff=True):
    """
    Plots the predicted heating demand (both main and backup) alongside the temperature difference (Diff).

    Parameters:
        time_intervals_main (numpy array): Time intervals for main heating demand.
        predictions_main (numpy array): Predicted primary heating demand values.
        time_intervals_backup (numpy array): Time intervals for backup heating demand.
        predictions_backup (numpy array): Predicted backup heating demand values.
        Diff (numpy array): Temperature difference values.
        archetype_label (str): Label indicating the archetype being analyzed.
    """

    # Create figure and axes
    fig, ax1 = plt.subplots(figsize=(10, 5))

    # Plot main heating demand
    ax1.plot(time_intervals_main, predictions_main, label="Primary Heating Demand", color='b', linewidth=2)

    # Plot backup heating demand
    ax1.plot(time_intervals_backup, predictions_backup, label="Backup Heating Demand", color='g', linestyle='dashed', linewidth=2)

    # Set labels for primary y-axis
    ax1.set_xlabel("Time of Day (15-min period)")
    ax1.set_ylabel("Heating Demand (kW)", color='black')
    ax1.tick_params(axis='y', labelcolor='black')

    # Create a secondary y-axis for Diff
    ax2 = ax1.twinx()
    
    # Plot Temperature Difference (Diff)
    ax2.plot(time_intervals_main, Diff, label="Temperature Difference (Diff)", color='r', linestyle='--', linewidth=2)
    ax2.set_ylabel("Temperature Difference (°C)", color='r')
    ax2.tick_params(axis='y', labelcolor='r')

    # Add legends
    fig.legend(loc="upper left", bbox_to_anchor=(0.1, 0.9))

    # Grid and title
    ax1.grid(True, linestyle='--', alpha=0.7)
    plt.title(f"Heating Demand Predictions & Temperature Difference\nArchetype: {archetype_label}")

    # Generate a valid filename by replacing spaces and special characters
    save_dir = './plot_results/'
    os.makedirs(save_dir, exist_ok=True)
    safe_filename = archetype_label.replace(" ", "_").replace("/", "-").replace(":", "-") + ".png"
    save_path = os.path.join(save_dir, safe_filename)

    # Check if file already exists (to track overwrites)
    is_overwrite = os.path.exists(save_path)

    log_file="plot_results/plot_log.txt"
    # Log the save operation
    with open(log_file, "a") as log:
        if is_overwrite:
            log.write(f"OVERWRITE: {save_path}\n")
            print(f"⚠️ Overwritten: {save_path}")
        else:
            log.write(f"NEW FILE: {save_path}\n")

    # Save the figure
    plt.savefig(save_path, dpi=300, bbox_inches="tight")
    print(f"Saved plot: {save_path}")

    # Show plot
    plt.show()


# Function 2: to load scalers dynamically
def load_scalers(scaler_path):
    """Load scaler mean and scale values dynamically from a .joblib file."""
    scaler = load(scaler_path)
    return scaler.mean_, scaler.scale_


# Function 3: to generate input data for a single home, for testing
def generate_test_data(outdoor_temp):
    """Generate a simple test dataset with fixed input parameters."""
    # Generate outdoor temperature and Diff for a colder day
    indoor_temp = 22  # Indoor temperature in °C

    Diff = indoor_temp - outdoor_temp

    test_data = pd.DataFrame({
        'Diff': Diff,
        'sqft': [1800] * 96,
        'representative_income': [60000] * 96,
        'occupants': [3] * 96,
        'usage_level': [2.0] * 96,
        'rmb_heating_primary': [60.0] * 96,
        'rmb_heating_backup': [0.0] * 96,
        'vintage': [1995] * 96,
        'rmb_window_area': [200] * 96,
        'duct_leakage_and_insulation': [15.0] * 96,
        'heating_setpoint': [70.0] * 96,
        'hvac_heating_efficiency': [85.0] * 96
    })

    return test_data


# Function 4: to apply model, after being loaded.
def apply_model(model_path, scaler_path, test_data):
    """Apply the model to the test dataset after scaling."""
    # Load scalers
    scaler_mean, scaler_scale = load_scalers(scaler_path)

    # Ensure test_data is a proper DataFrame with float values
    test_data = test_data.astype(float)

    # Scale the data
    scaled_data = (test_data - scaler_mean) / scaler_scale

    # Convert to NumPy array to avoid feature name warning
    scaled_data_np = scaled_data.to_numpy()

    # Load model and predict
    model = load(model_path)
    predictions = model.predict(scaled_data_np)
    adjusted_predictions = predictions * 4

    return test_data.index, adjusted_predictions, scaler_mean, scaler_scale


# Function 5: to generate the archetypes using Latin Hypercube Sampling
def generate_archetype_centroids(n_archetypes, scaler_mean, scaler_scale, seed=42, STD_MULTIPLIER=2.0, main_or_back='main'):
    """
    Generates distinct archetype centroids using Latin Hypercube Sampling (LHS) and stores them in a dictionary.
    
    Parameters:
        n_archetypes (int): Number of distinct archetypes to generate.
        scaler_mean (list): Mean values of the scaler.
        scaler_scale (list): Scale (std deviation if StandardScaler) values of the scaler.
        seed (int): Random seed for reproducibility.
        STD_MULTIPLIER (float): Scaling factor for standard deviation (default = 2.0 for widely spaced archetypes).
        
    Returns:
        dict: A dictionary of archetype DataFrames.
    """
    np.random.seed(seed)
    num_features = len(scaler_mean) - 1  # Exclude Diff

    # Define centroid search space: ±2 * STD to ensure they are very distinct
    lower_bounds = np.array(scaler_mean[1:]) - (STD_MULTIPLIER * np.array(scaler_scale[1:]))
    upper_bounds = np.array(scaler_mean[1:]) + (STD_MULTIPLIER * np.array(scaler_scale[1:]))

    # Use Latin Hypercube Sampling to ensure structured diversity
    sampler = qmc.LatinHypercube(d=num_features, seed=seed)
    centroid_samples = sampler.random(n_archetypes)  # Generate structured LHS centroids
    centroids = qmc.scale(centroid_samples, lower_bounds, upper_bounds)  # Scale to bounds

    # Define feature names
    feature_names = [
        "sqft", "representative_income", "occupants", "usage_level",
        "rmb_heating_primary", "rmb_heating_backup", "vintage", "rmb_window_area",
        "duct_leakage_and_insulation", "heating_setpoint", "hvac_heating_efficiency"
    ]
    if main_or_back == 'back':
        feature_names = [
            'out.electricity.heating.energy_consumption'] + feature_names
    else:
        pass

    # Store centroids in a dictionary
    archetype_centroids = {}

    for i, centroid in enumerate(centroids):
        # Convert to DataFrame (1 row per archetype)
        df = pd.DataFrame([centroid], columns=feature_names)

        # Store in dictionary
        arch_name = f"Archetype_{i+1}"
        archetype_centroids[arch_name] = df

        # Print separated values for better readability
        print(f"\n===========================")
        print(f"Feature values of {arch_name}:")
        print("===========================")
        for feature, value in zip(feature_names, centroid):
            print(f"{feature}: {value:.2f}")

    return archetype_centroids


# Function 6: add the temperature difference to the archetypes
def apply_diff_to_archetype(arch_df, outdoor_temp, main_or_back='main'):
    """
    Applies a common outdoor temperature profile and computes Diff for an individual archetype.
    Also ensures that no feature (except Diff) has negative values.

    Parameters:
        arch_df (pd.DataFrame): DataFrame of an archetype (without Diff).
        outdoor_temp (numpy array): The common outdoor temperature time series.
        main_or_back (str): Indicates whether the predictor is "main" or "backup". 
                            "main" places Diff as the first column, "backup" places Diff as the second.

    Returns:
        pd.DataFrame: Updated DataFrame with Diff positioned accordingly.
    """
    # Convert heating_setpoint from Fahrenheit to Celsius
    heating_setpoint_C = (arch_df["heating_setpoint"] - 32) * (5 / 9)

    # Compute Diff for each home
    diff_values = []
    for i in range(len(arch_df)):
        # Apply the thermostat multiplier to the setpoint
        modified_setpoint = heating_setpoint_C.iloc[i] * thermostat_multiplier  # Element-wise multiplication

        diff_values.append(modified_setpoint - outdoor_temp)

    '''
    # BELOW WE PRINT THE THERMOSTAT CONTROLS TO ENSURE THE CHANGE
    print('1', heating_setpoint_C.iloc[i])
    print('2', thermostat_multiplier)
    print('3', modified_setpoint)
    print('4', outdoor_temp)
    print('5', diff_values)
    #'''

    # Convert Diff values into a properly formatted DataFrame
    diff_array = np.vstack(diff_values)  # Shape: (num_samples, time_intervals)

    # Store new dataframe with Diff
    df_with_diff = arch_df.copy()
    df_with_diff = df_with_diff.loc[df_with_diff.index.repeat(len(outdoor_temp))].reset_index(drop=True)  # Repeat rows
    df_with_diff["Diff"] = diff_array.flatten()  # Flatten the Diff array into a column

    # Clip all other columns (except Diff) to be non-negative
    for col in df_with_diff.columns:
        if col != "Diff":
            df_with_diff[col] = df_with_diff[col].clip(lower=0)

    # Define column order based on main_or_back argument
    remaining_columns = [col for col in df_with_diff.columns if col != "Diff"]

    if main_or_back == 'main':
        column_order = ["Diff"] + remaining_columns
    elif main_or_back == 'backup':
        column_order = [remaining_columns[0], "Diff"] + remaining_columns[1:] if remaining_columns else ["Diff"]

    df_with_diff = df_with_diff[column_order]

    return df_with_diff


# Function 7: add a column to the test_data
def update_arch_with_main_preds(test_data_back, predictions_main):
    """
    Adds primary heating predictions as the first column in the backup heating dataset.
    
    Parameters:
        test_data_back (pd.DataFrame): The dataset for backup heating.
        predictions_main (numpy array): Predicted main heating demand values.
    
    Returns:
        pd.DataFrame: Updated `test_data_back` with `predictions_main` as the first column.
    """
    # Ensure predictions_main is a DataFrame
    predictions_df = pd.DataFrame({
        'out.electricity.heating.energy_consumption': predictions_main})

    # Reset index to align both DataFrames properly
    test_data_back = test_data_back.reset_index(drop=True)
    predictions_df = predictions_df.reset_index(drop=True)

    # Concatenate predictions as the first column
    updated_test_data_back = pd.concat([predictions_df, test_data_back], axis=1)

    return updated_test_data_back


# Function 8: load the configuration file for the permuations
def load_config(yaml_path):
    """
    Reads the YAML configuration file.

    Parameters:
        yaml_path (str): Path to the YAML file.

    Returns:
        dict: Parsed YAML configuration as a Python dictionary.
    """
    with open(yaml_path, 'r') as file:
        config = yaml.safe_load(file)
    return config


# Function 9: aggregate the loads such that we get a 
def agg_perm_loads(all_perm_loads, config, outdoor_temp, pass_arch_label):
    """
    Aggregates heating loads across all permutations for each year.

    Parameters:
        all_perm_loads (dict): Nested dictionary containing load data for each permutation and year.
        config (dict): YAML configuration containing years.
        outdoor_temp (numpy array): Outdoor temperature time series.

    Returns:
        dict: Aggregated loads across permutations for each year.
    """
    # Initialize the aggregation structure
    aggregated_loads = {year: {"main": np.zeros(len(outdoor_temp)), "backup": np.zeros(len(outdoor_temp))}
                        for year in config["years"]}

    # Sum loads across all permutations
    for perm_name, yearly_data in all_perm_loads.items():
        for year, load_data in yearly_data.items():
            aggregated_loads[year]["main"] += load_data["main"]
            aggregated_loads[year]["backup"] += load_data["backup"]

    # Plot aggregated results
    for year in config["years"]:
        print(f"\nAggregated results for {year}...")

        plot_results_with_diff(
            np.arange(len(aggregated_loads[year]["main"])),  # Time index
            aggregated_loads[year]["main"], 
            np.arange(len(aggregated_loads[year]["backup"])),  # Time index
            aggregated_loads[year]["backup"], 
            outdoor_temp,  # Use outdoor temp as Diff reference
            archetype_label=pass_arch_label + f" ({year})",
            use_diff=False
        )

    return aggregated_loads


###############################################################################
# START THE PROGRAM HERE:

# 0.a) => Define two archetypical outdoor temperature profiles (Low & Mild)
low_temp_base = -20  # Cold winter day
mild_temp_base = -5  # Milder winter day
amplitude = 5  # Daily temperature variation
time_intervals = np.arange(0, 24, 0.25)  # 96 intervals per day

# Generate outdoor temperatures for two cases
outdoor_temp_low = low_temp_base + amplitude * (1 - np.cos(2 * np.pi * time_intervals / 24))
outdoor_temp_mild = mild_temp_base + amplitude * (1 - np.cos(2 * np.pi * time_intervals / 24))

# Store in a list for iteration
outdoor_temp_cases = [("LowTemp", outdoor_temp_low), ("MildTemp", outdoor_temp_mild)]

# 0.b) => Add the compoents for the thermostat modeling:
# Generate a baseline (no adjustment) and preheat multiplier
multiplier_baseline = np.ones(len(time_intervals))  # No change

multiplier_preheat = np.ones(len(time_intervals))  # Initialize
multiplier_preheat[(time_intervals >= 15) & (time_intervals <= 17)] = 1.1  # Raise setpoint for preheat
multiplier_preheat[(time_intervals >= 18) & (time_intervals <= 20)] = 0.9  # Lower setpoint after peak hours

# Store in a list for iteration
multiplier_cases = [("Baseline", multiplier_baseline), ("Preheat", multiplier_preheat)]

'''
General note: the results will have as many data points as there are temperatures
'''

# print('STAY UP UNTIL HERE')
# sys.exit()

# 1️) Load YAML configuration and initialize storing the results for printing
config_path = "viz_config.yaml"
with open(config_path, 'r') as file:
    config = yaml.safe_load(file)
all_results_print = []

# 2) Here we will iterate through the outdoor temperature and the multipliers:
# Iterate over thermostat multipliers (Baseline & Preheat)
for multiplier_label, thermostat_multiplier in multiplier_cases:
    print(f"\nProcessing Thermostat Strategy: {multiplier_label}...")

    # Iterate over outdoor temperature cases (Low & Mild)
    for temp_label, outdoor_temp in outdoor_temp_cases:
        print(f"\nProcessing Temperature Profile: {temp_label}...")

        # 3) Results storage, per outdoor temp and thermostart multiplier
        all_permutation_loads = {}

        # Iterate over permutations first
        for perm in config["permutations"]:
        
            perm_name = perm["name"]
            upgrade_id = perm["upgrade_id"]
            yearly_percentage = perm["yearly_percentage"]  # Get percentages
            print(f"\nProcessing {perm_name}...")

            # Initialize storage for this permutation
            all_permutation_loads[perm_name] = {}

            # 4) Generate archetypes for this permutation
            print(f"\nGenerating archetypes for {perm_name}...")
            start_time = time.time()
            archetype_centroids = generate_archetype_centroids(
                n_archetypes=3,

                scaler_mean=load_scalers(
                    perm["main_scaler_path"])[0],  # Get scaler mean from perm

                scaler_scale=load_scalers(
                    perm["main_scaler_path"])[1],  # Get scaler scale from perm

                STD_MULTIPLIER=1.0,

                main_or_back='main'
            )
            end_time = time.time()
            per_len = end_time - start_time
            print(f"\n -> ETf generate_archetype_centroids: {per_len:.4f} s")
            # ETf: "Execution time for..."

            # 5) Iterate over years for this permutation
            for year in config["years"]:
                # Get total households
                households = config["households_by_year"][year]

                print(f"\nProcessing {perm_name} for {year}...")
                print(f"({temp_label}, {multiplier_label})...")
                print(f"Households: {households}, Fraction: {fraction})...")

                # The totals take into account the number of houses and their distribution in the neighborhood.
                total_main = np.zeros(len(outdoor_temp))
                total_backup = np.zeros(len(outdoor_temp))

                # 6) Iterate over archetypes
                fraction_accum = 0
                households_accum = 0
                for arch_name, arch_df in archetype_centroids.items():
                    print(f"\nProcessing Archetype: {arch_name}...")

                    # Access the entire dictionary
                    archetype_fraction = config["archetype_shares"][arch_name]

                    # Get percentage of households using this permutation
                    fraction = yearly_percentage[year] * archetype_fraction

                    # Apply Diff and Filter Negatives
                    start_time = time.time()
                    archetype_with_diff = apply_diff_to_archetype(
                        arch_df, outdoor_temp, main_or_back='main'
                    )
                    end_time = time.time()
                    per_len = end_time - start_time
                    print(f"ETf apply_diff_to_archetype: {per_len:.4f} s")

                    # Apply Main Model
                    start_time = time.time()
                    time_intervals_main, predictions_main, _, _ = apply_model(
                        perm["main_model_path"], perm["main_scaler_path"], 
                        archetype_with_diff
                    )
                    end_time = time.time()
                    per_len = end_time - start_time
                    print(f"ETf apply_model (main): {per_len:.4f} s")

                    # Prepare Backup Model Input
                    archetype_with_diff_back = update_arch_with_main_preds(
                        archetype_with_diff, predictions_main
                    )

                    # Apply Backup Model if upgrade != 4
                    if upgrade_id != 4:
                        start_time = time.time()
                        time_intervals_back, predictions_back, _, _ = \
                            apply_model(perm["backup_model_path"], 
                                        perm["backup_scaler_path"], 
                                        archetype_with_diff_back
                        )
                        end_time = time.time()
                        per_len = end_time - start_time
                        print(f"ETf apply_model (backup): {per_len:.4f} s")
                    else:
                        predictions_back = np.zeros_like(predictions_main)
                        time_intervals_back = time_intervals_main.copy()

                    # 7) Plot Results (individual)
                    pass_archetype_label = \
                        f"Individual {perm_name} {temp_label} " + \
                        f"{multiplier_label} ({year})"
                    start_time = time.time()
                    plot_results_with_diff(
                        time_intervals_main, predictions_main, 
                        time_intervals_back, predictions_back, 
                        archetype_with_diff["Diff"], 
                        archetype_label=pass_archetype_label,
                        use_diff=True
                    )
                    end_time = time.time()
                    per_len = end_time - start_time
                    print(f"ETf plot_results_with_diff: {per_len:.4f} s")

                    # 8) Apply scaling factor: (Households * Fraction)
                    scaling_factor = households * fraction
                    fraction_accum += fraction
                    households_accum += scaling_factor
                    total_main += predictions_main * scaling_factor
                    total_backup += predictions_back * scaling_factor

                # Store results for this permutation and year
                all_permutation_loads[perm_name][year] = {
                    "main": total_main,
                    "backup": total_backup
                }

                # 🔟 Store results **PER TIME STEP** instead of total sum
                for i in range(len(time_intervals_main)):
                    all_results_print.append({
                        "temperature_scenario": temp_label,
                        "thermostat_scenario": multiplier_label,
                        "year": year,
                        "permutation": perm_name,
                        "upgrade_id": upgrade_id,
                        "num_houses": households_accum,
                        "percentage_adoption": fraction_accum,
                        "time_interval": time_intervals_main[i],
                        "main_demand_kW": float(total_main[i]),
                        "backup_demand_kW": float(total_backup[i])
                    })

                # 9) Plot Results (aggregate)
                pass_archetype_label = \
                    f"Aggregate {perm_name} {temp_label} " + \
                    f"{multiplier_label} ({year})"
                start_time = time.time()
                plot_results_with_diff(
                    time_intervals_main, total_main, 
                    time_intervals_back, total_backup, 
                    outdoor_temp, 
                    archetype_label=pass_archetype_label,
                    use_diff=False
                )
                end_time = time.time()
                per_len = end_time - start_time
                print(f"ETf plot_results_with_diff: {per_len:.4f} s")


        # 10) Aggregate all the results per permutation
        pass_arch_label = f"Aggregated Load {temp_label} {multiplier_label} "
        aggregated_loads = agg_perm_loads(
            all_permutation_loads, config, outdoor_temp, pass_arch_label
        )

# 11) Convert to DataFrame and Export
results_df = pd.DataFrame(all_results_print)
csv_filename = "aggregated_heating_demand_periodic.csv"
results_df.to_csv(csv_filename, index=False)

print(f"\n Periodic results saved to {csv_filename}")

END_PROCESS = time.time()
TIME_ELAPSED = -START_PROCESS + END_PROCESS
print('\n TIME ELAPSED:')
print(str(TIME_ELAPSED) + ' seconds /', str(TIME_ELAPSED/60) + ' minutes.')