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



# ## Dataset
#
# Sources: [UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/climate+model+simulation+crashes)
#
#
# Abstract: Given Latin hypercube samples of 18 climate model input parameter values, predict climate model simulation crashes and determine the parameter value combinations that cause the failures.
#
# This dataset contains records of simulation crashes encountered during climate model uncertainty quantification (UQ) ensembles. 
#
# Ensemble members were constructed using a Latin hypercube method in LLNL's UQ Pipeline software system to sample the uncertainties of 18 model parameters within the Parallel Ocean Program (POP2) component of the Community Climate System Model (CCSM4). 
#
# Three separate Latin hypercube ensembles were conducted, each containing 180 ensemble members. 46 out of the 540 simulations failed for numerical reasons at combinations of parameter values. 
#
# The goal is to use classification to predict simulation outcomes (fail or succeed) from input parameter values, and to use sensitivity analysis and feature selection to determine the causes of simulation crashes. 
#
# Further details about the data and methods are given in the publication 'Failure Analysis of Parameter-Induced Simulation Crashes in Climate Models,' Geoscientific Model Development ([Web Link]).
#
# **Attribute Information:**
#
# The goal is to predict climate model simulation outcomes (column 21, fail or succeed) given scaled values of climate model input parameters (columns 3-20). 
#
# Column 1: Latin hypercube study ID (study 1 to study 3) 
#
# Column 2: simulation ID (run 1 to run 180) 
#
# Columns 3-20: values of 18 climate model parameters scaled in the interval [0, 1] 
#
# Column 21: simulation outcome (0 = failure, 1 = success)
#
# **NB: For this dataset there is no necessity of scaling the data and converting to any other form**

# ## Import and overview on the dataset

# +
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt 
import seaborn as sns

df_big = pd.read_table("pop_failures.dat", sep="\s+")
#np.loadtxt('pop_failures.dat', unpack = True, skiprows=1)
df_big.head()
# -

# ## Correlation matrix
# Note: forced in the range [-0.03,0.03]

# +
data = df_big.drop(columns=['Study', 'Run', 'outcome'])
features = data.columns
target = df_big['outcome']

fig = plt.figure(figsize=(10,10))
df_big.isnull().sum() # check for missing values in all the columns
sns.set(rc={'figure.figsize':(11.7,8.27)})
correlation_matrix = data.corr().round(2)
sns.heatmap(data=correlation_matrix, annot=True, vmin = -0.03, vmax = 0.03, cmap="RdBu_r") # annot = True to print the values inside the square
plt.title("Correlation matrix", size=20)
fig.tight_layout()
plt.savefig('Figures/Dataset/correlationmatrix.png', dpi=100)
plt.show()
# -

# ## Distribution of success and failure - 1D

# +
fig, axes = plt.subplots(int(len(features)),2,figsize=(10,len(features)*3))
failures = data[target == 0]
successes = data[target == 1]
ax = axes.ravel()
colors = np.where(target==0,'tomato','seagreen')

c=0
for i,col in enumerate(features):
    ax[c].scatter(data[col], target, marker='o', c=colors, alpha =0.5)
    ax[c].set_title(col)
    ax[c].set_xlabel(col)
    ax[c].set_xlabel("Feature magnitude")
    ax[c].set_yticks([0.,1.])
    ax[c].set_yticklabels(["failure", "success"])
    c +=1
    _, bins = np.histogram(data[col], bins =50)
    ax[c].hist(successes[col], bins = bins, alpha = 0.5, color="seagreen")
    ax[c].hist(failures[col], bins = bins, alpha = 0.5, color="tomato")
    ax[c].set_xlabel("Feature magnitude")
    ax[c].set_title(col)
    ax[c].set_yticks(())
    ax[c].legend(["success", "failure"])#, loc ="best")
    c +=1

plt.suptitle("Frequency and distributions of success-failure", size=20, y=1)
fig.tight_layout()
plt.grid(False)
plt.savefig('Figures/Dataset/distributions.png', dpi=100)
plt.show()
# -

# # Selection of most sensible parameters
#
# As we can see form the previous plots the most sensible paramters are 'vconst_corr', 'vconst_2', 'convect_corr', 'bckgrnd_vdc1' (also said in the paper [Lucas D.D. at al.](https://gmd.copernicus.org/articles/6/1157/2013/)). So a selection will be made just on this four variables to create relevant plots for the report.

# +
features_sel = [ 'vconst_corr', 'vconst_2', 'convect_corr', 'bckgrnd_vdc1']

fig, axes = plt.subplots(int(len(features_sel)),2,figsize=(10,len(features_sel)*3))
ax = axes.ravel()

c=0
for i,col in enumerate(features_sel):
    ax[c].scatter(data[col], target, marker='o', c=colors, alpha =0.5)
    ax[c].set_title(col)
    ax[c].set_xlabel(col)
    ax[c].set_xlabel("Feature magnitude")
    ax[c].set_yticks([0.,1.])
    ax[c].set_yticklabels(["failure", "success"])
    c +=1
    _, bins = np.histogram(data[col], bins =50)
    ax[c].hist(successes[col], bins = bins, alpha = 0.5, color="seagreen")
    ax[c].hist(failures[col], bins = bins, alpha = 0.5, color="tomato")
    ax[c].set_xlabel("Feature magnitude")
    ax[c].set_title(col)
    ax[c].set_yticks(())
    ax[c].legend(["success", "failure"])#, loc ="best")
    c +=1

plt.suptitle("Frequency and distributions of success-failure\n(4/18 selected features)", size=20, y=1)
fig.tight_layout()
plt.grid(False)
plt.savefig('Figures/Dataset/distributions_sel.png',dpi=100)
plt.show()

# +
fig, axes = plt.subplots(2,3,figsize=(12,8))
ax = axes.ravel()

fs=features_sel
c=0
for i in range(len(fs)):

    for j in range(i+1,len(fs)):

        ax[c].scatter(successes[fs[i]], successes[fs[j]], marker='o', s=50,color="seagreen", alpha =0.5, label ="success")
        ax[c].scatter(failures[fs[i]], failures[fs[j]], marker='o', s=80,color="tomato", alpha =0.6, label = "failure")
        #ax[c].set_title(col)
        ax[c].set_xlabel(fs[i])
        ax[c].set_ylabel(fs[j])
        ax[c].legend(loc =3)
        c += 1

plt.suptitle("Joint frequency of success-failure\n(4/18 selected)", size=20, y=1)
fig.tight_layout()
plt.grid(False)
plt.savefig('Figures/Dataset/distributions_2D_sel.png',dpi =100)
plt.show()
# -

# ## Principal Component Analysis
#
# Since there is this little correlation, first step is to have an insight on how much we can condence the information, setting different values of the wanted cumulative variance condensed in the $n$ variables. The more $n$ is large the more the PCA looses its effectiveness.

# +
from sklearn.decomposition import PCA, KernelPCA
from sklearn.preprocessing import StandardScaler

x = data.values
x = StandardScaler().fit_transform(x) # normalizing the features

pca = PCA(n_components=0.95)
x_reduced = pca.fit_transform(x)

print("Variance of 95%:", np.shape(x_reduced)[1], "PC")

pca = PCA(n_components=0.70)
x_reduced = pca.fit_transform(x)

print("Variance of 70%:", np.shape(x_reduced)[1], "PC")

pca = PCA(n_components=0.55)
x_reduced = pca.fit_transform(x)

print("Variance of 55%:", np.shape(x_reduced)[1], "PC")
# -

# With a cumulative variance of 95%, all the 18 features are kept. Even lowering the variance the result doesn't get any better. There is no point to do the PCA when the fautures are so little correlated, otherwise we loose too much information.

StarcoderdataJupyterScriptsDedupFiltered