<jupyter_start><jupyter_text>Double-feature Models<jupyter_code>import datetime
from pathlib import Path
import pickle
import sys
import time

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.compose import make_column_selector, make_column_transformer, ColumnTransformer
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import GridSearchCV

dir_generic_files = Path("/home/enrico/shared_virtualbox/phd_projects_Enrico_Gandini/phd_project_similarity_prediction/generic_input_files/")
sys.path.append(dir_generic_files.as_posix())
import machine_learning_helpers as mlh<jupyter_output><empty_output><jupyter_text>Define directories and files<jupyter_code>#The current notebook should be
#in the main directory containing queried results.
dir_results = Path.cwd()<jupyter_output><empty_output><jupyter_text>Select the date when the survey ended, and define the directory containing survey results up to that date.<jupyter_code>date_end_survey = datetime.date(year=2021, month=6, day=28)
dir_queries = Path(dir_results, f"queried_heroku_{date_end_survey}")<jupyter_output><empty_output><jupyter_text>Load DataFrame containing aggregated survey answer. The dataset was produced by the "retrieve answers" script, using SQLAlchemy.<jupyter_code>file_agg = Path(dir_queries, "aggregated_survey_answers.csv")
df_agg = pd.read_csv(file_agg, index_col='id_chosenPair')
df_agg
df_agg.columns<jupyter_output><empty_output><jupyter_text>Directory that will contain fitted models.<jupyter_code>dir_models = Path(dir_results, "models_2D_3D")
dir_models.mkdir(exist_ok=True)<jupyter_output><empty_output><jupyter_text>Directory that will contain analyses and visualization of models.<jupyter_code>dir_models_analysis = Path(dir_models, "analysis")
dir_models_analysis.mkdir(exist_ok=True)<jupyter_output><empty_output><jupyter_text>Define variablesContained in the input data, and necessary to create nice figures.<jupyter_code>colname_score_2d = "tanimoto_cdk_Extended"
colname_score_3d = "TanimotoCombo"

colname_dist = "pchembl_distance"

colname_target = "target_name"

colname_subset = "pair_type"

colname_pair = "id_surveyPair"

colname_n_ans = "n_answers"

colname_n_simil = "n_similar"

colname_frac_simil = "frac_similar"
categories_subset = ["dis2D,dis3D",
                     "dis2D,sim3D",
                     "sim2D,dis3D",
                     "sim2D,sim3D",
                     ]
targets_names = ["HERG",
                 "5HT2B",
                 "CYP2D6",
                 ]
nicename_score_2d = (colname_score_2d
                     .replace("_", " ")
                     .title()
                     .replace("Cdk", "CDK")
                     )
nicename_score_3d = colname_score_3d
nicenames_scores = {colname_score_2d: nicename_score_2d,
                    colname_score_3d: nicename_score_3d,
                    }

nicename_dist = "pChEMBL Distance"

nicename_target = "Target"

nicename_subset = "Pair Type"

nicename_similar = "Similar"

nicename_experience = "Academic Qualification"

nicename_n_ans = "Number of Answers"

nicename_ans_percent = "Answer Percentage"

nicename_simil_percent = "Similarity Percentage"

nicename_pair = "Pair ID"

nicename_n_pairs = "Number of Pairs"

nicename_percent_pairs_subset = "Pair Percentage in each subset"
FIGSIZE_GOLD = (9.556, 5.906) # Golden Rectangle / (in, in)
FIGSIZE_PAGE = (8.27, 11.69) #A4 page without margins.
FIGSIZE_SQUARE = (FIGSIZE_GOLD[0], FIGSIZE_GOLD[0])

FONTSIZE = 14

kwargs_fig_basic = {"constrained_layout": True, "figsize": FIGSIZE_GOLD}
n_pairs_each_subset = 25
lim_score_2d = (0, 1)
lim_score_3d = (0, 2)
lims_scores = {colname_score_2d: lim_score_2d,
               colname_score_3d: lim_score_3d,
               }

lim_percent = (0, 100)<jupyter_output><empty_output><jupyter_text>Define Training set<jupyter_code>colnames_features = [colname_score_2d, colname_score_3d]

X_train = df_agg.loc[:, colnames_features].copy()
n_train = X_train.shape[0]<jupyter_output><empty_output><jupyter_text>Training set `y` labels will also be added to the original DataFrame, for easier analyses.<jupyter_code>human_similarity = df_agg[colname_frac_simil] # fraction of human experts that considered a pair of molecules to be similar

thres_human = 0.5 # / fraction
y_train = (human_similarity >= thres_human).astype(int)
colname_judged = "judged_similar"
y_train.name = colname_judged
df_agg[colname_judged] = y_train<jupyter_output><empty_output><jupyter_text>Define all Logistic Regression models Define `kwargs` of Logistic Regression models with all penalty types.<jupyter_code>seed_rand = 1
max_iter = 2000
solver = "saga" # Supports all logistic regression penalties.


kwargs_base = {"random_state": seed_rand,
               "max_iter": max_iter,
               "solver": solver,
               }

kwargs_noreg = kwargs_base.copy()
kwargs_noreg["penalty"]: "none"

kwargs_l1 = kwargs_base.copy()
kwargs_l1["penalty"] = "l1"

kwargs_l2 = kwargs_base.copy()
kwargs_l2["penalty"] = "l2"

kwargs_enet = kwargs_base.copy()
kwargs_enet["penalty"] = "elasticnet"
#Elasticnet also requires specification of `l1_ratio` hyperparameter.
ratio_default = 0.5
kwargs_enet_default = kwargs_enet.copy()
kwargs_enet_default["l1_ratio"] = ratio_default<jupyter_output><empty_output><jupyter_text>Define Decision Tree and Random Forest Models<jupyter_code>kwargs_decisiontree = {"random_state": seed_rand,
                       "class_weight": "balanced",
                       }

kwargs_randomforest = kwargs_decisiontree.copy()<jupyter_output><empty_output><jupyter_text>Define models that will be trained with default hyperparameters<jupyter_code>initial_models_default = {"noreg": LogisticRegression(**kwargs_noreg),
                          "l1": LogisticRegression(**kwargs_l1),
                          "l2": LogisticRegression(**kwargs_l2),
                          "enet": LogisticRegression(**kwargs_enet_default),
                          "decisiontree": DecisionTreeClassifier(**kwargs_decisiontree),
                          "randomforest": RandomForestClassifier(**kwargs_randomforest),
                          }<jupyter_output><empty_output><jupyter_text>Fit Logistic Regression models with default hyperparameters<jupyter_code>fitted_models_default = {}

for name, model in initial_models_default.items():
    model.fit(X_train, y_train)
    
    fitted_models_default[name] = model
    print(f"fitted {name}")<jupyter_output>fitted noreg
fitted l1
fitted l2
fitted enet
fitted decisiontree
fitted randomforest<jupyter_text>Evaluate models based on default hyperparametersOn training-set.<jupyter_code>df_evals_default = [] #will contain results of all scores for default models
for name, model in fitted_models_default.items():
    pred = model.predict(X_train)
    df_agg[f"pred_default_{name}"] = pred
    
    correct = (y_train == pred).astype(int)
    n_correct = correct.sum()
    print(f"Model `{name}` correctly predicted {n_correct}/{n_train} pairs.")
    df_agg[f"correct_default_{name}"] = correct
    
    tmp = {"model": name,
           "n_correct": n_correct,
           }
    metrics_dict = mlh.various_metrics_binary_classification(model=model,
                                                             metrics=mlh.METRICS_BINARY_PROBA,
                                                             X=X_train,
                                                             y_true=y_train,
                                                             )
    tmp.update(metrics_dict)
    df_evals_default.append(tmp)


df_evals_default = pd.DataFrame(df_evals_default)
df_evals_default
file_evals_default = Path(dir_models_analysis, "metrics_double_feature_default_models.csv")
df_evals_default.to_csv(file_evals_default, index=False)<jupyter_output><empty_output><jupyter_text>Save default models<jupyter_code>dir_default = Path(dir_models, "double_feature_default_models")
dir_default.mkdir(exist_ok=True)
for name, model in fitted_models_default.items():
    file_model = Path(dir_default, f"{name}.pickle")
    with open(file_model, "wb") as f:
        pickle.dump(model, f)<jupyter_output><empty_output><jupyter_text>Define grids of hyperparameters, that will be later used for hyperparameter optimization.All hyperparameters of Logistic Regression are about penalty, so the model with no penalty does not have hyperparameters, and hyperparameter optimization is not needed.<jupyter_code>reg_strengths = np.geomspace(0.01, 100, num=20)
ratios = np.geomspace(0.01, 1, num=25)


grid_l1 = {"C": reg_strengths}

grid_l2 = grid_l1

grid_enet = {"C": reg_strengths,
             "l1_ratio": ratios,
             }


#Hyperparameters for Decision Tree and Random Forest
grid_decisiontree = {"min_samples_leaf": [1, 2, 5, 10],
                     "max_depth": [2, 3, 5, 10],
                     "min_samples_split": [2, 5, 10],
                     }
grid_randomforest = grid_decisiontree.copy()
grid_randomforest["n_estimators"] = [5, 10, 50, 100]


grids_hyperparams = {"l1": grid_l1,
                     "l2": grid_l2,
                     "enet": grid_enet,
                     "decisiontree": grid_decisiontree,
                     "randomforest": grid_randomforest,
                     }
kwargs_for_optim = {"l1": kwargs_l1,
                    "l2": kwargs_l2,
                    "enet": kwargs_enet,
                    "decisiontree": kwargs_decisiontree,
                    "randomforest": kwargs_randomforest,
                    }

n_folds_cv = 10

n_jobs = 5
initial_optimizers = {}
for name, grid in grids_hyperparams.items():
    kwargs_model = kwargs_for_optim[name]
    
    #Select estimator
    if name == "decisiontree":
        estimator = DecisionTreeClassifier
    elif name == "randomforest":
        estimator = RandomForestClassifier
    elif name in ["l1", "l2", "enet"]:
        estimator = LogisticRegression
    else:
        raise ValueError(f"`{name}` model not valid!")
    
    optimizer = GridSearchCV(estimator=estimator(**kwargs_model),
                             param_grid=grid,
                             cv=n_folds_cv,
                             n_jobs=n_jobs,
                             )
    initial_optimizers[name] = optimizer
fitted_optimizers = {}
for name, optimizer in initial_optimizers.items():
    start = time.perf_counter()
    optimizer.fit(X_train, y_train)
    delta_time = time.perf_counter() - start
    delta_time = datetime.timedelta(seconds=delta_time)
    print(f"Completed hyperparameter search of `{name}` "
          f"using {n_folds_cv}-fold CV in {delta_time} hours with {n_jobs} CPUs.")
    fitted_optimizers[name] = optimizer
df_evals_optim = [] #will contain results of all scores for all optim models
for name, optimizer in fitted_optimizers.items():
    pred = optimizer.predict(X_train)
    df_agg[f"pred_optim_{name}"] = pred
    
    correct = (y_train == pred).astype(int)
    n_correct = correct.sum()
    print(f"Optimized `{name}` correctly predicted {n_correct}/{n_train} pairs.")
    df_agg[f"correct_optimized_{name}"] = correct
    
    tmp = {"model": name,
           "n_correct": n_correct,
           }
    metrics_dict = mlh.various_metrics_binary_classification(model=optimizer,
                                                             metrics=mlh.METRICS_BINARY_PROBA,
                                                             X=X_train,
                                                             y_true=y_train,
                                                             )
    tmp.update(metrics_dict)
    df_evals_optim.append(tmp)


df_evals_optim = pd.DataFrame(df_evals_optim)
df_evals_optim
file_evals_optim = Path(dir_models_analysis, "metrics_double_feature_optimized_models.csv")
df_evals_optim.to_csv(file_evals_optim, index=False)<jupyter_output><empty_output><jupyter_text>Concatenate the two evaluation DataFrames for an easier visual inspection.<jupyter_code>#Before concatenation, add prefix to `model` column,
#to specify if models in concatenated DataFrame are default or optimized.
df_evals_default["model"] = "default_" + df_evals_default["model"]
df_evals_optim["model"] = "optimized_" + df_evals_optim["model"]


pd.concat([df_evals_default, df_evals_optim], ignore_index=True)<jupyter_output><empty_output><jupyter_text>***Notes***:* For default models, L1 is the best.* For optimized models, L1 is better for ROC AUC and Average Precision, whereas L2 is better for the other metrics.* For optimized models, Elastic Net became identical to L1 (I also checked parameters and hyperparameters) $\Rightarrow$ **do not use optimized Elastic Net**!* L2 and Elastic Net metrics on training set improved after optimization, whereas L1 metrics on training set slightly worsened after optimization. This is not an issue: the new parameters and hyperparameters after optimization should be more robust and less overfit to the training set, due to Cross Validation.  Save Optimized Models<jupyter_code>dir_optim = Path(dir_models, "double_feature_optimized_models")
dir_optim.mkdir(exist_ok=True)
for name, optimizer in fitted_optimizers.items():
    #Select final model with best hyperparameters,
    #refitted on whole dataset.
    model = optimizer.best_estimator_
    file_model = Path(dir_optim, f"{name}.pickle")
    with open(file_model, "wb") as f:
        pickle.dump(model, f)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09052740129640542, 'fraction_numerical': 0.011196228638774307, 'mean_word_length': 3.5377673796791442, 'pattern_counts': {'":': 35, '<': 72, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '27969', 'license': 'MIT', 'n_tokens_mistral': 4414, 'n_tokens_neox': 4212, 'n_words': 912, 'path': 'webapp/results_survey_molecular_similarity/.ipynb_checkpoints/create_double_feature_models_2D_3D_02-checkpoint.ipynb', 'repo_name': 'enricogandini/paper_similarity_prediction', 'size': '13576'}