<jupyter_start><jupyter_text>Setup<jupyter_code># Standard imports
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import time
import dill
import json
import os

from scipy.stats import pearsonr

import torch
from torch.autograd import Variable
from torch.nn import functional as F
from torch.utils.data import DataLoader
from torch import nn

import ipywidgets as widgets
from ipywidgets import interact, interactive, fixed, interact_manual
from IPython.display import display

from sklearn import metrics as skmetrics
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA


from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

from captum.attr import IntegratedGradients
from captum.attr import LayerConductance
from captum.attr import NeuronConductance
# Custom utilities imports
sys.path.append("/media/krzysztof/Nowy/Doktorat - Modelling drug efficacy in cancer/Recommender System Approach/Scripts/Modules")
from modeling import Dataset, Results, Model, ModelWithAutoencoders, ModelAnalysis, EarlyStopping
import architectures<jupyter_output><empty_output><jupyter_text>Load dataframe with setups<jupyter_code># setups_df = pd.read_csv("./Results/setups.csv")
# print(setups_df.shape)
# setups_df<jupyter_output><empty_output><jupyter_text>Train and evaluate new model  Load dataset<jupyter_code>filepath = "../../Data/Preprocessed Datasets/"
with open(filepath + "GDSC-KINOMEscan_proteins_intersection_+_remaining_GDSC_target_genes_dataset.pkl", "rb") as f:
    full_dataset = dill.load(f)
print(full_dataset.name, type(full_dataset))
print()
print(full_dataset.description)
print(full_dataset.full_cell_lines_data.shape, full_dataset.drugs_data.shape)<jupyter_output>Kinases Dataset + Remaning GDSC drug's putative targets <class 'modeling.Dataset'>

Dataset containing 74 common drugs of GDSC and HMS LINCS Kinome scan dataset.
Cell lines data types: expression, coding variant and tissue type. Expressions and coding variants are 
present only for proteins present in both GDSC and KINOMEscan data, resulting in expression of 188 genes and
nmutations in 18 genes.
In addition, expressions and mutations (17 new features) of remaining target genes from GDSC are included
Tissue types are dummy encoded GDSC Tissue Descriptions 1 (18 features).
Drugs representation: inhibition scores (% control) of 294 proteins. Set of proteins is the intersection of 
proteins screened for each of 74 drugs.
Drug response data: drug reponse data contains AUC metrics across cell lines for 74 drugs considered.
(922, 242) (74, 294)<jupyter_text>Data preprocessing<jupyter_code># Establish response data for samples (drug-cell line pairs)
response_df = full_dataset.response_data.copy()

# Establish cell line features data
cell_line_data_original_df = full_dataset.full_cell_lines_data.copy()

# Search for cell lines present in response data, but missing the genomic features
missing_cell_lines = []
for cosmic_id in response_df.COSMIC_ID.unique():
    if cosmic_id not in cell_line_data_original_df.cell_line_id.unique():
        missing_cell_lines.append(cosmic_id)
# Put cell line IDs into index and drop cell line IDs columns
cell_line_data_original_df.index = cell_line_data_original_df.cell_line_id
cell_line_data_original_df = cell_line_data_original_df.drop("cell_line_id", axis=1)

# Extract response only for cell lines for which features are present
response_df = response_df[~response_df.COSMIC_ID.isin(missing_cell_lines)]

# Establish drug features data
drug_data_original_df = full_dataset.drugs_data.copy()

# Convert drug index from LINCS name to GDSC drug ID
drug_data_original_df.index = drug_data_original_df.index.map(full_dataset.kinomescan_name_to_gdsc_id_mapper)

# Establish how many cells go to train/val/test sets
num_val_cell_lines = 100
num_test_cell_lines = 100

print(drug_data_original_df.shape, cell_line_data_original_df.shape)
# response_df["AUC"].describe()
# Establish train, validation, and test sets
split_seed = 40   # ESTABLISH SPLIT SEED
# Split data into train/val/test sets
samples_train, samples_val, samples_test, cell_lines_test, cell_lines_val = Dataset.samples_train_test_split(
                                                                        response_df,
                                                                        num_val_cell_lines,
                                                                        num_test_cell_lines,
                                                                        split_seed,
                                                                        shuffle=True)
# Merge training and validation samples
samples_train_and_val = pd.concat([samples_train, samples_val], axis=0)
                                                                       
# Normalize the data
# Cell line data
cols_subset = [col for col in list(cell_line_data_original_df) if col.endswith("_exp")]
rows_subset = [x for x in cell_line_data_original_df.index if x not in cell_lines_test + cell_lines_val]

cell_line_data_df = Dataset.standardize_data(cell_line_data_original_df, cols_subset=cols_subset,
                                            rows_subset=rows_subset)
# Drug data
drug_data_df = Dataset.standardize_data(drug_data_original_df)<jupyter_output><empty_output><jupyter_text>Exclude drugs if needed<jupyter_code># drugs_to_exclude = [110, 1038, 1158, 1614, 1632, 54, 59, 1052, 253]
# print(len(drugs_to_exclude))
# print(samples_train["DRUG_ID"].nunique(), samples_val["DRUG_ID"].nunique(),
#      samples_train_and_val["DRUG_ID"].nunique(), samples_test["DRUG_ID"].nunique())
# samples_train = Dataset.exclude_drugs(samples_train, drugs_to_exclude)
# samples_val = Dataset.exclude_drugs(samples_val, drugs_to_exclude)
# samples_train_and_val = Dataset.exclude_drugs(samples_train_and_val, drugs_to_exclude)
# samples_test = Dataset.exclude_drugs(samples_test, drugs_to_exclude)

# print(samples_train["DRUG_ID"].nunique(), samples_val["DRUG_ID"].nunique(),
#      samples_train_and_val["DRUG_ID"].nunique(), samples_test["DRUG_ID"].nunique())<jupyter_output><empty_output><jupyter_text>Definition of new PyTorch network if needed  Instantiate the new system<jupyter_code>state_dict = None
torch.manual_seed(11)
def establish_new_system(drug_dim, cell_line_dim, code_dim, architecture_type,
                    num_layers, drug_hidden_dim1, drug_hidden_dim2, drug_hidden_dim3,
                    cell_line_hidden_dim1, cell_line_hidden_dim2, cell_line_hidden_dim3,
                     dropout, dropout_rate, code_activation, activation_func, out_activation,
                        batch_norm):
    # Establish specs
    specs = {"drug_dim": drug_dim,
            "cell_line_dim": cell_line_dim,
            "code_dim": code_dim,
            "architecture_type": architecture_type,
            "num_layers": num_layers,
            "drug_hidden_dim1": drug_hidden_dim1,
            "drug_hidden_dim2": drug_hidden_dim2,
            "drug_hidden_dim3": drug_hidden_dim3,
            "cell_line_hidden_dim1": cell_line_hidden_dim1,
            "cell_line_hidden_dim2": cell_line_hidden_dim2,
            "cell_line_hidden_dim3": cell_line_hidden_dim3,
            "dropout": dropout,
            "dropout_rate": dropout_rate,
            "code_activation": code_activation,
            "activation_func": activation_func,
            "out_activation": out_activation,
            "batch_norm": batch_norm}
    
    
    network = ModelAnalysis.instantiate_system(specs, state_dict)
    return network

establish_new_system_widget = interactive(establish_new_system,
                                     drug_dim=[294],
                                     cell_line_dim=[241, 1144],
                                     code_dim=widgets.IntSlider(value=10, min=2, max=30),
                                     architecture_type=["autoencoders", "linear"],
                                     num_layers=[1, 2, 3, None],
                                     drug_hidden_dim1=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     drug_hidden_dim2=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     drug_hidden_dim3=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     cell_line_hidden_dim1=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     cell_line_hidden_dim2=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     cell_line_hidden_dim3=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     dropout=[False, True],
                                     dropout_rate=widgets.FloatSlider(value=0.5, min=0, max=1, step=0.05),
                                     code_activation=[False, True],
                                     activation_func=[nn.ReLU, nn.Tanh, nn.SELU],
                                     out_activation=[torch.sigmoid, None],
                                     batch_norm=[False, True])
display(establish_new_system_widget)
# MUST EXECUTE ALSO
network = establish_new_system_widget.result
print(network.out_activation)
print(next(network.parameters()).sum())
network<jupyter_output><built-in method sigmoid of type object at 0x7fd2a6488640>
tensor(4.2539, grad_fn=<SumBackward0>)<jupyter_text>Add code forwarding if needed<jupyter_code># Modify network to add forward propagation of codes if needed
def forward_network_specs(model_type, num_layers, code_interactions, code_dim,
                         out_activation, forward_net_hidden_dim1,
                         forward_net_hidden_dim2, drop_rate):
    if model_type == "linear":
        if code_interactions:
            m = architectures.ForwardLinearRegression(code_dim ** 2 + 2 * code_dim, 
                                            out_activation=out_activation)
        else:
            m = architectures.ForwardLinearRegression(2 * code_dim, 
                                                      out_activation=out_activation)
    else:
        if code_interactions:
            if num_layers == 1:
                m = architectures.ForwardNetworkOneHiddenLayer(code_dim ** 2 + 2 * code_dim,
                                                          forward_net_hidden_dim1,
                                                          dropout_rate=drop_rate,
                                                          out_activation=out_activation)
            elif num_layers == 2:
                m = architectures.ForwardNetworkTwoHiddenLayers(code_dim ** 2 + 2 * code_dim,
                                                          forward_net_hidden_dim1,
                                                          forward_net_hidden_dim2,
                                                          dropout_rate=drop_rate,
                                                          out_activation=out_activation)
        else:
            if num_layers == 1:
                m = architectures.ForwardNetworkOneHiddenLayer(2 * code_dim,
                                                          forward_net_hidden_dim1,
                                                          dropout_rate=drop_rate,
                                                          out_activation=out_activation)
            elif num_layers == 2:
                m = architectures.ForwardNetworkTwoHiddenLayers(2 * code_dim,
                                                          forward_net_hidden_dim1,
                                                          forward_net_hidden_dim2,
                                                          dropout_rate=drop_rate,
                                                          out_activation=out_activation)
    return m, code_interactions

forward_net_specs_widget = interactive(forward_network_specs,
                                       model_type=["linear", "nn"],
                                       num_layers=[1, 2, 3],
                                       code_interactions=[False, True],
                                       code_dim=widgets.IntSlider(value=10, min=2, max=30),
                                       out_activation=[None, torch.sigmoid],
                                       forward_net_hidden_dim1=widgets.BoundedIntText(value=128, min=16, max=2048),
                                       forward_net_hidden_dim2=widgets.BoundedIntText(value=128, min=16, max=2048),
                                       drop_rate=widgets.FloatSlider(value=0.5, min=0, max=1, step=0.05))
                                       
display(forward_net_specs_widget)<jupyter_output><empty_output><jupyter_text>Pre-train autoencoders<jupyter_code>print(next(network.drug_autoencoder.parameters()).sum())
print(next(network.cell_line_autoencoder.parameters()).sum())
drug_autoencoder = network.drug_autoencoder
cell_line_autoencoder = network.cell_line_autoencoder
print(next(drug_autoencoder.parameters()).sum())
print(next(cell_line_autoencoder.parameters()).sum())
# Batch training for drug autencoder
learning_rate = 0.001
num_epochs = 100
optimizer = torch.optim.Adam(drug_autoencoder.parameters(), lr=learning_rate)
reconstruction_criterion = nn.MSELoss()

for epoch in range(1, num_epochs + 1):
    drug_autoencoder.train()
    drug_input = torch.from_numpy(drug_data_df.values)
    drug_codes, drug_reconstruction = drug_autoencoder(drug_input.float())
    loss = reconstruction_criterion(drug_reconstruction, drug_input.float())
    optimizer.zero_grad()
    # Get the gradients w.r.t. the parameters
    loss.backward()
    # Update the parameters
    optimizer.step()
    
    print("Epoch: {}, loss: {}".format(epoch, loss.item()))

print(drug_input.shape, drug_reconstruction.shape)
print("Whole correlation:", pearsonr(drug_input.numpy().flatten(),
                  drug_reconstruction.detach().numpy().flatten()))
# Per drug corr
corrs = []
for i in range(drug_input.shape[0]):
    drug_original = drug_input[i]
    drug_rec = drug_reconstruction[i]
    corrs.append(pearsonr(drug_original.numpy(), drug_rec.detach().numpy())[0])
print("Median corr per drug:", np.median(corrs))
# Batch training for cell line autencoder
learning_rate = 0.001
num_epochs = 100
optimizer = torch.optim.Adam(cell_line_autoencoder.parameters(), lr=learning_rate)

cell_lines_train = list(set(cell_line_data_original_df.index) - set(cell_lines_test) - set(cell_lines_val))
cl_input_df = cell_line_data_df.loc[cell_lines_train]

for epoch in range(1, num_epochs + 1):
    cell_line_autoencoder.train()
    cl_input = torch.from_numpy(cl_input_df.values)
    cl_codes, cl_reconstruction = cell_line_autoencoder(cl_input.float())
    loss = reconstruction_criterion(cl_reconstruction, cl_input.float())
    optimizer.zero_grad()
    # Get the gradients w.r.t. the parameters
    loss.backward()
    # Update the parameters
    optimizer.step()
    
    print("Epoch: {}, loss: {}".format(epoch, loss.item()))
    
print(cl_input.shape, cl_reconstruction.shape)
print("Whole correlation:", pearsonr(cl_input.numpy().flatten(),
                  cl_reconstruction.detach().numpy().flatten()))
# Per drug corr
corrs = []
for i in range(cl_input.shape[0]):
    cl_original = cl_input[i]
    cl_rec = cl_reconstruction[i]
    corrs.append(pearsonr(cl_original.numpy(), cl_rec.detach().numpy())[0])
print("Median corr per cell line:", np.median(corrs))
print(next(network.drug_autoencoder.parameters()).sum())
print(next(network.cell_line_autoencoder.parameters()).sum())

print(next(drug_autoencoder.parameters()).sum())
print(next(cell_line_autoencoder.parameters()).sum())<jupyter_output>tensor(-35.8564, grad_fn=<SumBackward0>)
tensor(81.2122, grad_fn=<SumBackward0>)
tensor(-35.8564, grad_fn=<SumBackward0>)
tensor(81.2122, grad_fn=<SumBackward0>)<jupyter_text>Train<jupyter_code># Training params
num_epochs = 100
batch_size = 256
no_batches = samples_train.shape[0] // batch_size + 1
# Establish the device
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
learning_rate = 0.00005
reconstruction_term_drug = 0.0
reconstruction_term_cl = 0.
independence_term_drug = 0.
independence_term_cl = 0.
reg_lambda = 0.0

augment_data = False

if augment_data:
    balanced_augmentation = True
    if balanced_augmentation:
        auc_threshold = 0.85
    augmented_data_factor = 2
    cell_line_exp_idx = 202
    cell_line_noise_mean = 0.0
    cell_line_noise_std = 0.6
    target_variable_noise_mean = 0.
    target_variable_noise_std = 0.15


optimizer = torch.optim.Adam(network.parameters(), lr=learning_rate)

weigh_by_true_auc = False
if weigh_by_true_auc:
    output_criterion = nn.MSELoss(reduction="none")
else:
    output_criterion = nn.MSELoss()
    
reconstruction_criterion = nn.MSELoss()

print(no_batches, device)
early_stopping = EarlyStopping(patience=num_epochs, delta=0.000)
# Establish train samples
train_samples = samples_train.copy()
print(train_samples.shape)

np.random.seed(11)

# model = ModelWithAutoencoders("Test", network)
train_losses = []
train_y_losses = []
drug_independence_train_losses = []
cl_independence_train_losses = []
train_rmses = []
train_corrs = []
train_median_corrs = []
train_drug_rec_corrs = []
train_cl_rec_corrs = []
val_rmses = []
val_corrs = []
val_median_corrs = []
val_drug_rec_corrs = []
val_cl_rec_corrs = []

# Establish the device
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
print(device)
# Move the network into device
network.to(device)
# Training the model
network.train()

for epoch in range(1, num_epochs + 1):
    # train_samples = train_samples.sample(frac=1., random_state = 11)
    # Initialize lists for storing training losses during batches
    epoch_main_losses = 0
    epoch_y_losses = 0
    epoch_drug_indep_losses = 0
    epoch_cl_indep_losses = 0
    for batch in range(no_batches):
        # Separate response variable batch
        if batch != no_batches:
            samples_batch = train_samples.iloc[batch * batch_size:(batch + 1) * batch_size]
        else:
            samples_batch = train_samples.iloc[batch * batch_size:]

        # Extract output variable batch
        y_batch = samples_batch["AUC"].values.reshape(-1, 1)
        
        # Extract cell lines IDs for which data shall be extracted
        cl_ids = samples_batch["COSMIC_ID"].values
        # Extract corresponding cell line data
        cell_line_input_batch = cell_line_data_df.loc[cl_ids].values

        # Extract drug IDs for which data shall be extracted
        drug_ids = samples_batch["DRUG_ID"].values
        # Extract corresponding drug data
        drug_input_batch = drug_data_df.loc[drug_ids].values
        
        # Augment data if specified
        if augment_data:
            original_drug_batch = drug_input_batch.copy()
            original_cell_line_batch = cell_line_input_batch.copy()
            original_y_batch = y_batch.copy()
            
            if balanced_augmentation:
                indexes = (original_y_batch < auc_threshold).reshape(-1)
                
                original_drug_batch = original_drug_batch[indexes, :]
                original_cell_line_batch = original_cell_line_batch[indexes, :]
                original_y_batch = original_y_batch[indexes, :]

            for i in range(augmented_data_factor):
                augmented_cell_line_data = original_cell_line_batch[:, :cell_line_exp_idx]
                cell_line_noise = np.random.normal(loc=cell_line_noise_mean,
                                                  scale=cell_line_noise_std,
                                                  size=augmented_cell_line_data.shape)
                augmented_cell_line_data = augmented_cell_line_data + cell_line_noise
                augmented_cell_line_data = np.concatenate((augmented_cell_line_data, 
                                                           original_cell_line_batch[:, cell_line_exp_idx:]),
                                                          axis=1)
                augmented_response_data = original_y_batch + np.random.normal(loc=target_variable_noise_mean,
                                                            scale=target_variable_noise_std,
                                                            size=original_y_batch.shape)
                cell_line_input_batch = np.concatenate((cell_line_input_batch, augmented_cell_line_data),
                                                      axis=0)
                y_batch = np.concatenate((y_batch, augmented_response_data), axis=0)
                drug_input_batch = np.concatenate((drug_input_batch, original_drug_batch), axis=0)
                

        # Clear gradient buffers because we don't want to accummulate gradients 
        optimizer.zero_grad()
        
        # Put data batches into torch and device
        cell_line_input_batch = torch.from_numpy(cell_line_input_batch).to(device)
        drug_input_batch = torch.from_numpy(drug_input_batch).to(device)
        y_batch = torch.from_numpy(y_batch).to(device)

        # Perform forward pass
        batch_output, batch_drug_reconstruction, batch_cl_reconstruction = network(
            drug_input_batch.float(), cell_line_input_batch.float())

        # L2 regularization
        reg_sum = 0
        for param in network.parameters():
            reg_sum += 0.5 * (param ** 2).sum()  # L2 norm

        # Compute the loss for this batch, including the drug and cell line reconstruction losses
        # Variant wit weightning the samples
        if weigh_by_true_auc:
            output_loss = output_criterion(batch_output, y_batch.float())
            output_loss = torch.mul(output_loss, y_batch.float() ** 0.5)
            output_loss = output_loss.mean() + reg_lambda * reg_sum
        # Variant with normal loss
        else:
            output_loss = output_criterion(batch_output, y_batch.float()) + reg_lambda * reg_sum
        # Reconstruction losses
        drug_recounstruction_loss = reconstruction_criterion(batch_drug_reconstruction, drug_input_batch.float())
        cl_reconstruction_loss = reconstruction_criterion(batch_cl_reconstruction, cell_line_input_batch.float())

        # Compute independence loss
        # Covariance matrices
        drug_codes_batch = network.drug_autoencoder.encoder(drug_input_batch.float())
        cl_codes_batch = network.cell_line_autoencoder.encoder(cell_line_input_batch.float())
        drug_cov = ModelWithAutoencoders.covariance_matrix_torch(drug_codes_batch)
        cl_cov = ModelWithAutoencoders.covariance_matrix_torch(cl_codes_batch)
        
        # Actual dependence losses
        drug_independence_loss = (drug_cov * drug_cov).sum() - torch.trace(drug_cov * drug_cov)
        cl_independence_loss = (cl_cov * cl_cov).sum() - torch.trace(cl_cov * cl_cov)

        # Combine the losses in the final cost function
        loss = output_loss + reconstruction_term_drug * drug_recounstruction_loss + \
                reconstruction_term_cl * cl_reconstruction_loss + \
                independence_term_drug * drug_independence_loss + independence_term_cl * cl_independence_loss
        
        # Get the gradients w.r.t. the parameters
        loss.backward()
        # Update the parameters
        optimizer.step()
        
        epoch_main_losses += loss.item()
        epoch_y_losses += output_loss.item()
        epoch_drug_indep_losses += drug_independence_loss.item()
        epoch_cl_indep_losses += cl_independence_loss.item()
    
    f = lambda x: np.round(x, 3)
    print("Epoch: {}, main loss: {}, output loss: {} drug independence loss: {}, cell line independence loss: {}".format(
        epoch, f(epoch_main_losses / no_batches), f(epoch_y_losses / no_batches),
        f(epoch_drug_indep_losses / no_batches), f(epoch_cl_indep_losses / no_batches)))
    
    # Update metrics for learning curve
    train_losses.append(epoch_main_losses / no_batches)
    train_y_losses.append(epoch_y_losses / no_batches)
    drug_independence_train_losses.append(epoch_drug_indep_losses / no_batches)
    cl_independence_train_losses.append(epoch_cl_indep_losses / no_batches)
    
    # Evaluate on training data
    model = ModelWithAutoencoders("Example", network)
    predicted, y_true, drug_input, cl_input = model.predict(train_samples, cell_line_data_df, drug_data_df)
    preds, drug_reconstruction, cl_reconstruction = predicted
    try:
        train_rmse, train_corr = Model.evaluate_predictions(y_true, preds.cpu().numpy().reshape(-1))
    except ValueError:
        train_rmse, train_corr = np.nan, (np.nan, np.nan)
    # Drug reconstruction training error
    try:
        train_drug_rec_rmse, train_drug_rec_corr = Model.evaluate_predictions(drug_input.flatten(),
                                        drug_reconstruction.cpu().numpy().flatten())
    except ValueError:
        train_drug_rec_rmse, train_drug_rec_corr = np.nan, (np.nan, np.nan)

    # Cell line reconstruction training error
    try:
        train_cl_rec_rmse, train_cl_rec_corr = Model.evaluate_predictions(cl_input.flatten(),
                                        cl_reconstruction.cpu().numpy().flatten())
    except ValueError:
        train_cl_rec_rmse, train_cl_rec_corr = np.nan, (np.nan, np.nan)

    performance_df = Model.per_drug_performance_df(train_samples, preds.cpu())
    train_median_rmse = performance_df["Model RMSE"].median()
    train_median_corr = performance_df["Model correlation"].median()
    
    train_rmses.append(train_rmse)
    train_corrs.append(train_corr[0])
    train_median_corrs.append(train_median_corr)
    train_drug_rec_corrs.append(train_drug_rec_corr[0])
    train_cl_rec_corrs.append(train_cl_rec_corr[0])
    
    # Evaluate on validation data
    predicted, y_true, drug_input, cl_input = model.predict(samples_val, cell_line_data_df, drug_data_df)
    preds, drug_reconstruction, cl_reconstruction = predicted
    
    # Main validation error
    try:
        val_rmse, val_corr = Model.evaluate_predictions(y_true, preds.cpu().numpy().reshape(-1))
    except ValueError:
        val_rmse, val_corr = np.nan, (np.nan, np.nan)

    # Drug reconstruction validation error
    try:
        val_drug_rec_rmse, val_drug_rec_corr = Model.evaluate_predictions(drug_input.flatten(),
                                        drug_reconstruction.cpu().numpy().flatten())
    except ValueError:
        val_drug_rec_rmse, val_drug_rec_corr = np.nan, (np.nan, np.nan)

    # Cell line reconstruction validation error
    try:
        val_cl_rec_rmse, val_cl_rec_corr = Model.evaluate_predictions(cl_input.flatten(),
                                        cl_reconstruction.cpu().numpy().flatten())
    except ValueError:
        val_cl_rec_rmse, val_cl_rec_corr = np.nan, (np.nan, np.nan)
    # Per-drug validation evaluation
    performance_df = Model.per_drug_performance_df(samples_val, preds.cpu())
    val_median_corr = performance_df["Model correlation"].median()
        
    val_rmses.append(val_rmse)
    val_corrs.append(val_corr[0])
    val_drug_rec_corrs.append(val_drug_rec_corr[0])
    val_cl_rec_corrs.append(val_cl_rec_corr[0])
    val_median_corrs.append(val_median_corr)
    
    early_stopping(val_rmse, network)
        
    if early_stopping.early_stop:
        print("Early stopping")
        break
cell_line_input_batch.shape
def learning_curve(metric1, metric2):
    """Plot learning curve related to a given metric."""
    dic = {"Whole training loss": train_losses,
           "Y loss": train_y_losses,
          "Drug independence train loss": drug_independence_train_losses,
          "Cl independence train loss": cl_independence_train_losses,
           "Train RMSE": train_rmses,
           "Train corr.": train_corrs,
           "Train median corr.": train_median_corrs,
           "Train drug rec corr.": train_drug_rec_corrs,
           "Train cl rec corr.": train_cl_rec_corrs,
          "Val RMSE": val_rmses,
          "Val corr.": val_corrs,
          "Val median corr.": val_median_corrs,
          "Val drug rec corr.": val_drug_rec_corrs,
          "Val cl rec corr.": val_cl_rec_corrs}
    
    x = range(1, len(train_losses) + 1)
    print("Max {}: {}".format(metric1, max(dic[metric1])))
    sns.lineplot(x, dic[metric1], label=metric1)
    if metric2:
        print("Max {}: {}".format(metric2, max(dic[metric2])))
        sns.lineplot(x, dic[metric2], label=metric2)
    
    plt.grid()
    plt.show()
    
metrics = ["Whole training loss", "Y loss", "Drug independence train loss", "Cl independence train loss",
           "Train RMSE", "Train corr.", "Train median corr.", "Train drug rec corr.", "Train cl rec corr.",
          "Val RMSE", "Val corr.", "Val median corr.", "Val drug rec corr.", "Val cl rec corr."]

interact(learning_curve,
        metric1=metrics,
        metric2=[None] + metrics)<jupyter_output><empty_output><jupyter_text>Evaluate response predictions<jupyter_code># Evaluation on training set
samples = samples_train.copy()
print("Training data")
print()
# Make predictions
# Linear model case
if model.__class__ == Model:
    response_preds, y_true = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    print(response_preds.shape, y_true.shape)
# Model with autoencoders case
elif model.__class__ == ModelWithAutoencoders:
    predicted, y_true, drug_input, cl_input = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    response_preds, drug_rec_preds, cl_rec_preds = predicted
    print(response_preds.shape, y_true.shape, drug_input.shape, cl_input.shape)

# Create a dataframe with samples and corresponding predictions
samples_with_predictions_train_df = samples.copy()
samples_with_predictions_train_df.insert(3, "Predicted AUC", response_preds.reshape(-1))
print(samples_with_predictions_train_df.shape)

# Evaluate predictions numerically
rmse, corr = Model.evaluate_predictions(y_true, response_preds.numpy().reshape(-1))
print("Response prediction performance:")
print("RMSE: {}, Corr.: {}".format(rmse, corr))

per_drug_performance_train_df = Model.per_drug_performance_df(samples, response_preds)
print(per_drug_performance_train_df.shape)
per_drug_performance_train_df.insert(1, "Drug name", per_drug_performance_train_df["Drug ID"].map(
    full_dataset.gdsc_drug_id_to_gdsc_name_mapper))
print(per_drug_performance_train_df.shape)
per_drug_performance_train_df = per_drug_performance_train_df.sort_values("Model correlation", ascending=False)
print("Per drug median RMSE and correlation:", per_drug_performance_train_df["Model RMSE"].median(),
     per_drug_performance_train_df["Model correlation"].median())
print()


# Evaluate reconstruction predictions
# Metrics calculated on whole matrix
rmse, corr = Model.evaluate_predictions(drug_input.ravel(), drug_rec_preds.numpy().ravel())
print("Reconstruction performance:")
print("Whole drug data RMSE and correlation:", np.round(rmse, 3), corr)
rmse, corr = Model.evaluate_predictions(cl_input.ravel(), cl_rec_preds.numpy().ravel())
print("Whole cell line data RMSE and correlation:", np.round(rmse, 3), corr)

per_drug_performance_train_df.head(2)
samples_with_predictions_train_df["Predicted AUC"].describe()
# Evaluation on validation set
samples = samples_val.copy()
print("Validation data")
print()
# Make predictions
# Linear model case
if model.__class__ == Model:
    response_preds, y_true = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    print(response_preds.shape, y_true.shape)
# Model with autoencoders case
elif model.__class__ == ModelWithAutoencoders:
    predicted, y_true, drug_input, cl_input = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    response_preds, drug_rec_preds, cl_rec_preds = predicted
    print(response_preds.shape, y_true.shape, drug_input.shape, cl_input.shape)

# Create a dataframe with samples and corresponding predictions
samples_with_predictions_val_df = samples.copy()
samples_with_predictions_val_df.insert(3, "Predicted AUC", response_preds.reshape(-1))
print(samples_with_predictions_val_df.shape)

# Evaluate predictions numerically
rmse, corr = Model.evaluate_predictions(y_true, response_preds.numpy().reshape(-1))
print("Response prediction performance:")
print("RMSE: {}, Corr.: {}".format(rmse, corr))

per_drug_performance_val_df = Model.per_drug_performance_df(samples, response_preds)
print(per_drug_performance_val_df.shape)
per_drug_performance_val_df.insert(1, "Drug name", per_drug_performance_val_df["Drug ID"].map(
    full_dataset.gdsc_drug_id_to_gdsc_name_mapper))
print(per_drug_performance_val_df.shape)
per_drug_performance_val_df = per_drug_performance_val_df.sort_values("Model correlation", ascending=False)
print("Per drug median RMSE and correlation:", per_drug_performance_val_df["Model RMSE"].median(),
     per_drug_performance_val_df["Model correlation"].median())
print()


# Evaluate reconstruction predictions
# Metrics calculated on whole matrix
rmse, corr = Model.evaluate_predictions(drug_input.ravel(), drug_rec_preds.numpy().ravel())
print("Reconstruction performance:")
print("Whole drug data RMSE and correlation:", np.round(rmse, 3), corr)
rmse, corr = Model.evaluate_predictions(cl_input.ravel(), cl_rec_preds.numpy().ravel())
print("Whole cell line data RMSE and correlation:", np.round(rmse, 3), corr)

per_drug_performance_val_df.head(2)
samples_with_predictions_val_df["Predicted AUC"].describe()
np.percentile(samples_with_predictions_val_df["Predicted AUC"], 85)<jupyter_output><empty_output><jupyter_text>Investigate weights and predictions  Compute reconstructions and code embeddings for drugs and cell lines<jupyter_code># Compute reconstructions for drug
drug_codes, drug_recs = network.drug_autoencoder(torch.from_numpy(drug_data_df.values).float())
print(drug_codes.shape, drug_recs.shape)
# Create a DataFrame with drugs and their code embeddings
drug_codes_df = pd.DataFrame(drug_codes.detach().numpy(), index=drug_data_df.index)
print(drug_codes_df.shape)
drug_codes_df.head(2)
# Per-drug reconstruction performance
drug_rec_performance_df = ModelAnalysis.per_object_reconstruction_evaluations(drug_data_df,
                                                                             drug_recs)
print(drug_rec_performance_df.shape, drug_rec_performance_df["Correlation"].median())
drug_rec_performance_df.head(3)
# Compute reconstructions for cell lines
cl_codes, cl_recs = network.cell_line_autoencoder(torch.from_numpy(
    cell_line_data_df.values).float())
print(cl_codes.shape, cl_recs.shape)
cl_codes_df = pd.DataFrame(cl_codes.detach().numpy(), index=cell_line_data_df.index)
print(cl_codes_df.shape)
cl_codes_df.head(2)
# Per-cell line recosntruction performance
cl_rec_performance_df = ModelAnalysis.per_object_reconstruction_evaluations(cell_line_data_df,
                                                                             cl_recs)
print(cl_rec_performance_df.shape, cl_rec_performance_df["Correlation"].median())
cl_rec_performance_df.head(3)
response_df.shape<jupyter_output><empty_output><jupyter_text>Covariance matrices of drug and cell lines code embeddings<jupyter_code># Compute covariance matrices for drug and cell lines code representations
print(drug_codes.shape, cl_codes.shape)
drug_codes_cov = np.cov(drug_codes.detach().numpy(), rowvar=False)
cl_codes_cov = np.cov(cl_codes.detach().numpy(), rowvar=False)
print(drug_codes_cov.shape, cl_codes_cov.shape)

# Visualize drug representations covariance matrix
sns.set(font_scale=1.1)
plt.figure(figsize=(5, 4))
ax = sns.heatmap(np.abs(drug_codes_cov), annot=False)
ax.tick_params(labelsize=12)

plt.title("Dependence penalty: 0", fontsize=14)

plt.savefig("cov_matrix_drugs_exp_1_no_penalty.pdf")
plt.show()
# Visualize cell line representations covariance matrix
sns.set(font_scale=1.1)
plt.figure(figsize=(5, 4))
ax = sns.heatmap(np.abs(cl_codes_cov), annot=False)
ax.tick_params(labelsize=12)

plt.title("Dependence penalty: 0", fontsize=14)

plt.savefig("cov_matrix_cl_exp_1_no_penalty.pdf")
plt.show()
plt.show()<jupyter_output><empty_output><jupyter_text>Visualize code embeddings in 2D<jupyter_code># Visualize drug code representations with PCA
pca = PCA(n_components=2)
codes_embedded = pca.fit_transform(drug_codes.detach().numpy())
# Create a DataFrame with drugs and their 2D embeddings
drug_embeddings_df = pd.DataFrame(codes_embedded, index=drug_codes_df.index)
drug_embeddings_df["TP"] = drug_embeddings_df.index.map(full_dataset.drug_gdsc_id_to_target_pathway_mapper)
print(drug_embeddings_df.shape)
# Visualize in 2D
sns.scatterplot(x=0, y=1, data=drug_embeddings_df, hue="TP")

plt.show()
sns.scatterplot(drug_codes.detach().numpy()[:, 6], drug_codes.detach().numpy()[:, 9])

plt.show()
# Visualize cell lines code representations with PCA
pca = PCA(n_components=2)
codes_embedded = pca.fit_transform(cl_codes.detach().numpy())
# Create a DataFrame with drugs and their 2D embeddings
cl_embeddings_df = pd.DataFrame(codes_embedded, index=cl_codes_df.index)
print(cl_embeddings_df.shape)
# Visualize in 2D
sns.scatterplot(x=0, y=1, data=cl_embeddings_df)

plt.show()<jupyter_output>(922, 2)<jupyter_text>Relationship between code dimensions and input features (Neuron Attributions)  For drugs<jupyter_code>def net(drug_input):
    return model.network.drug_autoencoder(drug_input)[0]

k = 10  # How many samples use to compute importances
drug_input = torch.from_numpy(drug_data_df.values).float()[:k]
drug_input.requires_grad_()

# Integrated Gradients
ig = IntegratedGradients(net)

drug_code_attributions, delta = ig.attribute(drug_input, target=9,  return_convergence_delta=True)
drug_code_attributions = drug_code_attributions.detach().numpy()
drug_code_attributions = np.mean(drug_code_attributions, axis=0)

# Compute drug code attributions for every dimension in code space
drug_code_attributions_df = pd.DataFrame()
for dim in range(model.network.drug_autoencoder(drug_input[0])[0].shape[0]):
    code_attributions, delta = ig.attribute(drug_input, target=dim,  return_convergence_delta=True)
    code_attributions = code_attributions.detach().numpy()
    code_attributions = np.mean(code_attributions, axis=0)
    
    drug_code_attributions_df[dim] = code_attributions
drug_code_attributions_df.index = drug_data_df.columns
print(drug_code_attributions_df.shape)
drug_code_attributions_df.head()
drug_code_attributions_df.head(6)
# Visualize the codes attributions on heatmap
plt.figure(figsize=(20, 18))
sns.heatmap(drug_code_attributions_df, yticklabels=5)

plt.savefig("heatmap.pdf")
plt.show()
drug_data_df.tail()<jupyter_output><empty_output><jupyter_text>Model sensitivity analysis<jupyter_code>from captum.attr import IntegratedGradients
from captum.attr import LayerConductance
from captum.attr import NeuronConductance
# Analyze whole rec system
def net(drug_input, cl_input):
    return model.network(drug_input, cl_input)[0]
model.network
# Integrated Gradients
ig = IntegratedGradients(net)
# Establish samples on which to compute importances
k = 2   # How many samples use to compute importances
samples  = samples_train.copy()[:k]
cl_input = torch.from_numpy(cell_line_data_df.loc[samples["COSMIC_ID"].values].values).float()
drug_input = torch.from_numpy(drug_data_df.loc[samples["DRUG_ID"].values].values).float()

cl_input.requires_grad_()
drug_input.requires_grad_()

attributions, delta = ig.attribute((drug_input, cl_input), return_convergence_delta=True)
drug_attributions, cl_attributions = attributions
drug_attributions = drug_attributions.detach().numpy()
cl_attributions = cl_attributions.detach().numpy()

print(drug_attributions.shape, drug_attributions.sum())
print(cl_attributions.shape, cl_attributions.sum())

drug_attributions = np.mean(drug_attributions, axis=0)
cl_attributions = np.mean(cl_attributions, axis=0)

print(drug_attributions.shape, drug_attributions.sum())
print(cl_attributions.shape, cl_attributions.sum())

drug_attributions = [x for x in zip(list(drug_data_df.columns), drug_attributions)]
cl_attributions = [x for x in zip(list(cell_line_data_df.columns), cl_attributions)]
# Sort by importance
drug_attributions = sorted(drug_attributions, key=lambda x: abs(x[1]), reverse=True)
cl_attributions = sorted(cl_attributions, key=lambda x: abs(x[1]), reverse=True)
# Visualize k most important drug feats
k = 20
sns.barplot(x=[x[0] for x in drug_attributions[:k]], y=[x[1] for x in drug_attributions[:k]])
plt.xticks(rotation=90)
plt.show()
# Visualize k most important cl feats
k = 20
sns.barplot(x=[x[0] for x in cl_attributions[:k]], y=[x[1] for x in cl_attributions[:k]])
plt.xticks(rotation=90)
plt.show()
samples<jupyter_output><empty_output><jupyter_text>Train and evaluate new model with response variable transformation  Load dataset<jupyter_code>filepath = "../../Data/Preprocessed Datasets/"
with open(filepath + "GDSC-KINOMEscan_proteins_intersection_+_remaining_GDSC_target_genes_dataset.pkl", "rb") as f:
    full_dataset = dill.load(f)
print(full_dataset.name, type(full_dataset))
print()
print(full_dataset.description)
print(full_dataset.full_cell_lines_data.shape, full_dataset.drugs_data.shape)<jupyter_output>Kinases Dataset + Remaning GDSC drug's putative targets <class 'modeling.Dataset'>

Dataset containing 74 common drugs of GDSC and HMS LINCS Kinome scan dataset.
Cell lines data types: expression, coding variant and tissue type. Expressions and coding variants are 
present only for proteins present in both GDSC and KINOMEscan data, resulting in expression of 188 genes and
nmutations in 18 genes.
In addition, expressions and mutations (17 new features) of remaining target genes from GDSC are included
Tissue types are dummy encoded GDSC Tissue Descriptions 1 (18 features).
Drugs representation: inhibition scores (% control) of 294 proteins. Set of proteins is the intersection of 
proteins screened for each of 74 drugs.
Drug response data: drug reponse data contains AUC metrics across cell lines for 74 drugs considered.
(922, 242) (74, 294)<jupyter_text>Response variable transformation functions<jupyter_code>def one_minus_log_transform(x):
    return - np.log(1 - x)

def one_minus_log_transform_inverse(x):
    return 1 - np.exp( - x)

def log_transform(x, base="exp"):
    if base == "exp":
        return  - np.log(x)
    elif base == 2:
        return np.log2(x)
def log_transform_inverse(x):
    return  np.exp( - x)<jupyter_output><empty_output><jupyter_text>Data preprocessing<jupyter_code># Establish response data for samples (drug-cell line pairs)
response_df = full_dataset.response_data.copy()

# Transform the response variable
response_df["AUC"] = response_df["AUC"].map(log_transform)

# Establish cell line features data
cell_line_data_original_df = full_dataset.full_cell_lines_data.copy()

# Search for cell lines present in response data, but missing the genomic features
missing_cell_lines = []
for cosmic_id in response_df.COSMIC_ID.unique():
    if cosmic_id not in cell_line_data_original_df.cell_line_id.unique():
        missing_cell_lines.append(cosmic_id)
# Put cell line IDs into index and drop cell line IDs columns
cell_line_data_original_df.index = cell_line_data_original_df.cell_line_id
cell_line_data_original_df = cell_line_data_original_df.drop("cell_line_id", axis=1)

# Extract response only for cell lines for which features are present
response_df = response_df[~response_df.COSMIC_ID.isin(missing_cell_lines)]

# Establish drug features data
drug_data_original_df = full_dataset.drugs_data.copy()

# Convert drug index from LINCS name to GDSC drug ID
drug_data_original_df.index = drug_data_original_df.index.map(full_dataset.kinomescan_name_to_gdsc_id_mapper)

# Establish how many cells go to train/val/test sets
num_val_cell_lines = 100
num_test_cell_lines = 100

print(drug_data_original_df.shape, cell_line_data_original_df.shape)
response_df["AUC"].describe()
a = 0.01
print(a)
a_trans = log_transform(a)
print(a_trans)
print(log_transform_inverse(a_trans))
# Establish train, validation, and test sets
split_seed = 40   # ESTABLISH SPLIT SEED
# Split data into train/val/test sets
samples_train, samples_val, samples_test, cell_lines_test, cell_lines_val = Dataset.samples_train_test_split(
                                                                        response_df,
                                                                        num_val_cell_lines,
                                                                        num_test_cell_lines,
                                                                        split_seed,
                                                                        shuffle=True)
# Merge training and validation samples
samples_train_and_val = pd.concat([samples_train, samples_val], axis=0)
                                                                       
# Normalize the data
# Cell line data
cols_subset = [col for col in list(cell_line_data_original_df) if col.endswith("_exp")]
rows_subset = [x for x in cell_line_data_original_df.index if x not in cell_lines_test + cell_lines_val]

cell_line_data_df = Dataset.standardize_data(cell_line_data_original_df, cols_subset=cols_subset,
                                            rows_subset=rows_subset)
# Drug data
drug_data_df = Dataset.standardize_data(drug_data_original_df)<jupyter_output><empty_output><jupyter_text>Instantiate the new system<jupyter_code># Setup seed for PyTorch models
torch.manual_seed(11)
state_dict = None
def establish_new_system(drug_dim, cell_line_dim, code_dim, architecture_type,
                    num_layers, drug_hidden_dim1, drug_hidden_dim2, drug_hidden_dim3,
                    cell_line_hidden_dim1, cell_line_hidden_dim2, cell_line_hidden_dim3,
                     dropout, dropout_rate, code_activation, activation_func, out_activation,
                        batch_norm):
    # Establish specs
    specs = {"drug_dim": drug_dim,
            "cell_line_dim": cell_line_dim,
            "code_dim": code_dim,
            "architecture_type": architecture_type,
            "num_layers": num_layers,
            "drug_hidden_dim1": drug_hidden_dim1,
            "drug_hidden_dim2": drug_hidden_dim2,
            "drug_hidden_dim3": drug_hidden_dim3,
            "cell_line_hidden_dim1": cell_line_hidden_dim1,
            "cell_line_hidden_dim2": cell_line_hidden_dim2,
            "cell_line_hidden_dim3": cell_line_hidden_dim3,
            "dropout": dropout,
            "dropout_rate": dropout_rate,
            "code_activation": code_activation,
            "activation_func": activation_func,
            "out_activation": out_activation,
            "batch_norm": batch_norm}
    
    network = ModelAnalysis.instantiate_system(specs, state_dict)
    return network

establish_new_system_widget = interactive(establish_new_system,
                                     drug_dim=[294],
                                     cell_line_dim=[241, 1144],
                                     code_dim=widgets.IntSlider(value=10, min=2, max=30),
                                     architecture_type=["autoencoders", "linear"],
                                     num_layers=[1, 2, 3, None],
                                     drug_hidden_dim1=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     drug_hidden_dim2=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     drug_hidden_dim3=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     cell_line_hidden_dim1=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     cell_line_hidden_dim2=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     cell_line_hidden_dim3=widgets.BoundedIntText(value=128, min=16, max=2048),
                                     dropout=[False, True],
                                     dropout_rate=widgets.FloatSlider(value=0.5, min=0, max=1, step=0.05),
                                     code_activation=[False, True],
                                     activation_func=[nn.ReLU, nn.Tanh, nn.SELU],
                                     out_activation=[torch.sigmoid, None],
                                     batch_norm=[False, True])
display(establish_new_system_widget)
# MUST EXECUTE ALSO
network = establish_new_system_widget.result
network
print(network.out_activation)
print(next(network.parameters()).sum())
# Training params
num_epochs = 200
batch_size = 512
no_batches = samples_train.shape[0] // batch_size + 1
# Establish the device
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
learning_rate = 0.00001
reconstruction_term_drug = 0.05
reconstruction_term_cl = 0.1
independence_term_drug = 0.2
independence_term_cl = 0.2
reg_lambda = 0.00

optimizer = torch.optim.Adam(network.parameters(), lr=learning_rate)
criterion = nn.MSELoss()

print(no_batches, device)<jupyter_output>81 cpu<jupyter_text>Train with early stopping<jupyter_code>early_stopping = EarlyStopping(patience=10, delta=0.000)
model = ModelWithAutoencoders("Test", network)
train_losses = []
drug_independence_train_losses = []
cl_independence_train_losses = []
train_rmses = []
train_corrs = []
train_median_corrs = []
train_drug_rec_corrs = []
train_cl_rec_corrs = []
val_rmses = []
val_corrs = []
val_median_corrs = []
val_drug_rec_corrs = []
val_cl_rec_corrs = []

for epoch in range(1, num_epochs + 1):
    loss, drug_recounstruction_loss, cl_reconstruction_loss, drug_independence_loss, cl_independence_loss = model.train_with_independence_penalty(
    samples_train, cell_line_data_df, drug_data_df, batch_size, optimizer, criterion, reconstruction_term_drug,
    reconstruction_term_cl, independence_term_drug, independence_term_cl, reg_lambda)
    print("Epoch: {}, loss: {}, drug independence loss: {}, cell line independence loss: {}".format(
        epoch, loss, drug_independence_loss, cl_independence_loss))
    
    train_losses.append(loss.item())
    drug_independence_train_losses.append(drug_independence_loss.item())
    cl_independence_train_losses.append(cl_independence_loss.item())
    
    # Evaluate on training data
    predicted, y_true, drug_input, cl_input = model.predict(samples_train, cell_line_data_df, drug_data_df)
    preds, drug_reconstruction, cl_reconstruction = predicted
    t0 = time.time()
    try:
        train_rmse, train_corr = Model.evaluate_predictions(y_true, preds.cpu().numpy().reshape(-1))
    except ValueError:
        train_rmse, train_corr = np.nan, (np.nan, np.nan)
    # Drug reconstruction training error
    try:
        train_drug_rec_rmse, train_drug_rec_corr = Model.evaluate_predictions(drug_input.flatten(),
                                        drug_reconstruction.cpu().numpy().flatten())
    except ValueError:
        train_drug_rec_rmse, train_drug_rec_corr = np.nan, (np.nan, np.nan)

    # Cell line reconstruction training error
    try:
        train_cl_rec_rmse, train_cl_rec_corr = Model.evaluate_predictions(cl_input.flatten(),
                                        cl_reconstruction.cpu().numpy().flatten())
    except ValueError:
        train_cl_rec_rmse, train_cl_rec_corr = np.nan, (np.nan, np.nan)

    performance_df = Model.per_drug_performance_df(samples_train, preds.cpu())
    train_median_rmse = performance_df["Model RMSE"].median()
    train_median_corr = performance_df["Model correlation"].median()
    
    train_rmses.append(train_rmse)
    train_corrs.append(train_corr[0])
    train_median_corrs.append(train_median_corr)
    train_drug_rec_corrs.append(train_drug_rec_corr[0])
    train_cl_rec_corrs.append(train_cl_rec_corr[0])
    
    # Evaluate on validation data
    predicted, y_true, drug_input, cl_input = model.predict(samples_val, cell_line_data_df, drug_data_df)
    preds, drug_reconstruction, cl_reconstruction = predicted
    
    # Main validation error
    try:
        val_rmse, val_corr = Model.evaluate_predictions(y_true, preds.cpu().numpy().reshape(-1))
    except ValueError:
        val_rmse, val_corr = np.nan, (np.nan, np.nan)

    # Drug reconstruction validation error
    try:
        val_drug_rec_rmse, val_drug_rec_corr = Model.evaluate_predictions(drug_input.flatten(),
                                        drug_reconstruction.cpu().numpy().flatten())
    except ValueError:
        val_drug_rec_rmse, val_drug_rec_corr = np.nan, (np.nan, np.nan)

    # Cell line reconstruction validation error
    try:
        val_cl_rec_rmse, val_cl_rec_corr = Model.evaluate_predictions(cl_input.flatten(),
                                        cl_reconstruction.cpu().numpy().flatten())
    except ValueError:
        val_cl_rec_rmse, val_cl_rec_corr = np.nan, (np.nan, np.nan)
    # Per-drug validation evaluation
    performance_df = Model.per_drug_performance_df(samples_val, preds.cpu())
    val_median_corr = performance_df["Model correlation"].median()
        
    val_rmses.append(val_rmse)
    val_corrs.append(val_corr[0])
    val_drug_rec_corrs.append(val_drug_rec_corr[0])
    val_cl_rec_corrs.append(val_cl_rec_corr[0])
    val_median_corrs.append(val_median_corr)
    
    early_stopping(val_rmse, model.network)
        
    if early_stopping.early_stop:
        print("Early stopping")
        break
def learning_curve(metric1, metric2):
    """Plot learning curve related to a given metric."""
    dic = {"Training loss": train_losses,
          "Drug independence train loss": drug_independence_train_losses,
          "Cl independence train loss": cl_independence_train_losses,
           "Train RMSE": train_rmses,
           "Train corr.": train_corrs,
           "Train median corr.": train_median_corrs,
           "Train drug rec corr.": train_drug_rec_corrs,
           "Train cl rec corr.": train_cl_rec_corrs,
          "Val RMSE": val_rmses,
          "Val corr.": val_corrs,
          "Val median corr.": val_median_corrs,
          "Val drug rec corr.": val_drug_rec_corrs,
          "Val cl rec corr.": val_cl_rec_corrs}
    
    x = range(1, len(train_losses) + 1)
    print("Max {}: {}".format(metric1, max(dic[metric1])))
    sns.lineplot(x, dic[metric1], label=metric1)
    if metric2:
        print("Max {}: {}".format(metric2, max(dic[metric2])))
        sns.lineplot(x, dic[metric2], label=metric2)
    
    plt.show()
    
metrics = ["Training loss", "Drug independence train loss", "Cl independence train loss",
           "Train RMSE", "Train corr.", "Train median corr.", "Train drug rec corr.", "Train cl rec corr.",
          "Val RMSE", "Val corr.", "Val median corr.", "Val drug rec corr.", "Val cl rec corr."]

interact(learning_curve,
        metric1=metrics,
        metric2=[None] + metrics)<jupyter_output><empty_output><jupyter_text>Evaluate response predictions<jupyter_code># Evaluation on training set
samples = samples_train.copy()
print("Training data")
print()
# Make predictions
# Linear model case
if model.__class__ == Model:
    response_preds, y_true = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    print(response_preds.shape, y_true.shape)
# Model with autoencoders case
elif model.__class__ == ModelWithAutoencoders:
    predicted, y_true, drug_input, cl_input = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    response_preds, drug_rec_preds, cl_rec_preds = predicted
    print(response_preds.shape, y_true.shape, drug_input.shape, cl_input.shape)

# Create a dataframe with samples and corresponding predictions
samples_with_predictions_train_df = samples.copy()
samples_with_predictions_train_df.insert(3, "Predicted AUC", response_preds.reshape(-1))
print(samples_with_predictions_train_df.shape)

# Evaluate predictions numerically
rmse, corr = Model.evaluate_predictions(y_true, response_preds.numpy().reshape(-1))
print("Response prediction performance:")
print("RMSE: {}, Corr.: {}".format(rmse, corr))

per_drug_performance_train_df = Model.per_drug_performance_df(samples, response_preds)
print(per_drug_performance_train_df.shape)
per_drug_performance_train_df.insert(1, "Drug name", per_drug_performance_train_df["Drug ID"].map(
    full_dataset.gdsc_drug_id_to_gdsc_name_mapper))
print(per_drug_performance_train_df.shape)
per_drug_performance_train_df = per_drug_performance_train_df.sort_values("Model correlation", ascending=False)
print("Per drug median RMSE and correlation:", per_drug_performance_train_df["Model RMSE"].median(),
     per_drug_performance_train_df["Model correlation"].median())
print()


# Evaluate reconstruction predictions
# Metrics calculated on whole matrix
rmse, corr = Model.evaluate_predictions(drug_input.ravel(), drug_rec_preds.numpy().ravel())
print("Reconstruction performance:")
print("Whole drug data RMSE and correlation:", np.round(rmse, 3), corr)
rmse, corr = Model.evaluate_predictions(cl_input.ravel(), cl_rec_preds.numpy().ravel())
print("Whole cell line data RMSE and correlation:", np.round(rmse, 3), corr)

per_drug_performance_train_df.head(2)
# Evaluation on validation set
samples = samples_val.copy()
print("Validation data")
print()
# Make predictions
# Linear model case
if model.__class__ == Model:
    response_preds, y_true = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    print(response_preds.shape, y_true.shape)
# Model with autoencoders case
elif model.__class__ == ModelWithAutoencoders:
    predicted, y_true, drug_input, cl_input = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    response_preds, drug_rec_preds, cl_rec_preds = predicted
    print(response_preds.shape, y_true.shape, drug_input.shape, cl_input.shape)

# Create a dataframe with samples and corresponding predictions
samples_with_predictions_val_df = samples.copy()
samples_with_predictions_val_df.insert(3, "Predicted AUC", response_preds.reshape(-1))
print(samples_with_predictions_val_df.shape)

# Evaluate predictions numerically
rmse, corr = Model.evaluate_predictions(y_true, response_preds.numpy().reshape(-1))
print("Response prediction performance:")
print("RMSE: {}, Corr.: {}".format(rmse, corr))

per_drug_performance_val_df = Model.per_drug_performance_df(samples, response_preds)
print(per_drug_performance_val_df.shape)
per_drug_performance_val_df.insert(1, "Drug name", per_drug_performance_val_df["Drug ID"].map(
    full_dataset.gdsc_drug_id_to_gdsc_name_mapper))
print(per_drug_performance_val_df.shape)
per_drug_performance_val_df = per_drug_performance_val_df.sort_values("Model correlation", ascending=False)
print("Per drug median RMSE and correlation:", per_drug_performance_val_df["Model RMSE"].median(),
     per_drug_performance_val_df["Model correlation"].median())
print()


# Evaluate reconstruction predictions
# Metrics calculated on whole matrix
rmse, corr = Model.evaluate_predictions(drug_input.ravel(), drug_rec_preds.numpy().ravel())
print("Reconstruction performance:")
print("Whole drug data RMSE and correlation:", np.round(rmse, 3), corr)
rmse, corr = Model.evaluate_predictions(cl_input.ravel(), cl_rec_preds.numpy().ravel())
print("Whole cell line data RMSE and correlation:", np.round(rmse, 3), corr)

per_drug_performance_train_df.head(2)
# Evaluation on training set in the original response variable (AUC)
samples = samples_train.copy()
print("Training data")
print()
# Make predictions
# Linear model case
if model.__class__ == Model:
    response_preds, y_true = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    print(response_preds.shape, y_true.shape)
# Model with autoencoders case
elif model.__class__ == ModelWithAutoencoders:
    predicted, y_true, drug_input, cl_input = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    response_preds, drug_rec_preds, cl_rec_preds = predicted
    print(response_preds.shape, y_true.shape, drug_input.shape, cl_input.shape)
    
# Go back to orginal AUC
response_preds = torch.from_numpy(log_transform_inverse(response_preds.detach().numpy()))
y_true = log_transform_inverse(y_true)
samples["AUC"] = samples["AUC"].map(log_transform_inverse)

# Create a dataframe with samples and corresponding predictions
samples_with_predictions_train_df = samples.copy()
samples_with_predictions_train_df.insert(3, "Predicted AUC", response_preds.reshape(-1))
print(samples_with_predictions_train_df.shape)

# Evaluate predictions numerically
rmse, corr = Model.evaluate_predictions(y_true, response_preds.numpy().reshape(-1))
print("Response prediction performance:")
print("RMSE: {}, Corr.: {}".format(rmse, corr))

per_drug_performance_train_df = Model.per_drug_performance_df(samples, response_preds)
print(per_drug_performance_train_df.shape)
per_drug_performance_train_df.insert(1, "Drug name", per_drug_performance_train_df["Drug ID"].map(
    full_dataset.gdsc_drug_id_to_gdsc_name_mapper))
print(per_drug_performance_train_df.shape)
per_drug_performance_train_df = per_drug_performance_train_df.sort_values("Model correlation", ascending=False)
print("Per drug median RMSE and correlation:", per_drug_performance_train_df["Model RMSE"].median(),
     per_drug_performance_train_df["Model correlation"].median())
print()


# Evaluate reconstruction predictions
# Metrics calculated on whole matrix
rmse, corr = Model.evaluate_predictions(drug_input.ravel(), drug_rec_preds.numpy().ravel())
print("Reconstruction performance:")
print("Whole drug data RMSE and correlation:", np.round(rmse, 3), corr)
rmse, corr = Model.evaluate_predictions(cl_input.ravel(), cl_rec_preds.numpy().ravel())
print("Whole cell line data RMSE and correlation:", np.round(rmse, 3), corr)

per_drug_performance_train_df.head(2)
samples_with_predictions_train_df["Predicted AUC"].min()
# Evaluation on validation set in the original response variable (AUC)
samples = samples_val.copy()
print("Validation data")
print()
# Make predictions
# Linear model case
if model.__class__ == Model:
    response_preds, y_true = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    print(response_preds.shape, y_true.shape)
# Model with autoencoders case
elif model.__class__ == ModelWithAutoencoders:
    predicted, y_true, drug_input, cl_input = model.predict(samples, cell_line_data_df,
                                                           drug_data_df)
    response_preds, drug_rec_preds, cl_rec_preds = predicted
    print(response_preds.shape, y_true.shape, drug_input.shape, cl_input.shape)
    
# Go back to orginal AUC
response_preds = torch.from_numpy(log_transform_inverse(response_preds.detach().numpy()))
y_true = log_transform_inverse(y_true)
samples["AUC"] = samples["AUC"].map(log_transform_inverse)

# Create a dataframe with samples and corresponding predictions
samples_with_predictions_val_df = samples.copy()
samples_with_predictions_val_df.insert(3, "Predicted AUC", response_preds.reshape(-1))
print(samples_with_predictions_val_df.shape)

# Evaluate predictions numerically
rmse, corr = Model.evaluate_predictions(y_true, response_preds.reshape(-1))
print("Response prediction performance:")
print("RMSE: {}, Corr.: {}".format(rmse, corr))

per_drug_performance_val_df = Model.per_drug_performance_df(samples, response_preds)
print(per_drug_performance_val_df.shape)
per_drug_performance_val_df.insert(1, "Drug name", per_drug_performance_val_df["Drug ID"].map(
    full_dataset.gdsc_drug_id_to_gdsc_name_mapper))
print(per_drug_performance_val_df.shape)
per_drug_performance_val_df = per_drug_performance_val_df.sort_values("Model correlation", ascending=False)
print("Per drug median RMSE and correlation:", per_drug_performance_val_df["Model RMSE"].median(),
     per_drug_performance_val_df["Model correlation"].median())
print()


# Evaluate reconstruction predictions
# Metrics calculated on whole matrix
rmse, corr = Model.evaluate_predictions(drug_input.ravel(), drug_rec_preds.numpy().ravel())
print("Reconstruction performance:")
print("Whole drug data RMSE and correlation:", np.round(rmse, 3), corr)
rmse, corr = Model.evaluate_predictions(cl_input.ravel(), cl_rec_preds.numpy().ravel())
print("Whole cell line data RMSE and correlation:", np.round(rmse, 3), corr)

per_drug_performance_val_df.head(2)
per_drug_performance_val_df.head(20)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09212070669631826, 'fraction_numerical': 0.01116567828020756, 'mean_word_length': 3.5228050569253333, 'pattern_counts': {'":': 63, '<': 88, '<?xml version=': 0, '>': 87, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 165, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '274411', 'license': 'MIT', 'n_tokens_mistral': 20940, 'n_tokens_neox': 19567, 'n_words': 3885, 'path': 'Scripts/Model development/check_architectures_and_optimization.ipynb', 'repo_name': 'kkoras/rec-system-for-drug-response', 'size': '64752'}