<jupyter_start><jupyter_text>This notebook contains the codes to extract and process the necessary simulation data for constructing the machine learning model. The contents must be incorporated to the resilience_metrics.py once finalized.<jupyter_code>%load_ext autoreload
%autoreload 2

import warnings
warnings.filterwarnings('ignore')

from IPython.display import clear_output

import os
import pandas as pd
from pathlib import Path
from sklearn import metrics
import statistics
import copy

import infrarisk.src.network_sim_models.interdependencies as interdependencies
from infrarisk.src.network_sim_models.integrated_network import *
network_dir = Path('../../data/networks/micropolis')
water_folder = network_dir/'water/high'
power_folder = network_dir/'power/low'

micropolis_network = IntegratedNetwork(name = 'Micropolis', 
                                       water_folder= water_folder,
                                       power_folder = power_folder,
                                       water_sim_type = 'PDA',
                                       power_sim_type='1ph')
# Set scenarios folder
folder = Path('../../data/networks/micropolis/scenarios')
scenarios = [f for f in sorted(os.listdir(folder))]

#list of recovery strategies to be considered
strategies = ['capacity', 'centrality', 'zone']

#create the empty dataframe for ML dataset
ml_df = pd.DataFrame(columns =["scenario",
         "strategy",
         "meshedness",
         "water_crew",
         "power_crew",
         "transpo_crew",
         'water_perf_ecs',
         'water_perf_pcs',
         'power_perf_ecs',
         'power_perf_pcs',
         'water_mains',
         "power_lines", 
         "transpo_links",
         "all_compons"])

abnormal_results = []

junc_list = micropolis_network.wn.junction_name_list
base_water_demands = micropolis_network.base_water_node_supply
base_power_demands = micropolis_network.base_power_supply

for index, scenario in enumerate(scenarios):
    print(index, ". ", scenario)
    meshed_levels = os.listdir(f"{folder}/{scenario}")
    for mesh_level in meshed_levels:
        print(mesh_level)
        ml_df_new = {"scenario": scenario,
            "strategy": None,
            "meshedness" : mesh_level,
            "water_crew": 0,
            "power_crew": 0,
            "transpo_crew": 0,
            'water_perf_ecs': None,
            'water_perf_pcs': None,
            'power_perf_ecs': None ,
            'power_perf_pcs': None,
            "water_mains": 0,
            "power_lines": 0, 
            "transpo_links": 0, 
            "all_compons": 0}
        
        disruption_file = pd.read_csv(f"{folder}/{scenario}/{mesh_level}/disruption_file.csv")
        ml_df_new["all_compons"] = disruption_file.shape[0]
        for _, row in disruption_file.iterrows():
            if row['components'].startswith('P_L'):
                ml_df_new['power_lines'] += 1
            elif row['components'].startswith('W_PMA'):
                ml_df_new['water_mains'] += 1
            elif row['components'].startswith('T_L'):
                ml_df_new['transpo_links'] += 1
            else:
                print("Component type not detectable.")            
        
        for strategy in strategies: 
            try: 
                ml_df_new['strategy'] = strategy
                water_demands_file = f"{folder}/{scenario}/{mesh_level}/{strategy}/water_junc_demand.csv"
                power_demands_file = f"{folder}/{scenario}/{mesh_level}/{strategy}/power_load_demand.csv"
                water_pressure_file = f"{folder}/{scenario}/{mesh_level}/{strategy}/water_node_pressure.csv"
                crew_size_file = f"{folder}/{scenario}/{mesh_level}/{strategy}/crew_size.csv"
                
                if os.path.isfile(water_demands_file):
                    water_demands = pd.read_csv(water_demands_file, sep = "\t")
                    water_time_list = water_demands.time/60
                    water_time_list = water_time_list.tolist()
                    rel_time_list = water_demands['time'] % (24*3600)
                    index_list = [int(x/60) for x in rel_time_list if np.isnan(x) == False]
                    water_demands = water_demands[junc_list]
                    
                    # water_pressures = pd.read_csv(water_pressure_file, sep = "\t")
                    # water_pressures = water_pressures[junc_list]
                    # water_press_corrections = copy.deepcopy(water_pressures)
                    
                    power_demands = pd.read_csv(power_demands_file, sep = "\t")
                    power_time_list = power_demands.time/60
                    power_time_list= power_time_list.tolist()
                    
                    base_water_demands_new = base_water_demands.iloc[index_list].reset_index(drop=True)
                    base_water_demands_new = base_water_demands_new[junc_list]
                    
                    water_demands_ratio = water_demands/ base_water_demands_new
                    water_demands_ratio = water_demands_ratio.clip(upper=1, lower=0)
                    
                    water_ecs_list = water_demands_ratio.mean(axis = 1, skipna = True).tolist()
                    
                    if water_ecs_list[-1] < 0.9:
                        abnormal_results.append(scenario)
                        
                    water_pcs_list = pd.concat([water_demands, base_water_demands_new]).min(level=0).sum(axis=1, skipna = True)/base_water_demands_new.sum(axis=1, skipna = True)
                    water_pcs_list = water_pcs_list.tolist()

                    base_load_demands = pd.DataFrame(base_power_demands.load.p_mw.tolist() + base_power_demands.motor.pn_mech_mw.tolist()).transpose()
                    base_load_demands.columns = base_power_demands.load.name.tolist() + base_power_demands.motor.name.tolist()
                    base_load_demands = pd.concat([base_load_demands]*(power_demands.shape[0])).reset_index(drop=True)

                    power_demand_ratio = power_demands.iloc[:,1:] / base_load_demands
                    power_demand_ratio = power_demand_ratio.clip(upper=1)

                    power_ecs_list = power_demand_ratio.mean(axis = 1, skipna = True).tolist()
                    power_pcs_list = pd.concat([power_demands.iloc[:,1:], base_load_demands]).min(level=0).sum(axis=1, skipna = True)/base_load_demands.sum(axis=1, skipna = True)
                    power_pcs_list = power_pcs_list.tolist()
                    
                    crew_size = pd.read_csv(crew_size_file, sep = ",")
                    ml_df_new['power_crew'] = crew_size.iloc[0, 1]
                    ml_df_new['water_crew'] = crew_size.iloc[1, 1]
                    ml_df_new['transpo_crew'] = crew_size.iloc[2, 1]
                    
                    ml_df_new['water_perf_ecs'] = round(metrics.auc(water_time_list, [1 - x for x in water_ecs_list]), 3)
                    ml_df_new['water_perf_pcs'] = round(metrics.auc(water_time_list, [1 - x for x in water_pcs_list]), 3)
                    ml_df_new['power_perf_ecs'] = round(metrics.auc(power_time_list, [1 - x for x in power_ecs_list]), 3)
                    ml_df_new['power_perf_pcs'] = round(metrics.auc(power_time_list, [1 - x for x in power_pcs_list]), 3)
                    
                    power_auc_df = pd.DataFrame(data = {'time': power_time_list, 
                                                        'ecs': power_ecs_list, 
                                                        'pcs': power_pcs_list})
                    power_auc_df.to_csv(f"{folder}/{scenario}/{mesh_level}/{strategy}/power_auc.csv", index = False)
                    
                    water_auc_df = pd.DataFrame(data = {'time': water_time_list, 
                                                        'ecs': water_ecs_list, 
                                                        'pcs': water_pcs_list})
                    water_auc_df.to_csv(f"{folder}/{scenario}/{mesh_level}/{strategy}/water_auc.csv", index = False)
                    
                    ml_df = ml_df.append(ml_df_new, ignore_index=True)
                    print(ml_df.iloc[-1,:].tolist())
                    
                else:
                    pass
            except AttributeError:
                pass
    clear_output(wait=True)
    
abnormal_results = list(set(abnormal_results))
ml_df.drop_duplicates()
ml_df.shape
abnormal_results
ml_df.to_csv("auc_df.csv", index = False)
ml_df = pd.read_csv("auc_df.csv")
ml_df.head()
ml_df['ecs_auc'] = ml_df['water_perf_ecs']*0.5 + ml_df['power_perf_ecs']*0.5
ml_df['pcs_auc'] = ml_df['water_perf_pcs']*0.5 + ml_df['power_perf_pcs']*0.5
ml_df.head()
ml_df.columns

import numpy as np
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.feature_selection import SelectFromModel
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn import tree
from sklearn.linear_model import LinearRegression<jupyter_output><empty_output><jupyter_text>Decision Tree Model<jupyter_code>lab = 'ecs_auc'

features = ml_df[['strategy', 'meshedness', 
                  'water_crew', 'power_crew', 'transpo_crew', 
                  'water_mains' , 'power_lines', 'transpo_links', lab]]
features = features.replace('nan', np.NaN)
features = features.dropna()

labels = features[lab]
del features[lab]
features = pd.get_dummies(features)

feature_list = list(features.columns)
features = np.array(features)

clf = ExtraTreesRegressor(n_estimators=500)
clf = clf.fit(features, labels)

print(feature_list)
print(clf.feature_importances_)

pd.DataFrame(data = {'features': feature_list,
                     'feature importance': clf.feature_importances_})
model = SelectFromModel(clf, prefit=True, threshold = '0.5*median')
features_new = model.transform(features)
features_new.shape
feature_list_new = []

for i, feature in enumerate(feature_list):
    if model.get_support()[i]:
        feature_list_new.append(feature)

feature_list_new
train_features, test_features, train_labels, test_labels = train_test_split(features_new, labels, test_size = 0.25, random_state = 43)
print('Training Features Shape:', train_features.shape)
print('Training Labels Shape:', train_labels.shape)
print('Testing Features Shape:', test_features.shape)
print('Testing Labels Shape:', test_labels.shape)
parameters = {'max_depth':range(2,10), 'min_samples_leaf': range(3,10)}
clf = GridSearchCV(tree.DecisionTreeRegressor(), parameters, n_jobs=5, scoring = 'accuracy')
clf.fit(X=train_features, y=train_labels)
tree_model = clf.best_estimator_
print (clf.best_score_, clf.best_params_) 
clf = tree.DecisionTreeRegressor(max_depth = 2, random_state = 0, min_samples_leaf=3)
# Train the model on training data
clf.fit(train_features, train_labels)
y_pred=clf.predict(train_features)   

def measure_performance(y_test, y_pred, X_train):
    from sklearn.metrics import mean_squared_error, r2_score 

    rmse = np.sqrt(mean_squared_error(y_test,y_pred))
    r2 = r2_score(y_test,y_pred)

    # Scikit-learn doesn't have adjusted r-square, hence custom code
    n = y_pred.shape[0]
    k = X_train.shape[1]
    adj_r_sq = 1 - (1 - r2)*(n-1)/(n-1-k)

    print(rmse, r2, adj_r_sq)
    
        
measure_performance(train_labels,y_pred, train_features)<jupyter_output>273.0375941794084 0.5227670176615185 0.5221559639069238<jupyter_text>Linear regression<jupyter_code>linear_regressor = LinearRegression()  # create object for the class
linear_regressor.fit(train_features, train_labels)  # perform linear regression
y_pred = linear_regressor.predict(train_features)  # make predictions
print(feature_list_new)
print(linear_regressor.coef_)
linear_regressor.score(train_features, train_labels, sample_weight=None)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10410681088389387, 'fraction_numerical': 0.012337333107993917, 'mean_word_length': 2.5381165919282513, 'pattern_counts': {'":': 9, '<': 13, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '533029', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3901, 'n_tokens_neox': 3671, 'n_words': 611, 'path': 'infrarisk/experiments/micropolis_ml_model/final_data_for_ml.ipynb', 'repo_name': 'srijithbalakrishnan/dreaminsg-integrated-model', 'size': '11834'}